Cowell's formulation

For cases where we only study the gravitational forces, solving the Kepler's equation is enough to propagate the orbit forward in time. However, when we want to take perturbations that deviate from Keplerian forces into account, we need a more complex method to solve our initial value problem: one of them is Cowell's formulation.

In this formulation we write the two body differential equation separating the Keplerian and the perturbation accelerations:

$$\ddot{\mathbb{r}} = -\frac{\mu}{|\mathbb{r}|^3} \mathbb{r} + \mathbb{a}_d$$

For an in-depth exploration of this topic, still to be integrated in poliastro, check out https://github.com/Juanlu001/pfc-uc3m

An earlier version of this notebook allowed for more flexibility and interactivity, but was considerably more complex. Future versions of poliastro and plotly might bring back part of that functionality, depending on user feedback. You can still download the older version here.

First example

Let's setup a very simple example with constant acceleration to visualize the effects on the orbit.



In [1]:

    
# Temporary hack, see https://github.com/poliastro/poliastro/issues/281
from IPython.display import HTML
HTML('<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>')









    Out[1]:



In [2]:

    
import numpy as np
from astropy import units as u

from matplotlib import pyplot as plt
plt.ion()

from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.examples import iss

from poliastro.twobody.propagation import cowell
from poliastro.plotting import OrbitPlotter3D
from poliastro.util import norm

from plotly.offline import init_notebook_mode
init_notebook_mode(connected=True)

To provide an acceleration depending on an extra parameter, we can use closures like this one:



In [3]:

    
accel = 2e-5



In [4]:

    
def constant_accel_factory(accel):
    def constant_accel(t0, u, k):
        v = u[3:]
        norm_v = (v[0]**2 + v[1]**2 + v[2]**2)**.5
        return accel * v / norm_v

    return constant_accel



In [5]:

    
def custom_propagator(orbit, tof, rtol, accel=accel):
    # Workaround for https://github.com/poliastro/poliastro/issues/328
    if tof == 0:
        return orbit.r.to(u.km).value, orbit.v.to(u.km / u.s).value
    else:
        # Use our custom perturbation acceleration
        return cowell(orbit, tof, rtol, ad=constant_accel_factory(accel))



In [6]:

    
times = np.linspace(0, 10 * iss.period, 500)
times









    Out[6]:





$[0,~111.36212,~222.72424,~\dots,~55346.973,~55458.335,~55569.697] \; \mathrm{s}$



In [7]:

    
times, positions = iss.sample(times, method=custom_propagator)

And we plot the results:



In [8]:

    
frame = OrbitPlotter3D()

frame.set_attractor(Earth)
frame.plot_trajectory(positions, label="ISS")

frame.show()

Error checking



In [9]:

    
def state_to_vector(ss):
    r, v = ss.rv()
    x, y, z = r.to(u.km).value
    vx, vy, vz = v.to(u.km / u.s).value
    return np.array([x, y, z, vx, vy, vz])



In [10]:

    
k = Earth.k.to(u.km**3 / u.s**2).value



In [11]:

    
rtol = 1e-13
full_periods = 2



In [12]:

    
u0 = state_to_vector(iss)
tf = ((2 * full_periods + 1) * iss.period / 2).to(u.s).value

u0, tf









    Out[12]:





(array([ 8.59072560e+02, -4.13720368e+03,  5.29556871e+03,  7.37289205e+00,
         2.08223573e+00,  4.39999794e-01]), 13892.42425290754)



In [13]:

    
iss_f_kep = iss.propagate(tf * u.s, rtol=1e-18)



In [14]:

    
r, v = cowell(iss, tf, rtol=rtol)

iss_f_num = Orbit.from_vectors(Earth, r * u.km, v * u.km / u.s, iss.epoch + tf * u.s)



In [15]:

    
iss_f_num.r, iss_f_kep.r









    Out[15]:





(<Quantity [ -835.92108005,  4151.60692532, -5303.60427969] km>,
 <Quantity [ -835.92108005,  4151.60692532, -5303.60427969] km>)



In [16]:

    
assert np.allclose(iss_f_num.r, iss_f_kep.r, rtol=rtol, atol=1e-08 * u.km)
assert np.allclose(iss_f_num.v, iss_f_kep.v, rtol=rtol, atol=1e-08 * u.km / u.s)



In [17]:

    
assert np.allclose(iss_f_num.a, iss_f_kep.a, rtol=rtol, atol=1e-08 * u.km)
assert np.allclose(iss_f_num.ecc, iss_f_kep.ecc, rtol=rtol)
assert np.allclose(iss_f_num.inc, iss_f_kep.inc, rtol=rtol, atol=1e-08 * u.rad)
assert np.allclose(iss_f_num.raan, iss_f_kep.raan, rtol=rtol, atol=1e-08 * u.rad)
assert np.allclose(iss_f_num.argp, iss_f_kep.argp, rtol=rtol, atol=1e-08 * u.rad)
assert np.allclose(iss_f_num.nu, iss_f_kep.nu, rtol=rtol, atol=1e-08 * u.rad)

Numerical validation

According to [Edelbaum, 1961], a coplanar, semimajor axis change with tangent thrust is defined by:

$$\frac{\operatorname{d}\!a}{a_0} = 2 \frac{F}{m V_0}\operatorname{d}\!t, \qquad \frac{\Delta{V}}{V_0} = \frac{1}{2} \frac{\Delta{a}}{a_0}$$

So let's create a new circular orbit and perform the necessary checks, assuming constant mass and thrust (i.e. constant acceleration):



In [18]:

    
ss = Orbit.circular(Earth, 500 * u.km)
tof = 20 * ss.period

ad = constant_accel_factory(1e-7)

r, v = cowell(ss, tof.to(u.s).value, ad=ad)

ss_final = Orbit.from_vectors(Earth, r * u.km, v * u.km / u.s, ss.epoch + tof)



In [19]:

    
da_a0 = (ss_final.a - ss.a) / ss.a
da_a0









    Out[19]:





$2.989621 \times 10^{-6} \; \mathrm{\frac{km}{m}}$



In [20]:

    
dv_v0 = abs(norm(ss_final.v) - norm(ss.v)) / norm(ss.v)
2 * dv_v0









    Out[20]:





$0.0029960538 \; \mathrm{}$



In [21]:

    
np.allclose(da_a0, 2 * dv_v0, rtol=1e-2)









    Out[21]:





True

This means we successfully validated the model against an extremely simple orbit transfer with approximate analytical solution. Notice that the final eccentricity, as originally noticed by Edelbaum, is nonzero:



In [22]:

    
ss_final.ecc









    Out[22]:





$6.6621427 \times 10^{-6} \; \mathrm{}$

References

[Edelbaum, 1961] "Propulsion requirements for controllable satellites"