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.
Space Exploration Stack Exchange is a question and answer site for spacecraft operators, scientists, engineers, and enthusiasts. It only takes a minute to sign up.
Sign up to join this communityI 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.
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.
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.
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()