How do I obtain geocentric position vectors at three successive times before I use Gibbs' method?

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?

• @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! :) Nov 25, 2020 at 23:45
• @uhoh ok, I edited it a little bit! Now being more specific, thanks for pointing it out though! Nov 26, 2020 at 0:15
• I can figure my way around Python, although I have zero experience with "celestial calculations" and related libraries. :) What's on your mind? Nov 26, 2020 at 1:45
• Yes! Let me give it a try first. Thank you for the info! :) Nov 26, 2020 at 2:04
• 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? Nov 26, 2020 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

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]

• 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.
– uhoh
Nov 26, 2020 at 23:38
• 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.
– uhoh
Nov 26, 2020 at 23:48
• 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. Nov 27, 2020 at 0:49
• 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! Nov 27, 2020 at 0:52
• 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
– uhoh
Nov 27, 2020 at 5:59