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$
2
  • $\begingroup$ Maybe you should somehow emulate real TLE data for the library? $\endgroup$ – peterh Apr 2 '19 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 '19 at 6:58
1
$\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$
2
  • 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 '19 at 15:12
  • $\begingroup$ @astrojuanlu please feel free to add another answer, or to edit this one. Thanks! $\endgroup$ – uhoh Apr 2 '19 at 15:14
0
$\begingroup$

It is not obvious what orbital parameters you refer to, but both Skyfield and the underlying python implementation of SGP4 provide a way to generate a satellite object from some orbital parameters (equivalent to TLEs but not in the TLE format) and propagate an orbit for it. Here is some documentation:

In case the first reference goes offline, here's the example they provide to build a EarthSatellite object:

from skyfield.api import EarthSatellite, load
from sgp4.api import Satrec, WGS72

ts = load.timescale()
satrec = Satrec()
satrec.sgp4init(
    WGS72,           # gravity model
    'i',             # 'a' = old AFSPC mode, 'i' = improved mode
    5,               # satnum: Satellite number
    18441.785,       # epoch: days since 1949 December 31 00:00 UT
    2.8098e-05,      # bstar: drag coefficient (/earth radii)
    6.969196665e-13, # ndot: ballistic coefficient (revs/day)
    0.0,             # nddot: second derivative of mean motion (revs/day^3)
    0.1859667,       # ecco: eccentricity
    5.7904160274885, # argpo: argument of perigee (radians)
    0.5980929187319, # inclo: inclination (radians)
    0.3373093125574, # mo: mean anomaly (radians)
    0.0472294454407, # no_kozai: mean motion (radians/minute)
    6.0863854713832, # nodeo: right ascension of ascending node (radians)
)
sat = EarthSatellite.from_satrec(satrec, ts)

After the generation, the orbit can be propagated as described here in the Skyfield docs.

To understand the above notation, having a basic comprehension of orbital elements and TLEs is useful.

$\endgroup$

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.