ConeFormulas2D Subroutine

  SUBROUTINE ConeFormulas2D( z_xx, z_xy, z_yy, nBdry, xs, tradial, xray, BotTop )
  !! analytic formula for the conical seamount

    REAL (KIND=8),     INTENT( IN  ) :: xs( 3 ), tradial( 2 ), xray( 2 )
    CHARACTER (LEN=3), INTENT( IN  ) :: BotTop                  ! Flag indicating bottom or top reflection
    REAL (KIND=8),     INTENT( OUT ) :: z_xx, z_xy, z_yy, nBdry( 2 )
    REAL (KIND=8) :: nBdry3dCone( 3 ), theta, phi, Rray, Radius, x, y   ! for cone reflection

  IF ( BotTop == 'BOT' ) THEN
     phi   = 15 * DegRad   ! 15 degree angle of seamount
     Rray  = xray( 1 )   ! cylindrical range of ray from source

     ! get bearing from cone to ray using (x,y) coordinate of the ray
     x = xs( 1 ) + Rray * tradial( 1 )
     y = xs( 2 ) + Rray * tradial( 2 )
     theta = atan2( y, x )

     nBdry3dCone =  [ -cos( theta ) * sin( phi ), -sin( theta ) * sin( phi ), cos( phi ) ]

     ! nBdry3dCone = -[ -cos( theta ) * tan( phi ), -sin( theta ) * tan( phi ), 1.0D0 ]
     ! nBdry3dCone = -[ ray2D( is )%x( 1 ), ray2D( is )%x( 2 ), -1.0D0  ]
     ! nBdry3dCone = NBdry3dCone / NORM2( NBdry3dCone )

     nBdry( 1 ) = DOT_PRODUCT( nBdry3dCone( 1 : 2 ), tradial )
     nBdry( 2 ) = nBdry3dCone( 3 )

     ! curvatures
     Radius = SQRT( x ** 2 + y ** 2 )   ! radius of seamount at the bounce point

     ! z = z_xx / 2 * x^2 + z_xy * xy + z_yy / 2 * y^2; coefs are the kappa matrix

     z_xx =   y * y / Radius**3 * tan( phi )
     z_xy = - x * y / Radius**3 * tan( phi )
     z_yy =   x * x / Radius**3 * tan( phi )
  END IF
END SUBROUTINE ConeFormulas2D