Vis-Viva equation: It is a fundamental equation in astrodynamics that relates the speed of a satellite to its distance from the central body and the semi-major axis of its orbit. It is derived from the conservation of mechanical energy in an orbiting system. The total specific mechanical energy, \(\epsilon\), is the sum of the specific kinetic energy and specific potential energy:
\[
\epsilon = \frac{v^2}{2} - \frac{\mu}{r},
\]
where \(v\) is the orbital velocity, \(r\) is the distance from the central body, \(\mu\) is the standard gravitational parameter (\(\mu = G(M+m) \approx GM \approx 398600\ \text{km}^3/\text{s}^2\) for the Earth-Satellite system).
For any bound orbit, the total specific mechanical energy remains constant and can also be expressed in terms of the semi-major axis, \(a\):
Conic sections are the curves obtained by intersecting a plane with a double-napped cone. In orbital mechanics, they describe the possible paths of celestial bodies under the influence of gravity in the two-body problem. The four types of conic sections are:
Circle: A closed curve with eccentricity \(e = 0\); has negative specific orbital energy.
Ellipse: A closed curve with eccentricity \(0 < e < 1\) ; has negative specific orbital energy.
Parabola: An open curve with eccentricity \(e = 1\), representing a trajectory that just escapes the gravitational influence of the central body; has zero specific orbital energy.
Hyperbola: An open curve with eccentricity \(e > 1\), representing a trajectory that escapes the gravitational influence of the central body; has positive specific orbital energy.
Here is a python code that plots the orbits of different eccentricities with a common focus at the origin.
Code
# Converted from MATLAB code by Suyog Soman (MTech, AFM 2023-2025) with GPT-4oimport numpy as npimport matplotlib.pyplot as plt# List of eccentricity valuese_values = [0, 0.7, 1, 2.5]# Fixed periapsis distancer_p =1plt.figure(figsize=(10, 10))colors = plt.rcParams['axes.prop_cycle'].by_key()['color']for i, e inenumerate(e_values): col = colors[i %len(colors)]# p is determined from the periapsis distance: r_p = p/(1+e) => p = r_p*(1+e) p = r_p * (1+ e)if e <1:# Elliptical orbit: full 0 to 2pi is valid. theta = np.linspace(0, 2*np.pi, 500) r = p / (1+ e * np.cos(theta))elif np.isclose(e, 1):# Parabolic orbit: avoid theta = pi where r blows up. theta = np.linspace(-np.pi+0.1, np.pi-0.1, 500) r = p / (1+ e * np.cos(theta))else:# Hyperbolic orbit: valid for theta where denominator is positive theta_inf = np.arccos(-1/e) theta = np.linspace(-theta_inf, theta_inf, 500) r = p / (1+ e * np.cos(theta))# Convert polar to Cartesian coordinates x = r * np.cos(theta) y = r * np.sin(theta) plt.plot(x, y, label=f"e={e}", color=col)# Mark the common focus at the origin.plt.plot(0, 0, 'ko', label='Focus')plt.xlabel("x")plt.ylabel("y")plt.title("Orbits with Common Focus (rp = 1)")plt.legend(loc='upper right')plt.grid(True)plt.axis('equal')plt.xlim(-6, 6)plt.ylim(-6, 6)plt.show()
Classical Orbital Elements
Classical orbital elements are a set of six parameters that completely define the size, shape, orientation of an orbit, and the position of the satellite along the orbit. They are defined as:
Semi-major axis, \(a\): Defines the size of the orbit.
Eccentricity, \(e\): Describes the shape of the orbit (circular: \(e=0\), elliptical: \(0<e<1\), parabolic: \(e=1\), hyperbolic: \(e>1\)).
Inclination, \(i\): The tilt of the orbit with respect to the reference plane (usually the equatorial plane).
Right Ascension of the Ascending Node (RAAN), \(\Omega\): Specifies the horizontal orientation of the orbit by locating the ascending node.
Argument of Perigee, \(\omega\): Describes the orientation of the ellipse within the orbital plane.
True Anomaly, \(\nu\): Gives the satellite’s position along the orbit at a specific time.
The figure below illustrates the definitions of some of these these elements in the context of an elliptical orbit.
Figure 5.1: Classical Orbital Elements (By Lasunncty at the English Wikipedia, CC BY-SA 3.0, Link)
Eccentricity\(e\) describes how “stretched” the ellipse is. The inclination\(i\) is the angle between the orbital plane and the reference plane (the equatorial plane). The RAAN or the longitude of ascending node \(\Omega\) is the angle from a reference direction (the vernal equinox) to the ascending node, where the satellite crosses from south to north of the reference plane. The argument of perigee\(\omega\) is the angle from the ascending node to the point of closest approach (perigee) in the orbit. The true anomaly\(\nu\) is the angle between the direction of perigee and the current position of the satellite.
These elements form the basis for analyzing and predicting the motion of bodies in a two-body (Keplerian) system.
Orbit Equation:
The orbit equation, valid for all conic sections, relates the radial distance \(r\) to the true anomaly \(\theta\), the specific angular momentum \(h\), and the standard gravitational parameter \(\mu\): \[
\boxed{
r = \frac{h^2 / \mu}{1 + e \cos\theta}
}
\]
Converting Between ECI Coordinates and Classical Orbital Elements
The conversion between the Earth-Centered Inertial (ECI) frame (given by position vector \(\mathbf{r}\) and velocity vector \(\mathbf{v}\)) and classical orbital elements involves several steps:
From ECI to Classical Orbital Elements
Given: - r: Position vector in an inertial (ECI) frame
- v: Velocity vector in the same frame
- \(\mu\): Gravitational parameter
Follow these steps:
Compute the Magnitudes:\[
r = \|\mathbf{r}\|
\]
\[
v = \|\mathbf{v}\|
\]
Specific Angular Momentum:\[
\mathbf{h} = \mathbf{r} \times \mathbf{v}
\]
The conversion uses three successive rotations via the rotation matrices:
Rotation about the Z-axis by \(-\Omega\): \[
R_3(-\Omega) = \begin{bmatrix}
\cos\Omega & \sin\Omega & 0 \\
-\sin\Omega & \cos\Omega & 0 \\
0 & 0 & 1
\end{bmatrix}
\]
Rotation about the X-axis by \(-i\): \[
R_1(-i) = \begin{bmatrix}
1 & 0 & 0 \\
0 & \cos i & \sin i \\
0 & -\sin i & \cos i
\end{bmatrix}
\]
Rotation about the Z-axis by \(-\omega\): \[
R_3(-\omega) = \begin{bmatrix}
\cos\omega & \sin\omega & 0 \\
-\sin\omega & \cos\omega & 0 \\
0 & 0 & 1
\end{bmatrix}
\]
The overall rotation matrix is: \[
R = R_3(-\Omega)\,R_1(-i)\,R_3(-\omega)
\]
Then, the inertial position and velocity vectors are: \[
\mathbf{r} = R\,\mathbf{r}_{\text{pf}}
\]\[
\mathbf{v} = R\,\mathbf{v}_{\text{pf}}
\]
Here is a python code that converts the classical orbital elements to ECI coordinates.
Code
# Converted from MATLAB code by Suyog Soman (MTech, AFM 2023-2025) with GPT-4oimport numpy as npdef 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# Gravitational parameter of Earth [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 matrix from perifocal to ECI:# 1. Rotation about z-axis by RAAN R3_W = np.array([ [np.cos(RAAN), -np.sin(RAAN), 0], [np.sin(RAAN), np.cos(RAAN), 0], [0, 0, 1] ])# 2. Rotation about x-axis by inclination (i) R1_i = np.array([ [1, 0, 0], [0, np.cos(i), -np.sin(i)], [0, np.sin(i), np.cos(i)] ])# 3. Rotation about z-axis by argument of periapsis (omega) 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 # Using the "@" operator for matrix multiplication# Transform position and velocity vectors to the ECI frame r = R @ r_pqw v = R @ v_pqwreturn r, v# Example usage:if__name__=="__main__":# Provided orbital elements a =7199.239655216658# Semi-major axis [km] e =0.08294103697605933# Eccentricity i_deg =8.32282494084567# Inclination in degrees RAAN_deg =334.94055273017824# RAAN in degrees omega_deg =310.678594628741# Argument of periapsis in degrees theta_deg =78.72522050823235# True anomaly in degrees# Convert angles from degrees to radians i = np.deg2rad(i_deg) RAAN = np.deg2rad(RAAN_deg) omega = np.deg2rad(omega_deg) theta = np.deg2rad(theta_deg)# Compute the Cartesian coordinates r, v = class2cart(a, e, i, RAAN, omega, theta)print("Position vector (r) [km]:", r)print("Velocity vector (v) [km/s]:", v)
Again, if the resulting true anomaly is negative, it is adjusted by adding \(2\pi\).
4. Update the Orbital State
In the two-body problem, the other orbital elements remain unchanged during the propagation. The final orbital elements include the semi-major axis \(a\), eccentricity \(e\), inclination \(i\), right ascension of the ascending node (RAAN), argument of periapsis \(\omega\), and the computed true anomaly \(\theta_f\).
Code
# Converted from MATLAB code by Suyog Soman (MTech, AFM 2023-2025) with GPT-4oimport numpy as npdef 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 Time elapsed [s] Returns: a_final, e_final, i_final, RAAN_final, omega_final, theta_final : tuple of floats Final orbital elements. """ 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 = 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 = 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 = theta_final +2*np.pi# Other orbital elements remain unchanged in a two-body problem.return a, e, i, RAAN, omega, theta_final# Example usage:if__name__=="__main__":# Define initial orbital elements a =7200.0# Semi-major axis in km e =0.08# Eccentricity i = np.deg2rad(8.0) # Inclination in radians RAAN = np.deg2rad(335.0) # RAAN in radians omega = np.deg2rad(310.0) # Argument of periapsis in radians theta0 = np.deg2rad(80.0) # Initial true anomaly in radians dt =3600# Propagation time in seconds (e.g., 1 hour)# Propagate the orbit a_final, e_final, i_final, RAAN_final, omega_final, theta_final = propagate_orbit(a, e, i, RAAN, omega, theta0, dt)# Display final orbital elementsprint("Final Orbital Elements:")print("Semi-major axis (a):", a_final, "km")print("Eccentricity (e):", e_final)print("Inclination (i):", np.rad2deg(i_final), "degrees")print("RAAN:", np.rad2deg(RAAN_final), "degrees")print("Argument of Periapsis (omega):", np.rad2deg(omega_final), "degrees")print("True Anomaly (theta):", np.rad2deg(theta_final), "degrees")
Final Orbital Elements:
Semi-major axis (a): 7200.0 km
Eccentricity (e): 0.08
Inclination (i): 8.0 degrees
RAAN: 335.0 degrees
Argument of Periapsis (omega): 310.0 degrees
True Anomaly (theta): 275.15750711200366 degrees
Universal Variable Formulation
To extend the analytical solution to all conic sections—including circular, elliptical, parabolic, and hyperbolic trajectories — the universal variable formulation is used. This method introduces a single variable (often denoted by \(x\)) to solve the Kepler problem regardless of the orbit type.
Key Components
Stumpff Functions:
The Stumpff functions, \(C(z)\) and \(S(z)\), where \(z = \alpha x^2\) (with \(\alpha = 1/a\) for elliptical orbits and appropriately defined for parabolic/hyperbolic cases), help handle the various orbital regimes within one framework.
Universal Kepler’s Equation:\[
x = \sqrt{\mu}\,\Delta t + \left[1 - z\,S(z)\right] \frac{\mathbf{r}_0\cdot\mathbf{v}_0}{\sqrt{\mu}}\,x^2 + r_0\,C(z)\,x^3,
\]
where \(\Delta t\) is the time of flight and \(r_0 = |\mathbf{r}_0|\).
Algorithmic Process
Initial Setup:
Start with the initial state vectors \(\mathbf{r}_0\) and \(\mathbf{v}_0\), and compute the classical orbital elements.
Iterative Solution for\(x\):
Use an iterative method (such as Newton-Raphson) to solve the universal Kepler’s Equation for the variable \(x\).
Compute Lagrange Coefficients:
With \(x\) determined, compute the Lagrange coefficients \(f\), \(g\), \(\dot{f}\), and \(\dot{g}\) using: \[
\mathbf{r}(t) = f\,\mathbf{r}_0 + g\,\mathbf{v}_0,
\]
