Ground Track

Introduction

The ground track of a satellite is the path that the satellite traces on the surface of the Earth as it moves in its orbit. It is important for mission planning, communication, and surveillance purposes. The ground track can be visualized by projecting the satellite’s position onto the Earth’s surface as it orbits.

Algorithm

Given:

Orbital elements: Semi-major axis \(a\), eccentricity \(e\), inclination \(i\), right ascension of the ascending node (RAAN) \(\Omega\), and argument of perigee \(\omega\).

The gravitational parameter of Earth \(\mu\) and the Earth’s rotation rate \(\omega_{\text{Earth}}\) are also required.

To find:

The latitude, \(\phi\) and longitude \(\lambda\) of the projection of satellite at various points in its orbit.

1. Compute the Orbital Period

The orbital period is determined using Kepler’s Third Law:

\[ T = 2\pi \sqrt{\frac{a^3}{\mu}} \]

This period sets the time span over which the satellite’s ground track will be computed.

2. Time Discretization

A time array is created to sample the satellite’s orbit. At each time step, the mean anomaly \(M\) is calculated as:

\[ M = 2\pi \frac{t}{T} \]

3. Solve Kepler’s Equation and Compute the True Anomaly

Kepler’s Equation:

\[ M = E - e \sin(E) \]

is solved iteratively using the Newton-Raphson method for the eccentric anomaly \(E\). For low eccentricities, the initial guess is set to \(M\); otherwise, the guess starts at \(\pi\).

True Anomaly:
Once \(E\) is known, the true anomaly \(\nu\) is computed using: \[ \nu = 2 \arctan \left(\sqrt{\frac{1+e}{1-e}} \tan \left(\frac{E}{2}\right)\right) \]

4. Compute the Radial Distance

The radial distance \(r\) from the Earth’s center is given by:

\[ r = a(1 - e\cos(E)) \]

5. Transforming to ECI Coordinates

1. Perifocal Coordinates:
   The satellite’s position in the orbital plane is expressed in the perifocal (PQW) coordinate system as:

\[ \mathbf{r}_{\text{pqw}} = \begin{bmatrix} r\cos(\nu) \\ r\sin(\nu) \\ 0 \end{bmatrix} \]

2. Coordinate Transformations:
   To convert the position from the PQW frame to the ECI frame, three successive rotations are applied:

• Rotate by the argument of perigee (\(\omega\)) about the z-axis.
• Rotate by the inclination (\(i\)) about the x-axis.
• Rotate by the RAAN (\(\Omega\)) about the z-axis.

The combined transformation is:

\[ \mathbf{r}_{\text{ECI}} = R_z(\Omega) \; R_x(i) \; R_z(\omega) \; \mathbf{r}_{\text{pqw}} \]

6. Convert from ECI to ECEF Coordinates

Greenwich Mean Sidereal Time (GMST):

A simplified linear model computes GMST as:

\[ \theta = \omega_{\text{Earth}} \times t \]

Rotation:

Rotate the ECI coordinates by the GMST angle \(\theta\) about the z-axis:

\[ \mathbf{r}_{\text{ECEF}} = R_z(\theta) \, \mathbf{r}_{\text{ECI}} \]

7. Calculate latitude and longitude

Project the Earth Centered Earth fixed (ECEF) position onto a 2D Earth surface:

Longitude (\(\lambda\)): \[ \lambda = \tan^{-1}\left(\frac{y}{x}\right) \]

Latitude (\(\phi\)): \[ \phi = \tan^{-1}\left(\frac{z}{\sqrt{x^2+y^2}}\right) \]

where \((x,y,z)\) are the components of the ECEF position vector.

Here is a python code to plot the ground track of a satellite:

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# Earth and gravitational parameters
R_EARTH = 6378.0                # Earth's radius in km
OMEGA_EARTH = 7.2921159e-5      # Earth's rotation rate in rad/s
MU_EARTH = 398600.0             # Earth's gravitational parameter in km^3/s^2
THETA_0 = np.radians(280.46061837) # Value of GMST at J2000 epoch 

def solve_kepler(M, e, tol=1e-8):
    E = M if e < 0.8 else np.pi
    diff = 1
    while np.abs(diff) > tol:
        diff = (E - e * np.sin(E) - M) / (1 - e * np.cos(E))
        E -= diff
    return E

def true_anomaly_from_eccentric(E, e):
    return 2 * np.arctan2(np.sqrt(1 + e) * np.sin(E / 2),
                          np.sqrt(1 - e) * np.cos(E / 2))

def compute_eci_position(a, e, i, Omega, omega, M):
    E = solve_kepler(M, e)
    nu = true_anomaly_from_eccentric(E, e)
    r = a * (1 - e * np.cos(E))
    r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0])

    R3_omega = np.array([[np.cos(omega), -np.sin(omega), 0],
                         [np.sin(omega),  np.cos(omega), 0],
                         [0,              0,             1]])
    R1_i = np.array([[1, 0, 0],
                     [0, np.cos(i), -np.sin(i)],
                     [0, np.sin(i),  np.cos(i)]])
    R3_Omega = np.array([[np.cos(Omega), -np.sin(Omega), 0],
                         [np.sin(Omega),  np.cos(Omega), 0],
                         [0,              0,             1]])

    r_eci = R3_Omega @ R1_i @ R3_omega @ r_pqw
    return r_eci

def gmst(t):
    return THETA_0 + OMEGA_EARTH * t

def eci_to_ecef(r_eci, t):
    theta = gmst(t)
    R = np.array([[ np.cos(theta),  np.sin(theta), 0],
                  [-np.sin(theta),  np.cos(theta), 0],
                  [0,               0,             1]])
    return R @ r_eci

def ecef_to_latlon(r_ecef):
    x, y, z = r_ecef
    lon = np.arctan2(y, x)
    lat = np.arctan2(z, np.sqrt(x**2 + y**2))
    return lat, lon

# Orbital elements
a = 7000.0                        
e = 0.0                         
i = np.radians(40.0)             
Omega = np.radians(30.0)         
omega = np.radians(45.0)        
nu0 = np.radians(0.0)            

# Compute mean anomaly from true anomaly
E0 = 2 * np.arctan2(np.tan(nu0/2) * np.sqrt((1 - e)/(1 + e)), 1)
M0 = E0 - e * np.sin(E0)

# Compute orbital period
period = 2 * np.pi * np.sqrt(a**3 / MU_EARTH)

# Time array
num_points = 5000
t_array = np.linspace(0, 2*period, num_points)

# Store ground track points
lats = []
lons = []

for t in t_array:
    M = M0 + 2 * np.pi * t / period
    r_eci = compute_eci_position(a, e, i, Omega, omega, M)
    r_ecef = eci_to_ecef(r_eci, t)
    lat, lon = ecef_to_latlon(r_ecef)
    lats.append(np.degrees(lat))
    lons.append(np.degrees(lon))

# Fix longitude discontinuities by inserting NaNs
lons = np.array(lons)
lats = np.array(lats)

diffs = np.abs(np.diff(lons))
jump_indices = np.where(diffs > 180)[0]

for idx in reversed(jump_indices):
    lons = np.insert(lons, idx + 1, np.nan)
    lats = np.insert(lats, idx + 1, np.nan)

# Load world map image
world_map = Image.open("latex_img/Equirectangular_projection_MAP.jpg")

# Plot
fig, ax = plt.subplots(figsize=(12, 6))
ax.imshow(world_map, extent=[-180, 180, -90, 90], aspect='auto')
ax.plot(lons, lats, 'r-', linewidth=1.5, label="Ground Track")

ax.set_xlim([-180, 180])
ax.set_ylim([-90, 90])
ax.set_xlabel("Longitude (degrees)")
ax.set_ylabel("Latitude (degrees)")
ax.set_title("Satellite Ground Track")
ax.grid(True)
ax.legend()

plt.tight_layout()
plt.show()