The orbit of this satellite plotted using 2 body (Keplerian) dynamics is shown below:
Figure 7.3: Geosynchronous satellite orbit
Geostationary Orbit (GEO)
Altitude ~35,786 km
A geostationary orbit is a special case of geosynchronous orbit where the orbit is circular and lies exactly in the Earth’s equatorial plane, so the satellite appears fixed in the sky.
One possible exapmple of state of a satellite in geostationary orbit is given below:
The orbit of this satellite plotted using 2 body (Keplerian) dynamics is shown below:
Figure 7.4: Geostationary satellite orbit
Highly Elliptical Orbit (HEO) – Molniya Orbit
Molniya (“lightning”) was a military communications satellite system used by the Soviet Union. The satellites used highly eccentric elliptical orbits of approximately 63.4 degrees inclination and orbital period of about 12 hours (called Molniya orbits), which allowed them to be visible to near-polar regions for long periods.
One example of a Molniya satellite is, Molniya 3-50 (NORAD CAT ID 25847).
The approximate orbital parameters for this satellite at some time are: (From space-track.org)
The orbit of this satellite plotted using 2 body (Keplerian) dynamics is shown below:
Figure 7.5: HEO Satellite: Molniya 3-50 orbit
Codes used
The Python code used for plotting the orbits is given below: (With example of LEO satellite)
Code
# By Suyog Soman (MTech, AFM 2023-2025) using GPT-4oimport numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3D # necessary for 3d plottingdef class2cart(a, e, i, RAAN, omega, theta):""" Converts classical orbital elements to Cartesian coordinates (ECI frame). Parameters: a : float Semi-major axis [km] e : float Eccentricity i : float Inclination [rad] RAAN : float Right ascension of ascending node [rad] omega : float Argument of periapsis [rad] theta : float True anomaly [rad] Returns: r : numpy array (3,) Position vector in km v : numpy array (3,) Velocity vector in km/s """ mu =398600# Earth's gravitational parameter [km^3/s^2]# Distance from the central body in the perifocal frame r_pqw = (a * (1- e**2)) / (1+ e * np.cos(theta)) * np.array([np.cos(theta), np.sin(theta), 0])# Velocity in perifocal frame v_pqw = np.sqrt(mu / (a * (1- e**2))) * np.array([-np.sin(theta), e + np.cos(theta), 0])# Rotation matrices from perifocal to ECI frame: R3_W = np.array([ [np.cos(RAAN), -np.sin(RAAN), 0], [np.sin(RAAN), np.cos(RAAN), 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_w = np.array([ [np.cos(omega), -np.sin(omega), 0], [np.sin(omega), np.cos(omega), 0], [0, 0, 1] ])# Complete rotation matrix from perifocal to ECI frame R = R3_W @ R1_i @ R3_w# Transform position and velocity vectors to the ECI frame r = R @ r_pqw v = R @ v_pqwreturn r, vdef propagate_orbit(a, e, i, RAAN, omega, theta0, dt):""" Propagates the orbit using Kepler's classical method of analytical orbit propagation. Parameters: a : float Semi-major axis [km] e : float Eccentricity i : float Inclination [rad] RAAN : float Right Ascension of the Ascending Node [rad] omega : float Argument of periapsis [rad] theta0: float Initial true anomaly [rad] dt : float Propagation time [s] Returns: a, e, i, RAAN, omega, theta_final : tuple of floats Updated orbital elements (only theta changes for two-body propagation). """ mu =398600.0# Earth's gravitational parameter [km^3/s^2]# Mean motion (n) [rad/s] n = np.sqrt(mu/(a**3))# Compute the initial eccentric anomaly (E0) from theta0 E0 =2* np.arctan(np.sqrt((1- e)/(1+ e)) * np.tan(theta0/2))if E0 <0: E0 +=2*np.pi# Compute the initial mean anomaly (M0) M0 = E0 - e*np.sin(E0)# Compute the final mean anomaly (Mf) Mf = M0 + n*dt# Solve Kepler's equation: Ef - e*sin(Ef) = Mf using Newton-Raphson iteration Ef = Mf # initial guess tol =1e-8 error =1.0while np.abs(error) > tol: error = (Ef - e*np.sin(Ef) - Mf) / (1- e*np.cos(Ef)) Ef -= error# Compute the final true anomaly from Ef theta_final =2* np.arctan(np.sqrt((1+e)/(1-e)) * np.tan(Ef/2))if theta_final <0: theta_final +=2*np.pi# In a two-body problem, other orbital elements remain constant.return a, e, i, RAAN, omega, theta_finalif__name__=="__main__":# Define the initial classical orbital elements a =6718.827# Semi-major axis in km e =0.00019440# Eccentricity i = np.deg2rad(43.0000) # Inclination in radians RAAN = np.deg2rad(156.5261) # RAAN in radians omega = np.deg2rad(261.3203) # Argument of periapsis in radians theta0 = np.deg2rad(30.0) # Initial true anomaly in radians mu =398600.0# Earth's gravitational parameter [km^3/s^2]# Calculate orbital period: T = 2*pi * sqrt(a^3/mu) T =2* np.pi * np.sqrt(a**3/ mu)print(f"Orbital period (s): {T:.2f}")# Define propagation parameters: propagate over one orbital period dt =60# time step in seconds (1 minute) num_steps =int(np.ceil(T / dt))# Lists to store trajectory for plotting rx_list = [] ry_list = [] rz_list = []# Set the initial time and orbital state current_theta = theta0# Get initial state r, _ = class2cart(a, e, i, RAAN, omega, current_theta) rx_list.append(r[0]) ry_list.append(r[1]) rz_list.append(r[2])# Propagate the orbit step by stepfor step inrange(1, num_steps +1):# Propagate for dt seconds a, e, i, RAAN, omega, new_theta = propagate_orbit(a, e, i, RAAN, omega, current_theta, dt) current_theta = new_theta # update true anomaly for next step# Convert the updated orbital elements to Cartesian coordinates r, _ = class2cart(a, e, i, RAAN, omega, current_theta)# Save the state for plotting rx_list.append(r[0]) ry_list.append(r[1]) rz_list.append(r[2])# Plotting the 3D orbit with Earth at the centre fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d')# Plot the orbit trajectory ax.plot(rx_list, ry_list, rz_list, label='Orbit Trajectory', color='b')# Plot Earth as a sphere earth_radius =6378# Earth's radius in km u, v_grid = np.mgrid[0:2*np.pi:50j, 0:np.pi:25j] x_sphere = earth_radius * np.cos(u) * np.sin(v_grid) y_sphere = earth_radius * np.sin(u) * np.sin(v_grid) z_sphere = earth_radius * np.cos(v_grid) ax.plot_surface(x_sphere, y_sphere, z_sphere, color='c', alpha=0.5, linewidth=0)# Setting labels and title ax.set_xlabel('X (km)') ax.set_ylabel('Y (km)') ax.set_zlabel('Z (km)') ax.set_title('Orbit Propagation over One Orbital Period') ax.legend() ax.grid(True)# Set equal aspect ratio for all axes max_val = np.max(np.abs([rx_list, ry_list, rz_list])) ax.set_xlim([-max_val, max_val]) ax.set_ylim([-max_val, max_val]) ax.set_zlim([-max_val, max_val]) plt.show()
Orbital period (s): 5480.89
Plotting all in one fig.
Code
# By Suyog Soman (MTech, AFM 2023-2025) using GPT-4oimport numpy as npimport plotly.graph_objects as godef class2cart(a, e, i, RAAN, omega, theta):""" Converts classical orbital elements to Cartesian coordinates (ECI frame). """ mu =398600# Earth's gravitational parameter [km^3/s^2]# Position in perifocal frame r_pqw = (a * (1- e**2)) / (1+ e * np.cos(theta)) * np.array([np.cos(theta), np.sin(theta), 0])# Velocity in perifocal frame v_pqw = np.sqrt(mu / (a * (1- e**2))) * np.array([-np.sin(theta), e + np.cos(theta), 0])# Rotation matrices for perifocal to ECI conversion: R3_W = np.array([ [np.cos(RAAN), -np.sin(RAAN), 0], [np.sin(RAAN), np.cos(RAAN), 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_w = np.array([ [np.cos(omega), -np.sin(omega), 0], [np.sin(omega), np.cos(omega), 0], [0, 0, 1] ]) R = R3_W @ R1_i @ R3_w # overall rotation matrix r = R @ r_pqw v = R @ v_pqwreturn r, vdef propagate_orbit(a, e, i, RAAN, omega, theta0, dt):""" Propagates the orbit by dt seconds using Kepler's method. Only the true anomaly changes for a two-body orbit. """ mu =398600.0# Earth's gravitational parameter [km^3/s^2] n = np.sqrt(mu / (a**3))# Convert initial true anomaly to eccentric anomaly E0 =2* np.arctan(np.sqrt((1- e)/(1+ e)) * np.tan(theta0/2))if E0 <0: E0 +=2*np.pi M0 = E0 - e*np.sin(E0) Mf = M0 + n*dt# Solve Kepler's Equation using Newton-Raphson iteration Ef = Mf tol =1e-8 error =1.0while np.abs(error) > tol: error = (Ef - e*np.sin(Ef) - Mf) / (1- e*np.cos(Ef)) Ef -= error theta_final =2* np.arctan(np.sqrt((1+e)/(1-e)) * np.tan(Ef/2))if theta_final <0: theta_final +=2*np.pireturn a, e, i, RAAN, omega, theta_final# Define satellite orbits using a list of dictionaries.# Angles are given in degrees and will be converted to radians.satellites = [ {"name": "LEO: STARLINK-33612", "a": 6718.827, "e": 0.00019440, "i": np.deg2rad(43.0), "RAAN": np.deg2rad(156.5261), "omega": np.deg2rad(261.3203), "theta": np.deg2rad(30.0)}, {"name": "MEO: GPS NAVSTAR 76", "a": 26559.490, "e": 0.00874290, "i": np.deg2rad(55.2608), "RAAN": np.deg2rad(106.9860), "omega": np.deg2rad(242.0236), "theta": np.deg2rad(30.0)}, {"name": "Geosynchronous", "a": 42164.000, "e": 0.100000, "i": np.deg2rad(15.0), "RAAN": np.deg2rad(40.0), "omega": np.deg2rad(0.0), "theta": np.deg2rad(30.0)}, {"name": "Geostationary", "a": 42164.000, "e": 0.000000, "i": np.deg2rad(0.0), "RAAN": np.deg2rad(0.0), "omega": np.deg2rad(0.0), "theta": np.deg2rad(30.0)}, {"name": "HEO: Molniya 3-50", "a": 26553.584, "e": 0.68462900, "i": np.deg2rad(63.5126), "RAAN": np.deg2rad(201.1207), "omega": np.deg2rad(280.6099), "theta": np.deg2rad(30.0)}]mu =398600.0# km^3/s^2# Prepare a list to hold Plotly tracestraces = []colors = ['blue', 'green', 'red', 'magenta', 'black']for idx, sat inenumerate(satellites): a = sat["a"] e = sat["e"] i_val = sat["i"] RAAN = sat["RAAN"] omega = sat["omega"] theta0 = sat["theta"]# Calculate the orbital period for the satellite T =2* np.pi * np.sqrt(a**3/ mu) dt =60# time step in seconds num_steps =int(np.ceil(T / dt)) rx_list, ry_list, rz_list = [], [], [] current_theta = theta0# Save initial position r, _ = class2cart(a, e, i_val, RAAN, omega, current_theta) rx_list.append(r[0]) ry_list.append(r[1]) rz_list.append(r[2])# Propagate orbit over one periodfor step inrange(1, num_steps +1): a, e, i_val, RAAN, omega, new_theta = propagate_orbit(a, e, i_val, RAAN, omega, current_theta, dt) current_theta = new_theta r, _ = class2cart(a, e, i_val, RAAN, omega, current_theta) rx_list.append(r[0]) ry_list.append(r[1]) rz_list.append(r[2])# Create a trace for the satellite orbit trace = go.Scatter3d( x=rx_list, y=ry_list, z=rz_list, mode='lines', line=dict(color=colors[idx %len(colors)], width=2), name=sat["name"] ) traces.append(trace)# Create an Earth sphereearth_radius =6378# kmu = np.linspace(0, 2*np.pi, 50)v = np.linspace(0, np.pi, 25)u, v = np.meshgrid(u, v)x_sphere = earth_radius * np.cos(u) * np.sin(v)y_sphere = earth_radius * np.sin(u) * np.sin(v)z_sphere = earth_radius * np.cos(v)earth_trace = go.Surface( x=x_sphere, y=y_sphere, z=z_sphere, colorscale='Blues', opacity=0.5, showscale=False, name='Earth')# Create the figure with all tracesfig = go.Figure(data=traces + [earth_trace])fig.update_layout( title="Orbits of Satellites", scene=dict( xaxis_title='X (km)', yaxis_title='Y (km)', zaxis_title='Z (km)', aspectmode='data' ), legend=dict(x=0.02, y=0.98))fig.show()
