I try to plot satellite ground track. I found lots of formulas but hardly find algorithm which describe how to uses them one-by-one and step-by-step. As I understand the algorithm will be:
1) Get orbital elements:
a - semimajor axis i - inclination Ω - right ascension of the ascending node ω - argument of perigee v - true anomaly Note: I need to plot ground track knowing only this elements and I don’t get them from TLE.
2) Get time variables:
tCur – current time tBegin – begin time of propagation period tEnd – end time of propagation period deltaT – time step
3) Calculate supporting variables like:
Ra = Ha + Re – radius of semimajor axis Rp = Hp + Re – radius of semiminor axis e = (Ra-Rp)/(Ra+Rp) – eccentricity T = 2PI*sqrt(a^3/MU) – period n = 2PI/T – mean motion Where: Ha – satellite apogee height Hp – satellite perigee height MU – Earth gravity constant = 398602
4) Calculate deltas of Ω and ω:
dΩ = -(1.5*sqrt(MU)*J2*(R3*R3)/((1-e^2)* a^3.5))*cos(i) dω = dΩ *(2.5*(sin(i)*sin(i)) - 2)/cos(i) Where: J2 – constant harmonic coefficient = 0.0010826257
5) Before loop calculate init values of:
E = 2*atan(tan(v/2)*sqrt((1-e)/(1+e))) – eccentricity anomaly M = E - e*sin(E) – mean anomaly I find another formulas for E, but I don’t know which is better for accuracy.
6) Begin loop and continue while tCur < tEnd
E = KeplerEquatation(M,e,1E-09) – calculate eccentricity anomaly by solving Kepler equitation M = E-e*sin(E) – change step in mean anomaly, also I find formula: M = n*(tCur-t0) – Where: t0 – time when satellite cross perigee point v = atan2(sqrt(1-e^2)*sin(E), cos(E)-e) – calculate true anomaly u = v + ω – latitude argument? Show angle position respective to center of Earth, but I don’t actually know what is it. r = a*(1-e*cos(E))
8) Calculate state vector movement coordinates XYZeci (ECI – Earth Centered Inertial) – position of satellite in 3D ?:
X = (cos(Ω)*cos(u)-sin(Ω)*sin(u)*cos(i)) * r Y = (sin(Ω)*cos(u)+cos(Ω)*sin(u)*cos(i)) * r Z = (sin(u)*sin(i)) * r
9) Rotate XYZ in ECI to ECEF(Earth-Centered Eartg-Fixed) with rotation
matrix: XYZecef = Rot(GST)* XYZeci Where: cos(GST) sin(GST) 0 Rot = -sin(GST) cos(GST) 0 0 0 1 GST – greenwich sidereal time (in radians)
10) Translate ECEF to LLA(longitude-latitude-altitude) to find coordinates of satellite for plot ground track on 2D map.
11) Change time step:
tCur = tCur + deltaT
12) Repeat from 6.
Is this algorithm correct for calculating satellite ground track and orbit position?