Orbit Determination

Introduction

Orbit determination is the process of establishing the orbit of an object using observational data. This data may come from radar, optical observations, or telemetry, and it is used to estimate the six Keplerian orbital elements (or an equivalent set of parameters) that fully describe the orbit.

The accuracy of orbit determination is crucial for satellite tracking, mission planning, collision avoidance, and scientific studies.

Orbit Determination Methods

There are several methods for orbit determination depending on type of data available and the desired accuracy. Some common methods for preliminary orbit determination include:

  1. From Optical Sightings: This method uses Right Ascension (RA) and Declination (Dec) data at three different times to determine the orbit. The RA and Dec are the angles that describe the position of an object in the sky relative to the celestial sphere.
RA and Dec

Figure 9.1: RA and Dec (By Tfr000 (talk) 15:34, 15 June 2012 (UTC)), CC BY-SA 3.0, Link

RA is the angle measured from the vernal equinox along the celestial equator, and it can be considered as astronomical equivalent of longitude; while Dec is the angle measured from the celestial equator to the object and it can be considered as the astronomical equivalent of latitude. This method is especially useful for celestial objects like stars, planets, asteroids and comets.

  1. Gibbs Method: This is a method that uses three position vectors (\(\mathbf{r_1}\), \(\mathbf{r_2}\), \(\mathbf{r_3}\)) to determine the orbit.

  2. Lambert’s Problem: This approach uses two position vectors and the time of flight between them to determine the orbit.

Lambert’s Problem

Given:

Initial position vector, \(\mathbf{r}_1\) with magnitude \(r_1 = \|\mathbf{r}_1\|\); Final position vector,\(\mathbf{r}_2\) with magnitude \(r_2 = \|\mathbf{r}_2\|\); Time of flight, \(\Delta t = t_2 - t_1\)

To find: \(\mathbf{v}_1\) and \(\mathbf{v}_2\), the velocity vectors at the initial and final positions, respectively; thus determining the orbit with its complete state vector. Then these can be converted to the six classical orbital elements.

Step 1. Compute the Geometry

  1. Transfer Angle, \(\Delta\theta\):
    Compute \[ \cos\Delta\theta = \frac{\mathbf{r}_1\cdot\mathbf{r}_2}{r_1\,r_2} \]

(Select the “short‐way” or “long‐way”.)

  1. Chord, \(c\), and Semi-perimeter, \(s\): \[ c = \sqrt{r_1^2 + r_2^2 - 2\,r_1\,r_2\cos\Delta\theta} \]

\[ s = \frac{r_1 + r_2 + c}{2} \]

  1. Parameter \(A\): \[ A = \sin\Delta\theta\, \sqrt{\frac{r_1\,r_2}{1-\cos\Delta\theta}} \]

Step 2. Solve for the Universal Anomaly

Introduce a universal variable \(x\) (and let \(z=x^2\)). Use the Stumpff functions to account for all conic sections.

  1. Stumpff Functions: \[ C(z)= \begin{cases} \dfrac{1-\cos\sqrt{z}}{z}, & z>0, \\ \dfrac{1}{2}, & z=0, \\ \dfrac{1-\cosh\sqrt{-z}}{z}, & z<0, \end{cases} \] \[ S(z)= \begin{cases} \dfrac{\sqrt{z}-\sin\sqrt{z}}{\sqrt{z^3}}, & z>0, \\ \dfrac{1}{6}, & z=0, \\ \dfrac{\sinh\sqrt{-z}-\sqrt{-z}}{\sqrt{(-z)^3}}, & z<0. \end{cases} \]

  2. Define \(y\): \[ y = r_1 + r_2 + A\Bigl(x^2\,C(z) - 1\Bigr),\quad \text{with } z=x^2. \]

  3. Time-of-Flight Equation: The time of flight is given by: \[ \Delta t = \frac{x^3\,S(z) + A\sqrt{y}}{\sqrt{\mu}}\,. \]

Define the function: \[ F(x) \equiv \frac{x^3\,S(z) + A\sqrt{y}}{\sqrt{\mu}} - \Delta t\,. \]

We now solve \(F(x) = 0\) for \(x\).

  1. Iterative Solution (Newton–Raphson): Choose an initial guess \(x_0\) and iterate: \[ x_{n+1} = x_n - \frac{F(x_n)}{F'(x_n)}\,, \]

where \(F'(x)\) is the derivative with respect to \(x\). Continue iterating until \(|F(x)|\) is below a chosen tolerance.

Step 3. Compute Lagrange Coefficients and Velocity Vectors

After determining \(x\) (and hence \(z\) and \(y\)), compute the Lagrange coefficients.

  1. Lagrange Coefficients: \[ f = 1 - \frac{y}{r_1}\,, \] \[ g = A\sqrt{\frac{y}{\mu}}\,. \]

Their time derivatives are: \[ \dot{f} = \frac{\sqrt{\mu}}{r_1\,r_2}\Bigl(x^2\,S(z)-1\Bigr) \] \[ \dot{g} = 1 - \frac{y}{r_2} \]

  1. Compute the Velocities: Using the Lagrange–Gauss formulas: \[ \mathbf{v}_1 = \frac{1}{g}\Bigl(\mathbf{r}_2 - f\,\mathbf{r}_1\Bigr) \] \[ \mathbf{v}_2 = \frac{1}{g}\Bigl(\dot{g}\,\mathbf{r}_2 - \mathbf{r}_1\Bigr) \]

Here, \(\mathbf{v}_1\) is the velocity at the initial point \(P_1\) (time \(t_1\)) and \(\mathbf{v}_2\) is the velocity at the final point \(P_2\) (time \(t_2\)).

Here is a Python implementation of the Lambert’s problem:

Code 
# By Suyog Soman (MTech, AFM 2023-2025) using GPT-4o

import numpy as np

def stumpff_C(z):
    if z > 1e-8:
        sqrt_z = np.sqrt(z)
        return (1 - np.cos(sqrt_z)) / z
    elif z < -1e-8:
        sqrt_neg_z = np.sqrt(-z)
        return (1 - np.cosh(sqrt_neg_z)) / z
    else:
        return 1/2

def stumpff_S(z):
    if z > 1e-8:
        sqrt_z = np.sqrt(z)
        return (sqrt_z - np.sin(sqrt_z)) / (sqrt_z**3)
    elif z < -1e-8:
        sqrt_neg_z = np.sqrt(-z)
        return (np.sinh(sqrt_neg_z) - sqrt_neg_z) / (sqrt_neg_z**3)
    else:
        return 1/6

def solve_lambert(r1_vec, r2_vec, dt, mu, tol=1e-11, max_iter=1000):
    # Magnitudes of r1 and r2
    r1 = np.linalg.norm(r1_vec)
    r2 = np.linalg.norm(r2_vec)
    
    # Compute the transfer angle (choose the short way)
    dot_r1r2 = np.dot(r1_vec, r2_vec)
    cos_dtheta = dot_r1r2 / (r1 * r2)
    # Clamp value to avoid numerical issues
    cos_dtheta = np.clip(cos_dtheta, -1, 1)
    dtheta = np.arccos(cos_dtheta)
    
    # If cross product is negative, use the long way 
    if np.cross(r1_vec, r2_vec)[-1] < 0:
        dtheta = 2*np.pi - dtheta

    # Compute chord c and semi-perimeter s
    c = np.sqrt(r1**2 + r2**2 - 2 * r1 * r2 * np.cos(dtheta))
    s = (r1 + r2 + c) / 2

    # Compute parameter A
    A = np.sin(dtheta) * np.sqrt(r1 * r2 / (1 - np.cos(dtheta)))
    
    # Define the function F(x) whose root gives the correct time-of-flight
    def F(x):
        z = x**2
        C = stumpff_C(z)
        S = stumpff_S(z)
        y = r1 + r2 + A * (x**2 * C - 1)
        # The time-of-flight equation
        return (x**3 * S + A * np.sqrt(y)) / np.sqrt(mu) - dt

    # Numerical derivative of F(x) using a small finite difference
    def dF_dx(x, dx=1e-6):
        return (F(x + dx) - F(x - dx)) / (2 * dx)
    
    # Initial guess for x 
    x = 0.1  # initial guess; can be adjusted
    
    # Newton-Raphson iteration
    for i in range(max_iter):
        Fx = F(x)
        if abs(Fx) < tol:
            break
        dF = dF_dx(x)
        if dF == 0:
            raise RuntimeError("Zero derivative encountered during iteration")
        x = x - Fx/dF
    else:
        raise RuntimeError("Newton-Raphson did not converge")
    
    # With converged x, compute z, C, S and y
    z = x**2
    C = stumpff_C(z)
    S = stumpff_S(z)
    y = r1 + r2 + A * (x**2 * C - 1)
    
    # Lagrange coefficients
    f = 1 - y / r1
    g = A * np.sqrt(y / mu)
    g_dot = 1 - y / r2

    # Compute velocity vectors using Lagrange-Gauss formulas
    v1 = (r2_vec - f * r1_vec) / g
    v2 = (g_dot * r2_vec - r1_vec) / g

    return v1, v2

# Example usage:
if __name__ == "__main__":
    # Gravitational parameter for Earth [km^3/s^2]
    mu = 398600

    # Example position vectors (in km)
    r1_vec = np.array([5000.0, 10000.0, 2100.0])
    r2_vec = np.array([-14000.0, 2500.0, 7000.0])
    
    # Time of flight in seconds 
    dt = 3600.0

    v1, v2 = solve_lambert(r1_vec, r2_vec, dt, mu)
    print("Initial velocity (v1):", v1)
    print("Final velocity (v2):", v2)
Initial velocity (v1): [-1.07364577  6.32500613  3.12687842]
Final velocity (v2): [ 3.17049053 -3.59197962 -2.86303442]

Reference:

  1. Curtis, H. D. (2020). “Orbital Mechanics for Engineering Students, Fourth Edition”, Butterworth-Heinemann.

  2. Bate, Mueller, White. (1971). “Fundamentals of Astrodynamics”.