Calculating longitude, latitude and altitude using RAAN, argument of perigee etc

Question

I'm relatively new to orbital mechanics and I know there are a ton of resources and similar questions on the web, but I can't seem to get a straight answer/guidance from any of them.

I have a set of data regarding any particular satellite which consists of the following:

1. Mean motion
2. Mean anomaly
3. Argument of perigee
4. Eccentricity
5. Inclination
6. Right Ascension of Ascending Node (RAAN) AKA Longitude of Ascending Node

This data is extracted using PyEphem library. What I would like to calculate is the longitude, latitude and altitude of the satellite at any given time. I know what each parameter means but I don't know how they relate to my unknowns.

I tried reading papers, tried several simulations and checked out lots of articles about Keplerian elements conversion but I still can't figure out how to connect the dots. The main resources I have been studying are:

http://fgg-web.fgg.uni-lj.si/~/mkuhar/Pouk/SG/Seminar/Vrste_tirnic_um_Zemljinih_sat/Keplerian%20Elements%20Tutorial.html#epoch

https://au.mathworks.com/matlabcentral/fileexchange/54875-geostationary-satellites-tracking

https://au.mathworks.com/matlabcentral/fileexchange/982-gui-based-satellite-tracking-system

and I also read most of the related questions on this website. Would it be easier if I used real-time tracking data instead of using TLE to predict my unknown?

Thanks in advance!

Answer 1

1. Calculate eccentric anomaly. It can be done by solving Kepler's equation. Example Python code for Newton's method:

   M   = ...           # your mean anomaly
   ecc = ...           # your eccentricity
   maxIter  = 15       # maximum number of iterations
   maxError = 1e-11    # maximum error
   
   i = 0
   E = M if (ecc < 0.8) else math.pi
   F = E - ecc * math.sin(M) - M
   while ((abs(F) > maxError) and (i < maxIter)):
       E = E - F / (1 - ecc * math.cos(E))
       F = E - ecc * math.sin(E) - M
       i += 1
   
   return E
   

2. Use this paper to calculate the cartesian position from eccentric anomaly and other elements. You'll need to rotate some vectors, you can do that with rotation matrices described here.

3. Now you need to convert your cartesian position to spherical coordinates. The Pyhton expressions are as follows:

   r   = math.sqrt(x**2 + y**2 + z**2)
   lat = math.asin(z / r)
   lon = math.atan2(y, x)
   if lon < 0:
       lon += 2 * math.pi
   

Note that r here is the distance from the Earth barycenter. You should subtract Earth radius from it to calculate altitude. You can just take 6371 km as a mean value or use a more precise one from some ellipsoid models described here and here.

Also note that lon is relative to the equinox vector, not Earth's prime meridian, so you'll need to account for Earth rotation (get current Earth rotation angle about its axis relative to the equinox vector and subtract it from your derived longitude).

Answer 2

if you're using pyephem, you can get the latitude and longitude directly from the body object.

body = ephem.EarthSatellite() 
body._raan = ....  
body._ap = ....  

then get the lat, lon at datetime dt

body.compute(dt)  
body.sublat  
body.sublon