# 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.

-
Let me know if you need any more information! –  Mathematician Guy Sep 27 '12 at 17:31
If you want more than about 2.5 decimal digits of precision, then you need to model the earth as an oblate ellipsoid with inverse flattening of approximately 298.3. –  whuber Sep 27 '12 at 20:22
Cross-posted at math.stackexchange.com/questions/203558. Are you "Mathematician Guy"? If so, you should flag the math moderators to merge your accounts. –  whuber Sep 28 '12 at 21:11
What do "p1/2 p and p1/2 c/a values" mean? I hope they include information about the distance to the satellite, for otherwise this problem is not sufficiently well defined (the answer depends on the satellite's height). –  whuber Sep 28 '12 at 21:14
Meant to say L1 / L1 p and c/a. These do carry information that will allow us to get the distance from the observer to the satellite. –  Korvin Szanto Sep 29 '12 at 2:01

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}
]]

-
I'd try to translate this to cartesian coordinates: The observer is at $(R \cos \phi \cos \lambda, R\cos\phi R \sin \lambda, R \sin \phi)$ with $R=$ Earth radius, $\lambda=$ longitude, $\phi=$ latitude. Next, azimuth and altitude (measured against the tangent at the observer) determine a line and with the elevation you get a quadratic with two solutions (one is "behind" the observer in the wrong direction), thus the position $(x,y,z)$ of the satellite in space. Finally, $\frac z{H+R}$ is the sine of the latitude and $\frac y x$ is the tangent of the longitude.