osculating_elements values for an EarthSatellite don't make any sense

[bsugerman]

I'm using skyfield to try to map out where the ISS is as a function of time starting with a given TLE. The most basic data that I'm looking for is the true anomaly as a function of time. I've created a simple python method to do this, and the results are insane. I've included the code and result. Since the ISS is in almost circular orbit, the anomaly should vary almost linearly between 0 and 360 over an orbit, and just repeat. Instead, I get a crazy pattern. I can't figure out I'm using skyfield and the osculating elements wrong, or if there is a bug in the skyfield code. I expect the osculating_elements is just for objects orbiting the sun, not the earth?

My idea below is to compute the osculating_elements for the satellite at a given time, then increment the time and recompute.

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import EarthSatellite, load, Topos
from skyfield.elementslib import osculating_elements_of
import datetime

class Satellite:
    def __init__(self,TLE,name='ISS (ZARYA)',time=None):
        self.satellite=EarthSatellite(TLE[0], TLE[1], name=name)
        self.ts = load.timescale()
        if time is not None:
            self.set_time(time)
        else:
            self.tu=self.satellite.epoch
        self.update()

    def update(self):
        # geocentric measurements class
        self.geocentric=self.satellite.at(self.tu) 
        # orbital elements
        self.elements=osculating_elements_of(self.geocentric) 
        # true anomaly of orbit = angle of satellite on its orbit
        self.alpha=self.elements.true_anomaly.radians 

    def set_time(self,time): # time is a datetime structure
        self.tu = self.ts.tt(time.year, time.month, time.day, 
                              time.hour, time.minute, time.second)
        self.update()

    def add_time(self,dt):
        self.tu=self.ts.tt(jd=self.tu.tt+dt/86400.)  
        self.update()

# TLE for the ISS
TLE=['1 25544U 98067A   15012.37614128  .00015029  00000-0  24160-3 0  9994',
 '2 25544  51.6468 141.8063 0006103 243.0591 165.4479 15.53207751923791']

s=Satellite(TLE) 
t=[]
a=[]
dt=10
for i in np.arange(0,8000,dt):
    t.append(i)
    s.add_time(dt)
    a.append(s.alpha)

plt.plot(t,a)
plt.show()

[JoshPaterson]

I must admit this initially seemed like a big error to me too. It makes sense that an object in a circular orbit should have a true anomaly that increases or decreases roughly linearly. I realized though that we can test this against data from JPL Horizons. Here's the plot I made with that:

I think the reason this is happening is because the orbit is nearly circular. True anomaly is the angle between the object and periapsis, but for circular orbits the periapsis is not well defined. For nearly circular orbits small perturbations can have a large effect on the argument of periapsis, which means small perturbations would also have a large effect on true anomaly.

A solution to this problem might be to use the argument of latitude instead. This is the angle between the object and the ascending node (which is well defined because the inclination of the ISS is not near 0). Argument of latitude for the example you gave does increase approximately linearly.

[ghost]

Neptune has the same issue: https://astronomy.stackexchange.com/questions/28691/why-is-neptune-true-anomaly-decreasing/28704#28704

[JoshPaterson]

One other thing is that all of Skyfield's functions are vectorized, which the code above could take advantage of. You can create a Time object that contains an array of times, and then get arrays of results back. This code below produces the same plot as the example above:

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import EarthSatellite, load
from skyfield.elementslib import osculating_elements_of

ts = load.timescale()

# TLE for the ISS
TLE=['1 25544U 98067A   15012.37614128  .00015029  00000-0  24160-3 0  9994',
 '2 25544  51.6468 141.8063 0006103 243.0591 165.4479 15.53207751923791']

ISS = EarthSatellite(TLE[0], TLE[1], name='ISS (ZARYA)')

seconds = np.arange(0, 8000, 10)
time = ts.tt_jd(ISS.epoch.tt + seconds/86400)
elements = osculating_elements_of(ISS.at(time))

plt.plot(seconds, elements.true_anomaly.radians)
plt.show()

[bsugerman]

What I ended up doing is computing the normal vector through the arg of perigee given in the TLE (the first row of the Euler angle rotation for the orbital elements). Then I use your code to propagate the orbit forward satellite.at(t), and compute the angle between its position vector and the perigee vector. Might be worth making that kind of option available for orbits that are close to circular?

-ben

Pardon any type-o's. The iPhone gnomes were grumpy.

On Jun 14, 2019, at 11:23 PM, Josh Paterson notifications@github.com wrote:

I must admit this initially seemed like a big error to me too. It makes sense that an object in a circular orbit should have a true anomaly that increases or decreases roughly linearly. I realized though that we can test this against data from JPL Horizons. Here's the plot I made with that:

I think the reason this is happening is because the orbit is nearly circular. True anomaly is the angle between the object and periapsis, but for circular orbits the periapsis is not well defined. For nearly circular orbits small perturbations can have a large effect on the argument of periapsis, which means small perturbations would also have a large effect on true anomaly.

A solution to this problem might be to use the argument of latitude instead. This is the angle between the object and the ascending node (which is well defined because the inclination of the ISS is not near 0). Argument of latitude for the example you gave does increase approximately linearly.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[brandon-rhodes]

@bsugerman That's an interesting approach! Is it code that you can share here for us to consider as an alternative approach?

[MichaelRThompson]

I actually just explained a similar example here. True anomaly is hard to define for a near-circular orbit when perturbations are being taken into account because your perigee can move around quite a bit.

The approach by @bsugerman is interesting, something I haven't seen before.