3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Maybe you should somehow emulate real TLE data for the library? $\endgroup$ – peterh says reinstate Monica Apr 2 at 6:35
  • 1
    $\begingroup$ TLE is essentially Euler angles plus a couple non-essential parameters formatted into a rigidly defined text representation. If your orbital parameters are Euler angles, you should be able to turn them into TLE trivially. If you use a state vector (position+velocity), it gets a bit trickier. $\endgroup$ – SF. Apr 2 at 6:58
0
$\begingroup$

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

enter image description here

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()
$\endgroup$
  • 1
    $\begingroup$ Disclaimer: poliastro author Two minor details: matplotlib 3D capabilities are not enough to properly plot orbits (as can be seen in your image), so I recommend using Plotly instead (poliastro can do it). Also, OrbitalPy is unmaintained (last commit 2015) whereas poliastro is still active. $\endgroup$ – astrojuanlu Apr 2 at 15:12
  • $\begingroup$ @astrojuanlu please feel free to add another answer, or to edit this one. Thanks! $\endgroup$ – uhoh Apr 2 at 15:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.