Orbit Propagation and Mission Geometry: Methods and Reference Data

Scope and conventions

The simulator takes an orbit defined by its classical elements at a chosen epoch and advances it forward in time under gravity. From the resulting trajectory it derives the quantities a mission analyst usually wants: the sub-satellite ground track, the elevation of the spacecraft seen from a ground station and the intervals during which it is visible, the evolution of the orbital elements, and whether the spacecraft is sunlit or in eclipse.

Lengths are in kilometres and time in seconds. Internally every angle is carried in radians and wrapped to the interval [0, 2π); results are presented in degrees. The gravitational parameter μ is expressed in km³/s². Position and velocity are resolved in an Earth-centred inertial frame and converted to other frames only where a particular result requires it, for example the rotating Earth-fixed frame used for the ground track.

Reference frames

Four frames appear in the calculations. Keeping them distinct is what makes the geometry unambiguous.

The Earth-centred inertial (ECI) frame is fixed with respect to the stars, with its x-axis toward the vernal equinox, z along the rotation axis, and y completing the right-handed set. The equations of motion are integrated here, because Newton's second law holds in an inertial frame.

The Earth-fixed (ECEF) frame shares its origin and z-axis with ECI but rotates with the planet. A position is taken from ECI to ECEF by a single rotation about the polar axis through the Greenwich sidereal angle. This frame carries longitude and latitude, and therefore the ground track.

The perifocal (PQW) frame lies in the orbit plane, with its first axis toward perigee and its third along the angular momentum. The position and velocity of a Keplerian orbit take their simplest form here, which is why the element-to-vector conversion passes through it.

The orbital or RTN frame travels with the spacecraft: a radial axis along the position vector, a normal axis along the orbital angular momentum, and a transverse axis completing the triad. Small control or perturbing accelerations are most naturally described in this frame, since radial, along-track, and cross-track effects act differently on the orbit.

Time and the sidereal angle

The epoch is supplied as a UTC calendar date and time and converted to a Julian Date JD. The orientation of the rotating Earth at that instant is given by the Greenwich Mean Sidereal Time, the angle between the Greenwich meridian and the inertial x-axis. It is evaluated from the standard polynomial in the number of days and Julian centuries since the J2000 epoch,

d = JD − 2 451 545.0,    T = d / 36 525

θ_GMST = 280.460 618 37° + 360.985 647 366°·d + 3.8793×10⁻⁴·T² − T³/3.871×10⁷

reduced to the range 0 to 360°. The leading rate, 360.9857° per day, is the rotation of the Earth in one solar day; it exceeds a full turn because the Earth must turn slightly past 360° to bring the Sun back to the meridian. This value, rather than a time read from an ephemeris kernel, fixes the sidereal angle at the start epoch, and the angle is then carried forward in time at the planet's rotation rate (next section).

UT1 and UTC are treated as equal. The difference is under a second and shifts the sidereal angle by well under an arcsecond, which is negligible for the geometry computed here.

Physical constants and rotation rates

Each central body is described by four numbers: its gravitational parameter μ, its equatorial radius R, the second zonal harmonic J₂ that measures the equatorial bulge, and the sidereal rotation rate ω that turns the body-fixed frame. The values used are the standard reference figures.

The Earth rate, 7.292 115×10⁻⁵ rad/s, corresponds to a sidereal day of about 86 164 s, some four minutes shorter than the 86 400 s solar day. Two further constants enter the solar geometry: the astronomical unit, taken as 149 597 870.7 km, and the length of the day, 86 400 s.

The radius and J₂ set the size and the leading distortion of the gravity field; μ sets the time scale of the motion through the orbital period; and ω sets how fast the ground beneath the orbit slides eastward.

Orbital elements and state vectors

An orbit is specified by six classical elements: the semi-major axis a and eccentricity e for the size and shape, the inclination i and the right ascension of the ascending node Ω for the orientation of the plane, the argument of perigee ω for the orientation of the ellipse within that plane, and the true anomaly ν for the position along the orbit. When the orbit is circular the perigee is undefined and the position is given instead by the argument of latitude u = ω + ν; when it is also equatorial the node is undefined as well, and the true longitude locates the spacecraft directly. The simulator resolves these degenerate cases automatically.

From elements to position and velocity

The position and velocity are first written in the perifocal frame, where with the semi-latus rectum p = a(1 − e²),

r_PQW = [ p·cosν / (1 + e·cosν),  p·sinν / (1 + e·cosν),  0 ]

v_PQW = √(μ/p) · [ −sinν,  e + cosν,  0 ]

and then rotated into the inertial frame by the 3–1–3 sequence through Ω, i, and ω. Writing c and s for cosine and sine, the rotation matrix from perifocal to inertial is

R = [ cΩcω − sΩsωci   −cΩsω − sΩcωci   sΩsi ;   sΩcω + cΩsωci   −sΩsω + cΩcωci   −cΩsi ;   sωsi   cωsi   ci ]

From position and velocity to elements

The inverse step recovers the elements from a state vector. With r and v the magnitudes of position and velocity, the specific angular momentum h = r×v, the node vector n = &zaccent;×h, and the eccentricity vector

e_vec = [ (v² − μ/r)·r_vec − (r_vec·v_vec)·v_vec ] / μ

the semi-major axis follows from the energy, 1/a = 2/r − v²/μ. The inclination is the angle of h from the pole; the node right ascension comes from n; the argument of perigee from the angle between n and e_vec; and the true anomaly from the angle between e_vec and r. Each angle is placed in the correct half-plane using the sign of an appropriate component, and for circular or equatorial orbits the argument of latitude and true longitude take over from the undefined quantities.

Anomalies and Kepler's equation

Three angles describe position along an ellipse. The true anomaly ν is the geometric angle from perigee. The eccentric anomaly E follows from

tan(ν/2) = √[(1 + e)/(1 − e)] · tan(E/2)

and the mean anomaly M, which advances uniformly in time at the mean motion, comes from Kepler's equation

M = E − e·sin E

Propagating in time is trivial for M but the geometry needs E, so Kepler's equation is inverted numerically. A Newton iteration,

E_k+1 = E_k − (E_k − e·sin E_k − M) / (1 − e·cos E_k)

started from E = M for moderate eccentricities and from E = π otherwise, converges to better than 10⁻¹² rad in a handful of steps across the full range of bound orbits.

Equations of motion

The state is the inertial position and velocity, and its time derivative is the velocity together with the acceleration produced by the modelled forces. The dominant term is the point-mass attraction of the central body,

a_2−body = −μ · r_vec / r³

To this is added the acceleration from the equatorial bulge. The J₂ term is the largest departure from a spherical field and, in inertial components with r the radius and z the polar coordinate,

a_J2 = (3/2)·J₂·μ·R²/r⁴ · [ (5z²/r² − 1)·x/r,  (5z²/r² − 1)·y/r,  (5z²/r² − 3)·z/r ]

The bracket shows why the bulge matters: the in-plane and out-of-plane components differ, so the force is not central, and over many orbits this is what turns the orbit plane and the line of apsides (quantified in the section on secular rates).

An optional constant acceleration may be applied in the orbital frame to represent low-thrust control. Given its radial, transverse, and normal components, the radial direction is taken along the position, the normal along the angular momentum r×v, and the transverse as their cross product; the three components are summed onto the gravitational acceleration. Building the frame from the instantaneous position and velocity keeps the thrust consistent throughout the orbit without reference to the slowly changing elements.

Numerical integration

The six-element state is advanced with an explicit Runge–Kutta method of the Dormand–Prince 8(7) family, which carries an eighth-order solution with a seventh-order error estimate and adapts its internal step to hold the error within tolerance. The relative tolerance is set to 10⁻⁹ and the absolute tolerance to 10⁻¹², and the solution is sampled onto the uniform output grid chosen by the user through the step size.

The whole interval is integrated in a single pass rather than restarted at every output point, which avoids the small discontinuities a restarting scheme introduces and is faster. As a check on the method, a pure two-body orbit conserves specific energy to about one part in 10⁹ over a day. Under J₂ the two-body energy is no longer constant: it oscillates as the bulge exchanges energy with the orbit, which is physical rather than numerical, while the true total energy including the bulge potential remains conserved.

Ground-track computation

The ground track is the path of the sub-satellite point, the spot on the surface directly beneath the spacecraft. At each instant the inertial position is rotated into the Earth-fixed frame through the current sidereal angle θ,

θ(t) = ( θ_GMST + ω·t ) mod 2π

x_ECEF =  x·cosθ + y·sinθ    y_ECEF = −x·sinθ + y·cosθ    z_ECEF = z

and longitude and latitude are read off directly,

λ = atan2(y_ECEF, x_ECEF)     φ = asin(z_ECEF / r) geocentric

Because the rotation in θ is what shifts each successive pass westward, the spacing between ground tracks is governed by the rotation rate ω relative to the orbital period. The latitude of the track is bounded by the inclination: a prograde orbit reaches a maximum geographic latitude equal to its inclination, which is why a 51.6° orbit covers the globe up to that parallel and no further.

Station visibility and access

A ground station at longitude λ, latitude φ, and altitude h has an inertial position that turns with the Earth,

r_sta = (R + h)·[ cos(θ+λ)·cosφ,  sin(θ+λ)·cosφ,  sinφ ]

The vector from the station to the spacecraft is ρ = r_sat − r_sta. The spacecraft sits above the local horizon by the elevation angle, measured from the horizontal plane up to the line of sight,

el = 90° − arccos( r̂_sta · ρ̂ )

where r̂_sta is the local vertical and ρ̂ the unit line of sight. The spacecraft is in view whenever the elevation exceeds the station's minimum, a threshold usually set near 5° to 10° to clear terrain and horizon refraction. Each unbroken interval above that threshold is one access window, or pass. The simulator reports, for every pass, the time of acquisition when the elevation first rises through the minimum, the time of loss when it falls back below, the duration between them, and the peak elevation reached in between, which indicates how favourable the pass is.

Solar position and eclipse

Illumination requires the direction to the Sun. A low-precision analytic ephemeris, accurate to better than a hundredth of a degree, gives the Sun's inertial position from the same Julian-date time argument. With T the Julian centuries since J2000, the mean longitude and mean anomaly of the Sun are

λ_M = 280.460° + 36 000.771°·T     M_☉ = 357.5288° + 35 999.050°·T

the ecliptic longitude and Sun–Earth distance follow from the equation of centre,

λ_ecl = λ_M + 1.9147°·sin M_☉ + 0.0200°·sin 2M_☉

r_☉ = ( 1.000 141 − 0.016 709·cos M_☉ − 0.000 140·cos 2M_☉ )·AU

and, rotating from the ecliptic to the equator through the obliquity ε = 23.4393° − 0.0130°·T, the inertial unit vector to the Sun is

ŝ = [ cosλ_ecl,  cosε·sinλ_ecl,  sinε·sinλ_ecl ]

The spacecraft is judged to be in shadow with a cylindrical model. Let β be the angle at the Earth's centre between the spacecraft and the Sun, and define the limb angle from the spacecraft's radius,

γ = 90° + arccos( R / r )     β = arccos( r̂_sat · ŝ )

The spacecraft is sunlit while β ≤ γ and in eclipse once β exceeds γ, that is, once it has moved far enough around the night side that the Earth's disc cuts off the Sun. The fraction of each orbit spent in eclipse, reported with the run, drives the sizing of solar arrays and batteries and the thermal cycling of the spacecraft.

The model uses a spherical Earth and a single umbral boundary; it omits the partial shadow of the penumbra and any atmospheric effect at the limb.

Secular rates from J₂

Although the simulator obtains the element history by integration, it is worth knowing the average rates the bulge imposes, because they explain the slow drift seen in the plots and underpin several orbit designs. With the mean motion n = √(μ/a³) and the semi-latus rectum p = a(1 − e²), the node and perigee drift on average as

Ω̇ = −(3/2)·n·J₂·(R/p)²·cos i nodal regression

ω̇ = (3/4)·n·J₂·(R/p)²·(4 − 5·sin²i) apsidal rotation

The node rate is negative for a prograde orbit, so the plane drifts westward; for a typical low orbit it amounts to a few degrees per day. Two consequences are widely used. Setting the node rate equal to the Earth's mean motion around the Sun, about 0.9856° per day, makes the orbit plane keep a fixed angle to the Sun, the condition for a Sun-synchronous orbit, which for low altitudes calls for a slightly retrograde inclination near 98°. Setting the bracket 4 − 5·sin²i to zero freezes the line of apsides, which occurs at the critical inclinations of 63.4° and 116.6° used by highly eccentric communications and early-warning orbits.

Assumptions and references

The force model presently includes the point-mass attraction and the J₂ zonal term, with an optional constant orbital-frame thrust. Higher zonal and tesseral harmonics, atmospheric drag, third-body attraction from the Sun and Moon, and solar radiation pressure are not yet modelled, so results are most accurate over spans of days to weeks in low and medium orbits where J₂ dominates. The surface geometry uses a spherical body and geocentric latitude, the solar position is the low-precision series above, and the shadow is a single-boundary cylinder. Time is referred to UTC with UT1 taken equal to it.

The methods follow standard references, to which the reader is referred for derivations and for the higher-fidelity models that extend this work.

• Vallado, D. A. Fundamentals of Astrodynamics and Applications. Microcosm Press.
• Bate, R. R., Mueller, D. D., and White, J. E. Fundamentals of Astrodynamics. Dover.
• Montenbruck, O., and Gill, E. Satellite Orbits: Models, Methods and Applications. Springer.

Reference notes for the orbit-propagation simulator. Constants and equations above reflect the models actually evaluated by the tool.