4

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?

lawndownunder
  • 537
  • 4
  • 12
  • 1
    @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! :) – lawndownunder Nov 25 '20 at 23:45
  • 1
    @uhoh ok, I edited it a little bit! Now being more specific, thanks for pointing it out though! – lawndownunder Nov 26 '20 at 00:15
  • looks great now! – uhoh Nov 26 '20 at 00:28
  • Can you use some Python (programming language) or would you like to try? – uhoh Nov 26 '20 at 00:33
  • 1
    I can figure my way around Python, although I have zero experience with "celestial calculations" and related libraries. :) What's on your mind? – lawndownunder Nov 26 '20 at 01:45
  • If you'd like to try yourself, have a look, give it a try, and then go ahead and answer your own question! Else I'm happy to post a short script and you can modify it. https://rhodesmill.org/skyfield/positions.html#azimuth-and-altitude and https://rhodesmill.org/skyfield/earth-satellites.html This is the funnest and most satisfying Python package in the universe. It even has its own tag here and in Space SE – uhoh Nov 26 '20 at 01:55
  • 1
    Yes! Let me give it a try first. Thank you for the info! :) – lawndownunder Nov 26 '20 at 02:04
  • 1
    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? – lawndownunder Nov 26 '20 at 22:53
  • This is a pretty friendly site so I'd say just post an answer and begin with something like "Here's my tentative answer to my own question based on suggestions in comments, I'm posting it for review" or similar and it should be fine. – uhoh Nov 26 '20 at 23:02
  • 1
    Ok, nice! :)

    After the answer gets confirmed, I have a following question on that solution, which I'll post in this thread. This is fun. :)

     – lawndownunder Nov 26 '20 at 23:22

1 Answers1

4

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]
lawndownunder
  • 537
  • 4
  • 12
  • 1
    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 '20 at 23:38
  • 1
    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 '20 at 23:48
  • 1
    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. – lawndownunder Nov 27 '20 at 00:49
  • Great! If you feel this answers your question completely it's also okay to click accept. You could do that now, but if you wait a few days there's always a chance others will add additional helpful answers. Can you double check the title again, I think it reflects my partial mis-understanding of the scope of your question and needs to be adjusted as well. – uhoh Nov 27 '20 at 00:51
  • 1
    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!

     – lawndownunder Nov 27 '20 at 00:52
  • 1
    this is a little old, probably for Python 2 and needs tuning for 3, but it might be fun to try: https://space.stackexchange.com/a/25959/12102 – uhoh Nov 27 '20 at 05:59