Propagating satellites without TLE in python

Question

I know how to propagate satellites in python using ephem module (that uses TLE data). However, I am not sure how I would do the same for my own set of orbital parameters (i.e, without TLE data) in python.

Thanks.

Answer 1

note: Alas, PyEphem is deprecated, so to plot from TLE's use Skyfield, written/maintained by the same person as PyEphem. I show an example here.

---

Step 1: Convert your orbital parameters to a state vector

A state vector is a 6 dimensional vector which is the combination of a position vector $\mathbf{x}$ and a velocity vector $\mathbf{v}$. In a deterministic system like a simple two-body orbit calculation, it's all you need to determine mathematically the orbit and how the body will move at any given time afterward.

This answer explains how to do this and contains a python script.

There's also https://pythonhosted.org/OrbitalPy/examples/defining_orbits/defining_orbits/

and I think you can learn a lot from poliastro at http://docs.poliastro.space/en/latest/

Choose the method that you are most comfortable with.

Step 2: Roll your own orbit propagator:

I took the plotting from here and added a simple Python integration of an initial state vector from here or here or even here. I've only added the Earth's spherically symmetric force, if you want to add the effects due to oblateness for better accuracy, it's shown in two of those linked scripts. Look for the $J_2$ term.

The equation of motions are:

$$\dot{\mathbf{x}} = \mathbf{v}$$

$$\ddot{\mathbf{x}} = - \frac{\mathbf{GM}}{x^2} \mathbf{\hat{x}} = -GM \frac{\mathbf{x}}{x^3} $$

and those are scripted in the function deriv(X, t) which is independent of time. This is called by the integrator scipy.integrate.odeint. It uses variable step sizes internally, then re-interpolates to user-specified time points.

Output:

period 6307.12290204 seconds or 105.118715034 minutes
inclination 57.0 degrees
initial position [ 7378.137     0.        0.   ] km
initial velocity [    0.          4003.17019194  6164.34152276] m/s
initial speed 7350.13455625 m/s

def makecubelimits(axis, centers=None, hw=None):
    lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    if centers == None:
        centers = [0.5*sum(pair) for pair in lims] 

    if hw == None:
        widths  = [pair[1] - pair[0] for pair in lims]
        hw      = 0.5*max(widths)
        ax.set_xlim(centers[0]-hw, centers[0]+hw)
        ax.set_ylim(centers[1]-hw, centers[1]+hw)
        ax.set_zlim(centers[2]-hw, centers[2]+hw)
        print("hw was None so set to:", hw)
    else:
        try:
            hwx, hwy, hwz = hw
            print("ok hw requested: ", hwx, hwy, hwz)

            ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
            ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
            ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
        except:
            print("nope hw requested: ", hw)
            ax.set_xlim(centers[0]-hw, centers[0]+hw)
            ax.set_ylim(centers[1]-hw, centers[1]+hw)
            ax.set_zlim(centers[2]-hw, centers[2]+hw)

    return centers, hw

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -GMe * x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180

km  = 0.001
GMe = 3.986E+14 # m^3/s^2
Re  = 6378137.  # meters
alt = 1E+06     # meters
a   = Re + alt
T   = twopi * np.sqrt(a**3/GMe)

v0  = np.sqrt(GMe/a)
incdegs = 57.

inc = rads * incdegs
cinc, sinc = [f(inc) for f in (np.cos, np.sin)]

X0 = np.array([Re+alt, 0, 0] + [0, cinc*v0, sinc*v0])

print 'period {} seconds or {} minutes'.format(T, T/60.)
print 'inclination {} degrees'.format(incdegs)
print 'initial position {} km'.format(km * X0[:3])
print 'initial velocity {} m/s'.format(X0[3:])
print 'initial speed {} m/s'.format(v0)

times = np.linspace(0, T, 201)

answer, info = ODEint(deriv, X0, times, full_output=True)


theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = Re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
                     lon0[1]*cph + lon0[0]*sph,
                     lon0[2]) )
    lons.append(lon)

lat0 = Re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lat = Re*np.vstack((cth*cph, sth*cph, zth+sph))
    lats.append(lat)

if True:
    fig = plt.figure(figsize=[10, 8])  # [12, 10]

    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = answer.T[:3]
    ax.plot(km*x, km*y, km*z)
    for x, y, z in lons:
        ax.plot(km*x, km*y, km*z, '-k')
    for x, y, z in lats:
        ax.plot(km*x, km*y, km*z, '-k')

    centers, hw = makecubelimits(ax)

    print("centers are: ", centers)
    print("hw is:       ", hw)

    plt.show()