[GSoC2018|OrbitDeterminator|Arya] Week #7-8: Implementing Cowell propagator

Last time when we tried to use the SGP4 propagator in our Kalman Filter, we ran into problems. This was because you can’t stop the algorithm at a certain stage and continue from there. Suppose you want to propagate from $t_1$ to $t_2$ , you can directly propagate from $t_1$ to $t_2$ but you can’t propagate it from $t_1$ to $t_3$ to $t_2$ . So I set out to make our own propagator taking into account two major perturbations – oblateness of Earth and atmospheric drag. Instead of making an analytical propagator like SGP4, I went with a numerical one because that can be stopped and started at any point. So the first challenge was making a numerical integrator.

Numerical Integration

Numerical integration are used to integrate a differential equation from one point to another point. I studied and tested various integrators and finally settled on Runge-Kutta 4th order integrator (RK4).

Advantages of RK4:

Good accuracy for small to medium step sizes
Low computation cost
Easy to implement

Let the initial value problem be $\dot y = f(t,y)$ and $y(t_0) = y_0$ .

Choose a step size $h > 0$ and define

$\begin{align*} y_{n+1} &= y_n + \frac{1}{6}(k_1+2k_2+2k_3+k_4) \\ t_{n+1} &= t_n + h \end{align*}$

for $n = 0,1,2,3,\dots$ where,

$\begin{align*} k_1 &= h\,f\left( t_n,y_n \right) \\ k_2 &= h\,f\left( t_n + \frac{h}{2}, y_n + \frac{k_1}{2} \right) \\ k_3 &= h\,f\left( t_n + \frac{h}{2}, y_n + \frac{k_2}{2} \right) \\ k_4 &= h\,f\left( t_n + h, y_n + k_3 \right) \\ \end{align*}$

Note that these equations work equally well whether y is a scalar or a vector. In our case, the differential equation is $\mathbf{\ddot r}=-\frac{\mu}{r^3}\vec{\mathbf r} + \vec{\mathbf a}_p$ . It is not in a form suitable to be used directly in RK4. To use it in RK4 we have to vectorize the equation. First, let’s define the state of the system $\mathbf s$ as $[x, y, z, \dot x, \dot y, \dot z]^T$ . Let $r = \sqrt{x^2+y^2+z^2}$ . Then,

This form can be directly used in the RK4 integrator.

Modelling oblateness of the Earth

Spherical harmonics is used to model the gravity field around the Earth. The method is similar to how Taylor series works. Just like Taylor series, there are various coefficients for each degree of a term in the expansion. The most important coefficient is J₂. The next non-zero coefficient J₄ is 1000 times smaller than J₂. Considering only J₂, the gravitational potential around Earth is:

$U(r,\phi) = -\frac{\mu}{r}+\frac{1}{2} \frac{\mu}{r^3} J_2 R_{\oplus}^2(3\cos^2\phi -1)$

where,
$r$ is the radius from the centre
$\phi$ is the latitude of the point with respect to the equator
$J_2$ is a dimensionless constant equal to $0.001083$
$R_{\oplus}$ is the equatorial radius of the Earth, equal to $6378.1\, \text{km}$

To get the acceleration in xyz axes, find the gradient of U and transform coordinates to cartesian. Doing this, we get:

$\begin{align*} a_x &= \frac{3}{2}\frac{J_2 R_{\oplus}^2}{r^5}\mu x\left(5\frac{z^2}{r^2}-1\right ) \\ a_y &= \frac{3}{2}\frac{J_2 R_{\oplus}^2}{r^5}\mu y\left(5\frac{z^2}{r^2}-1\right ) \\ a_z &= \frac{3}{2}\frac{J_2 R_{\oplus}^2}{r^5}\mu z\left(5\frac{z^2}{r^2}-3\right ) \\ \end{align*}$

These equations implement J₂ perturbation. The main effect of J₂ perturbation is precession of the nodes and the argument of perigee.

Adding drag to the propagator

The simplest model for drag force is $\vec{\mathbf a} = -\frac{1}{2}\frac{C_DA}{m}\rho\,v_{rel}\vec{\mathbf v}_{rel}$ . We have to determine the air density, relative velocity and drag coefficients. For air density I used the model given on this website and adjusted the constants for the current solar cycle. The final air density model is:

$\rho = 0.6\exp\left(-\frac{1}{915}(h-175)(29.4-0.012h)\right) \text{kg}\, {\text km}^{-3}$

where h is the altitude of the satellite above the surface of the Earth.

To find the relative velocity, one must account for the rotation of the Earth. $\vec{\mathbf v}_{surface} = \vec{\mathbf \omega} \times \vec{\mathbf r}$ where $\vec{\mathbf \omega} = 7.29115 \times 10^{-5} \ \hat{\mathbf z} \ \text{rad/s}$ .

$\vec{\mathbf v}_{rel} = \vec{\mathbf v} - \vec{\mathbf v}_{surface}$

After putting together all this, I adjusted the drag coefficient so that satellites decay at approximately the current rate of decay. This finishes the implementation of drag perturbation.

Combining both the perturbations we get: $\ddot{\mathbf r} = -\frac{\mu}{r^3}\vec{\mathbf r} + \vec{\mathbf a}_{\text{oblateness}} + \vec{\mathbf a}_{\text{drag}}$
Putting it in the RK4 integrator completes the propagator.

Problems faced

I tried the algorithm on a recent TLE of the ISS (shown below).

ISS (ZARYA)             
1 25544U 98067A   18190.66798432  .00001104  00000-0  24004-4 0  9994
2 25544  51.6418 266.6543 0003456 290.0933 212.4518 15.54021918122018

According to this, it makes 15.54021 revolutions a day. This translates to a time period of 5559.8 secs. However, according to other simulators I found online, the actual period is around 5556 secs. A difference of 4 secs might not sound like much, but over the course of a day it adds up. The Wikipedia page of the ISS lists the time period to be 92.49 mins (5549.4 secs) but also states that it makes 15.54 orbits in a day (5559.8 secs). This is a contradiction. I searched for this issue and tried to correct it, but failed. Right now, the only method is to adjust the semi-major axis by a few kilometres to fix the time period. I will try to correct it later.

Future Work

Implementing the Kalman Filter
Testing against real world satellites

References

Autor: Arya Das

Computer science student, I like to code, fly planes, and learn about our beautiful universe. Alle Beiträge von Arya Das anzeigen