coordinate system - Identify the latitudinal and longitudinal point directly beneath an object in space

Given the azimuth and altitude / latitude and longitude from the observer to the satellite, how can I accurately find the latitude and longitude of the point of intersection directly under the satellite?

This problem has been plaguing me for quite some time now, thank you in advance for your help!

Edit: I also have the L1/2 p and L1/2 c/a values.

Answer

Assuming you have the distance from the observer to the satellite--without which the problem has no definite solution--then this amounts to solving three subproblems, using the strategy of computing the satellite's geocentric Cartesian (x,y,z) coordinates.

In the following, a is the semimajor axis (6,378,137.0 meters in WGS 84) and b is the semiminor axis (approximately 6,356,752.314 245 meters). All calculations should consistently use either radians or degrees--whatever is preferred for trigonometric functions by your software. Double-precision calculations will give sub-micron accuracy, assuming such accuracy is present in the original data!

The (working) examples use Mathematica.

Convert between geodetic and Cartesian coordinates

We may parameterize a vertical cross-section of the ellipsoid in the form (x, z) = (a * cos(t), b * sin(t)) for the latitude t, and then rotate it around the z axis. Inverting the equations enables us to recover (lon, lat) from (x,y,z).

ellipsoidToCartesian[{lon_, lat_}, {a_,b_}] := {a Cos[lat] Cos[lon], a Cos[lat] Sin[lon], b Sin[lat]};

cartesianToEllipsoid[{x_, y_, z_}, {a_,b_}] := {ArcTan[x, y], ArcTan[Norm[{x, y}]/a, z/b]};

Obtain a basis for the tangent plane at any point

This code returns a list whose first element is the origin and whose second element is the orthonormalized basis {east, north, up}. It computes the derivatives of the geodetic parameters and normalizes them to obtain a basis for the tangent plane, then uses their cross product to calculate the unit "up" vector.

tangent[{lon_, lat_}, {a_, b_}] := Module[{north, east, up},
  north = {-a Cos[lon] Sin[lat], -a Sin[lat] Sin[lon], b Cos[lat]};
  north = north / Norm[north];
  east = {-a Cos[lat] Sin[lon], a Cos[lat] Cos[lon], 0};
  east = east / Norm[east];
  up = Cross[east, north];

  {ellipsoidToCartesian[{lon, lat}, {a, b}], {east, north, up}}
  ]

Compute the satellite's Cartesian coordinates

The altitude and azimuth, along with the distance from the observer, are essentially local spherical coordinates. Converting them into a vector displacement gives the satellite's coordinates relative to the observer's location.

The altitude and azimuth are relative to a local coordinate system based on (east, north, up). The following presumes the azimuth is given east of north.

satellite[{azimuth_, altitude_}, {lon_, lat_}, elevation_, {a_, b_}] := 
  Module[{origin, east, north, up},
   {origin, {east, north, up}} = tangent[{lon, lat}, {a, b}];
   elevation (east Sin[azimuth] Cos[altitude] + 

       north Cos[azimuth] Cos[altitude] + up Sin[altitude]) + origin
   ];

The solution

Having obtained the satellite's Cartesian coordinates, we are done, because all that is needed is to convert them back to the geodetic coordinates.

position[{azimuth_, altitude_}, {lon_, lat_}, elevation_, {a_, b_}] := 
 Module[{xyz},
  xyz = satellite[{azimuth, altitude}, {lon, lat}, elevation, {a, b}];
  cartesianToEllipsoid[xyz, {a, b}]
  ]

A dynamic visualization confirms all this works as intended. A line segment drawn from the origin (gray) to the satellite (blue) indeed coincides with the computed satellite coordinates (red).

For the convenience of those who may have access to Mathematica, code for the visualization follows.

Module[{world, azimuth, altitude, lon, lat, elevation, a = 1, s, p, o},
 Manipulate[
  world = 
   ParametricPlot3D[
    ellipsoidToCartesian[{lon, lat}, {1, b}], {lon, -\[Pi], \[Pi]}, {lat, -\[Pi]/2, \[Pi]/2}, 
    PlotStyle -> Opacity[0.5], MeshFunctions -> {#4 &, #5 &}, 

    MeshStyle -> Directive[Opacity[0.5], Gray], Boxed -> False];
  s = satellite[{azimuth, altitude}, {lon, lat}, elevation, {a, b}];
  p = ellipsoidToCartesian[cartesianToEllipsoid[s, {a, b}], {a, b}];
  o = ellipsoidToCartesian[{lon, lat}, {a, b}];
  Show[Graphics3D[{
     Thick, Line[{{0, 0, 0}, s}], Line[{o, s}], 
     PointSize[0.02], Gray, Point[{0, 0, 0}], Green, Point[o], Blue, Point[s], Red, Point[p]
     }, Boxed -> False], world],
  {{azimuth, 0, "Azimuth"}, 0, 2 \[Pi]},
  {{altitude, \[Pi]/2, "Altitude"}, 0, \[Pi]/2},

  {{lon, 0, "Longitude"}, -\[Pi], \[Pi]},
  {{lat, \[Pi]/3, "Latitude"}, -\[Pi]/2, \[Pi]/2},
  {{elevation, 1, "Distance"}, 0, 2},
  {{b, 1, "b/a"}, 0, 1}
  ]]