I would like to try out Gibbs' method of preliminary orbit determination, which requires that from observations of a space object at the three successive times $t_1$, $t_2$ and $t_3$ ($t_1 < t_2 < t_3$) we obtain the geocentric position vectors $r_1$, $r_2$ and $r_3$.

Now, my question is: how do I obtain these geocentric position vectors? Could I for example obtain these from a computer program or a website that provides such information? And if yes, what would you recommend?

  • 1
    $\begingroup$ @uhoh thanks for the answer mate! while that sounds indeed cool, I'll have to skip those steps and jump right into calculating the orbital elements. I was thinking, maybe I can find a real time satellite tracker, and when the sat passes over, I could log the lon/lat, and somehow (magically :)) convert them to a position vector? I'm really new to astrodynamics, so this might not make sense at all! :) $\endgroup$ – lawndownunder Nov 25 '20 at 23:45
  • 1
    $\begingroup$ @uhoh ok, I edited it a little bit! Now being more specific, thanks for pointing it out though! $\endgroup$ – lawndownunder Nov 26 '20 at 0:15
  • 1
    $\begingroup$ I can figure my way around Python, although I have zero experience with "celestial calculations" and related libraries. :) What's on your mind? $\endgroup$ – lawndownunder Nov 26 '20 at 1:45
  • 1
    $\begingroup$ Yes! Let me give it a try first. Thank you for the info! :) $\endgroup$ – lawndownunder Nov 26 '20 at 2:04
  • 1
    $\begingroup$ Ok! I think I managed to find some results, but I'll need a confirmation. Can I post a comment here, or as a new answer below and then remove it if it's wrong, or edit it? $\endgroup$ – lawndownunder Nov 26 '20 at 22:53

This is my tentative answer to my own question (based on suggestions in the comments).

Obtaining three different geocentric position vectors (km) over 10 minutes

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

ts = load.timescale()
minutes = np.linspace(0, 10, 3)
t = ts.utc(2020, 11, 26, 0, minutes, 0)

l1 = '1 25544U 98067A   20331.71797510  .00003353  00000-0  68809-4 0  9995'
l2 = '2 25544  51.6456 264.2547 0001926  82.8378 352.3176 15.49071921257216'

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

geocentric = satellite.at(t)

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]]

print(f"r1 = {r1}")
print(f"r2 = {r2}")
print(f"r3 = {r3}")

Which will print out:

r1 = [1760.6277509111617, -6050.90821032532, 2539.4465361221073]
r2 = [2937.5864721748503, -4682.676035746447, 3946.062541936976]
r3 = [3781.709583052999, -2783.861629642808, 4904.312360626345]
  • 1
    $\begingroup$ Yep looks good to me! Since your TLE has a fixed epoch this might give bad results if you run it again in a year, so might instead use minutes=np.linspace(0, 10, 3) and t=ts.utc(2020, 11, 26, 0, minutes, 0) to get three vectors over ten minutes. Then positions will be a 2D numpy array so r1, r2, r3 = positions.T and/or x, y, z = positionswill parse it out. $\endgroup$ – uhoh Nov 26 '20 at 23:38
  • 1
    $\begingroup$ then try stdev=5 and shape= position_in_km.shape and noise = np.random.normal(scale=stdev, size=shape) and noisy_positions = position_in_km+noise and see what some Gaussian noise in your vectors can do. You will find that if your three measurements are too close to each other (say 1 minute apart) and your noise is large, Gibbs' method (or any other) may start getting some wild results! The farther apart they are i.e. the larger the orbit's arc you can sample, the better to compensate for real-world limitations on observations. $\endgroup$ – uhoh Nov 26 '20 at 23:48
  • 1
    $\begingroup$ Ok, while I play around with the Gaussian noise, I edited the answer. Amazing! I also removed elevation, since I think it's out of the scope of the question. $\endgroup$ – lawndownunder Nov 27 '20 at 0:49
  • 1
    $\begingroup$ No, I'll definitively wait for a few days! I'll give it a week. No rushing anywhere. Thank you very much for your guidance! Really helpful, I appreciate it! :) I'll edit the title right away! $\endgroup$ – lawndownunder Nov 27 '20 at 0:52
  • 1
    $\begingroup$ this is a little old, probably for Python 2 and needs tuning for 3, but it might be fun to try: space.stackexchange.com/a/25959/12102 $\endgroup$ – uhoh Nov 27 '20 at 5:59

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.