2
$\begingroup$

I am trying to calculate the Classical Orbital Elements of the ISS using Python, and then compare them with the "actual" ones provided on the Internet. However, my obtained data differs by a lot from the "actual" one.

I am using the ISS TLE from this website, which at the time I am writing this (2020-12-04 04:19) it says:

1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038
2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464

Epoch (UTC):    04 December 2020 19:59:17

Now, using a modified version of the code from the answer from this question: How do I obtain geocentric position vectors at three successive times before I use Gibbs' method?,

from skyfield.api import EarthSatellite, Topos, load
import numpy as np

ts = load.timescale()
minutes = np.linspace(0, 60, 3)

t = ts.utc(2020, 12, 4, 20, minutes, 0)

l1 = '1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038'
l2 = '2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464'

satellite = EarthSatellite(l1, l2, name='ISS (ZARYA)')

geocentric = satellite.at(t)

r1, r2, r3 = geocentric.position.km

print(f"x/y/z of r1: {r1}")
print(f"x/y/z of r2: {r2}")
print(f"x/y/z of r2: {r3}")

I was able to obtain these positional vectors:

x/y/z of r1: [-4755.01172716  4947.37677595   377.21687117] # at '2020-12-04 20:00:00 UTC'
x/y/z of r2: [-4851.24226498  -367.40786871  5177.39228056] # at '2020-12-04 20:30:00 UTC'
x/y/z of r2: [  268.17962108  4636.21164042 -4395.73178942] # at '2020-12-04 21:00:00 UTC'

After getting these, I used Gibbs' method to obtain the following v2

v2 = [3.0406; -6.1187; 2.3591]

At this point, I use poliastro's rv2coe function to find the COEs:

from poliastro.constants import GM_earth
from poliastro.core.elements import rv2coe
from astropy import units as u
import numpy as np

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

# r2 from the obtained positional vectors
r = np.array([-4851., -367., 5177.])

# v2 calculated from Gibbs' method
v = np.array([3.040, -6.118, 2.359])

p, ecc, inc, raan, argp, nu = rv2coe(k, r, v)

print("p:", p, "[km]")
print("ecc:", ecc)
print("inc:", np.rad2deg(inc), "[deg]")
print("raan:", np.rad2deg(raan), "[deg]")
print("argp:", np.rad2deg(argp), "[deg]")
print("nu:", np.rad2deg(nu), "[deg]")

Which resulted in:

p: 6613.627129899017 [km]
ecc: 0.06923948841179417
inc: 53.14716567039852 [deg]
raan: 131.4224364765848 [deg]
argp: 241.26105012632252 [deg]
nu: 184.34321702974734 [deg]

Now, these result differ a lot from the results shown on the website. At the time I am writing this, the results from the website mentioned above are:

Epoch (UTC):    04 December 2020 19:59:17
Eccentricity:   0.0001925
inclination:    51.6418°
perigee height:     418 km
apogee height:  420 km
right ascension of ascending node:  224.1133°
argument of perigee:    115.0430°
revolutions per day:    15.49154811
mean anomaly at epoch:  245.0920°
orbit number at epoch:  1846

My question being: why do my results differ so much? Did I miss something?

$\endgroup$
9
  • 1
    $\begingroup$ Yes you missed something easy to miss but very critical, and if you used any number of observation times besides three your script would have thrown an exception and you'd have figured it out right away! Look again carefully at this comment and when you notice the difference you can write an answer! :-) $\endgroup$ – uhoh Dec 4 '20 at 10:46
  • 2
    $\begingroup$ Two more suggestions; 1) if you had checked the altitudes of your three geocentric positions $r=\sqrt{x^2 + y^2 + z^2}$ you would have noticed something was way off right away. 2) If the Poliastro method doesn't return an error or flag or poor goodness of fit when the plane containing the three positions does not pass anywhere near the center of the Earth, then maybe you can ask about that separately! Perhaps that's a weakness of Gibbs' method itself, or perhaps the implementation. (I added the poliastro tag) $\endgroup$ – uhoh Dec 4 '20 at 11:03
  • 2
    $\begingroup$ Other problems: (1) The epoch time is 04 December 2020 19:59:17. SGP4 can be very time sensitive. (2) Too much truncation; issue #1 is part of this problem. You truncated the position and velocity to four places in your call to rv2coe. (3) You did not show how you calculated velocity. Maybe you made a mistake there. (4) Gibb's method doesn't quite apply here. A sanity check shows that you should expect two places of accuracy. (5) You are calculating Keplerian osculating elements 30+ minutes after the epoch. (6) This is the biggest: Two line elements are not Keplerian osculating elements. $\endgroup$ – David Hammen Dec 4 '20 at 11:41
  • 2
    $\begingroup$ What you should have done is to calculate the position at $n$ minutes before the epoch, at the epoch, and at $n$ minutes after the epoch, where $n$ is small but not too small. Your choice of 30 minutes is too large; this will introduce the non-Keplerian nature of the orbit. Gibb's method would then have yielded an estimate of the velocity at the epoch. While SGP4 is only good to five or so places of accuracy, you should have retained full double precision accuracy until the end. Do not truncate intermediate values. $\endgroup$ – David Hammen Dec 4 '20 at 11:56
  • 5
    $\begingroup$ And finally, you should not expect a perfect match, not even close. SGP4 accounts for Earth's oblateness and for atmospheric effects. Keplerian orbits (which is what rv2coe assumes) do not account for these effects. Two line elements are not osculating Keplerian elements. $\endgroup$ – David Hammen Dec 4 '20 at 11:58
1
$\begingroup$

(Ok, this is my attempt of answering my own question based on the suggestions in the comments, full credits to @uhoh and @David-Hammen)

The main issue

As pointed out by @uhoh, the main issue here is parsing of the positional vectors. This part:

r1, r2, r3 = geocentric.position.km

Parses r1 = [r1.x, r2.x, r3.x], r2 = [r1.y, r2.y, r3.y] and r3 = [r1.z, r2.z, r3.z], but you want to have r1 = [r1.x, r1.y, r1.z], etc..

So, the correct parsing would be:

x, y, z = geocentric.position.km

r1 = [x[0], y[0], z[0]]
r2 = [x[1], y[1], z[1]]
r3 = [x[2], y[2], z[2]]

How could this be spotted earlier earlier:

  • By checking the altitudes of the three geocentric positions.

What else could have gone wrong:

  • The epoch time is 04 December 2020 19:59:17. SGP4 can be very time sensitive.
  • Too much truncation
  • Calculating Keplerian osculating elements 30+ minutes after the epoch (this is way too large spacing!)
  • Gibb's method doesn't quite apply here. A sanity check shows that you should expect two places of accuracy.
  • TLEs are not Keplerian osculating elements, Keplerian orbits (which is what rv2coe assumes) do not account for effects such as Earth's oblateness and atmospheric effects.

So let's re-write this based on these suggestions:

  1. Obtaining three positional vectors over 10 minutes, where the satellite's epoch is the time of the second positional vector.
from skyfield.api import EarthSatellite, Topos, load
import numpy as np

# ISS
l1 = '1 25544U 98067A   20339.83283773  .00016717  00000-0  10270-3 0  9038'
l2 = '2 25544  51.6418 224.1133 0001925 115.0430 245.0920 15.49154811 18464'

ts = load.timescale()

t1 = ts.utc(2020, 12, 4, 19, 54, 17)
t2 = ts.utc(2020, 12, 4, 19, 59, 17) # Epoch
t3 = ts.utc(2020, 12, 4, 20, 4, 17)

time = [t1, t2, t3]

satellite = EarthSatellite(l1, l2, name='ISS (ZARYA)')

for index, t in enumerate(time):
    index = index + 1

    geocentric = satellite.at(t)
    x, y, z = geocentric.position.km
    mag = np.linalg.norm(np.array([x, y, z]))

    print(f"r{index} = [{x}; {y}; {z}];")
    print(f"r{index} magnitude: {mag}\n")

We get that:

r1 = [-5599.602701921242; -3436.3402635476973; -1756.6275176484512]; 
r1 magnitude: 6800.715040494088

r2 = [-4903.0245824528865; -4709.731428621669; 10.226801416237763]; # at '2020-12-04 19:59:17 UTC' <- Epoch
r2 magnitude: 6798.626682893481

r3 = [-3651.287058898108; -5449.812055367599; 1775.887224806899]; 
r3 magnitude: 6796.03737927766
  1. Using a GNU Octave implementation of the Gibbs' method, we obtain the velocity of the second positional vector with:
mu = 398600;

R1 = [-5599.602701921242; -3436.3402635476973; -1756.6275176484512];
R2 = [-4903.0245824528865; -4709.731428621669; 10.226801416237763];
R3 = [-3651.287058898108; -5449.812055367599; 1775.887224806899];

r1 = norm(R1);
r2 = norm(R2);
r3 = norm(R3);

c12 = cross(R1, R2);
c23 = cross(R2, R3);
c31 = cross(R3, R1);

N = r1*c23 + r2*c31 + r3*c12;
D = c12 + c23 + c31;

S = R1*(r2 - r3) + R2*(r3 - r1) + R3*(r1 - r2);

v = sqrt(mu./norm(N)./norm(D)).*(cross(D,R2)./r2 + S);

Which will result in:

v = [3.3074; -3.4186; 5.9972]

  1. Plugging the results in Poliastro's rv2coe
from poliastro.constants import GM_earth
from poliastro.core.elements import rv2coe

from astropy import units as u

import numpy as np

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

v = np.array([3.3073, -3.4186, 5.9972])
r = np.array([-4903.025, -4709.731, 10.227])

p, ecc, inc, raan, argp, nu = rv2coe(k, r, v)

print("Semi-latus Rectum of Parameter:", p, "[km]")
print("Eccentriity:", ecc)
print("Inclination", np.rad2deg(inc), "[deg]")
print("Right Ascension of Ascending Node:", np.rad2deg(raan), "[deg]")
print("Argument of Perigee:", np.rad2deg(argp), "[deg]")
print("True Anomaly:", np.rad2deg(nu), "[deg]")

Will results in:

Semi-latus Rectum of Parameter: 6794.203064821717 [km]
Eccentriity: 0.001220230040456721
Inclination 51.58120030298049 [deg]
Right Ascension of Ascending Node: 223.77968914802696 [deg]
Argument of Perigee: 122.33400836496219 [deg]
True Anomaly: 237.77599772309975 [deg]

When compared to the results obtained from the website mentioned in the question, i.e.:

Eccentricity:   0.0001925
inclination:    51.6418°
right ascension of ascending node:  224.1133°
argument of perigee:    115.0430°

Might seem "close", but it is not perfect. However, plotting the obtained results:

from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from poliastro.plotting.core import OrbitPlotter2D as op2d
from astropy.time import Time

epoch = Time("2020-12-4 19:59:17", scale="utc")

r = np.array([-4903.025, -4709.731, 10.227])
v = np.array([3.3073, -3.4186, 5.9972])

r = r * u.km
v = v * u.km / u.s

orbit = Orbit.from_vectors(Earth, r, v, epoch)

fig = op2d().plot(orbit)
fig.show()

Will give us:

enter image description here

Which is pretty close to the orbit plot shown at website from the question.

$\endgroup$
2
  • 1
    $\begingroup$ Looks great, you are on your way! Now you can refine further by experimenting with times centered near the TLE's epoch and how wide to spread them on either side, adding noise, etc. I think what you will see is that if you add a few km of noise to measurements that are only say 10 seconds apart (and repeat several times with different random numbers) you'll start seeing a lot of variation in final parameters, but if the measurements are 10 minutes apart the same amount of noise will add a lot less variability. $\endgroup$ – uhoh Dec 5 '20 at 4:59
  • 1
    $\begingroup$ I'm taking notes. That's a battle to fight another day. Which might be really soon. :-D $\endgroup$ – lawndownunder Dec 5 '20 at 5:05

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.