Calculating GEO Satellite Longitude from its (t; RA, Dec) observations

[mikolajkarawacki]

Hi :) I've been - unsuccessfully so far - trying to calculate satellite longitude given its observed location and time. After failing to obtain correct results using this formula https://space.stackexchange.com/questions/22850/computing-geo-satellites-longitude-from-elevation-azimuth-from-a-given-latitude , I realized I could use Astropy + Skyfield. Unfortunatelly, I cannot use SkyCoord(ra,dec,obstime,distance).transform_to(GCRS).cartesian (distance is of course converted from AU to km ;) as a geocentric input vector for wgs84.latlon_of() - as my code crashes with errors on .center attribute missing.
My attempts at solving it in "plain" Skyfield with position_of_radec(ra,dec,distance) also crash on the shores of .center - even with manually defining pos.center=399. I also tried forcing it with positionlib.Geocentric(pos) but all I get for my troubles is ValueError: setting an array element with a sequence.

I would really appreciate any and all help or suggestions as to how to solve this problem. The deadline for our department's paper for the AMOS conference is getting uncomfortably close :(

[brandon-rhodes]

Could you share a snippet of code that you would like to get working, along with the exact exception?

[brandon-rhodes]

I'm working with a DataFrame containing Tracking Data Messages (TDMs), in which I have both the observatory's location and the satellite's (t; RA, Dec) as well as e.g. the satellite's NORAD ID, thru which I can get the apo- and perigeum => SemiMajor Axis ~ distance.

So — let me think through this a little more. If the RA and dec are from the observatory's point of view, then, given the position of the observatory and the time t, we can draw a line: it starts at the observatory's position, then heads off towards deep space in the direction of the given RA and dec.

And the satellite is somewhere along that line.

Each position along that line is a different distance from the center of the Earth. And the TLE should let you define a rough range within which you'd expect the satellite's altitude to fall — as you say, somewhere between its typical perigee and apogee. So by combining the two, you could narrow down the line to just the segment bounded by the perigee and apogee distances, and say "the satellite at time t was somewhere along this line segment."

You'd have two ways to solve for that line segment, I'm thinking?

1. Brute force: generate a few hundred points out along the line, then compute the altitude of each one, then filter to find the ones that fall within the perigee-apogee range.
2. Geometry: sit down with paper and pencil and come up with a formula for altitude along an arbitrary line segment, then turn the formula around so that it takes altitude as a parameter and returns a position on the line. Calling the reverse formula twice, you could find the endpoints for perigee and apogee.

Skyfield does not have a solver like this built in, that lets you say "object X was at this RA and dec from my observatory, and I know its altitude varies between Y and Z, so what line segment bounds its possible positions", so it makes sense that you weren't able to find a quick Skyfield call that would perform your computation for you. But we can probably put something together. Let me know if what I've been thinking here so far makes sense and matches your problem!

[mikolajkarawacki]

Since my last comment, I tried to do this in a purely AstroPy way, following a helpful tip from the Python users in Astronomy FB group that ended up as

lat, lon, h = float(row['COMMENT LATITUDE']), float(
        row['COMMENT LONGITUDE']), float(row['COMMENT ALTITUDE'])
    with solar_system_ephemeris.set('jpl'):
        SatCoords=SkyCoord(ra=row['Ra']*u.deg,
                           dec=row['Dec']*u.deg,obstime=row['Time'],
                           distance=row['SMA']*u.km,
                           location=EarthLocation.from_geodetic(
                               lat=lat*u.deg, lon=lon*u.deg, height=h*u.m))
        satGCRS=SatCoords.transform_to(GCRS)
        satITRS=satGCRS.transform_to(ITRS)
        satlon=satITRS.earth_location.lon.value

However, comparing the results of such a calculation to data available at https://www.n2yo.com/satellites/?c=10&srt=2&dir=0 for satellites identified (?) on our frames, gives huge discrepancies for all satellites that are in both my and n2yo datasets.

Right after noticing the discrepancies between n2yo and the aforementioned formula from space.stackexchange.com , my colleague and I tried to calculate satellite longitude from first principles but arrived at a formula that doesn't seem to be analytically solvable, vide attached files.


in those calculations, lambda_s is the satellite's longitude while phi and lambda are geographic coordinates of the observer.
I'm at a total loss as to what to do next with this problem. And I'm ever more nagged by the project's P.I. who seems to think that this is easy and straightforward...

[brandon-rhodes]

What does SMA stand for?

[mikolajkarawacki]

Oh, sorry, I thought it was obvious. SMA == Semi-Major Axis.

[brandon-rhodes]

Understood. In that case, I think that the distance= you provide to the SkyCoord needs to be the "range" distance from the observer to the satellite, rather than the SMA?

Which I think might be what you were working towards on your own whiteboard, with the line between P and S that's not quite the same length as the line from the Earth's center to S?