I have a measurement of the RA and Dec for an Earth orbiting satellite, as measured from a sensor on the Earth's surface. I'm trying to calculate the satellite's position vector in the GCRF reference frame (so a vector from the centre of the earth).
Since the object is in earth orbit, I can't assume that the RA/Dec is the same from the centre of the earth as it is from the surface. I therefore can't use this example: https://rhodesmill.org/skyfield/examples.html#what-latitude-and-longitude-is-beneath-this-right-ascension-and-declination
I have tried the following code.
from skyfield.api import load, wgs84, utc
from skyfield.positionlib import position_of_radec
from skyfield.units import Distance
from datetime import datetime
ra = 90
dec = 5
sensorlat = -30
sensorlon = 150
sensoralt = 1000
range = 37000
timestring = "2022-11-18T00:00:00.0Z"
distance = Distance(km=range)
time = datetime.strptime(timestring,'%Y-%m-%dT%H:%M:%S.%fZ')
time = time.replace(tzinfo=utc)
ts = load.timescale()
t = ts.from_datetime(time)
eph = load('de421.bsp')
earth = eph['earth']
sensor = wgs84.latlon(sensorlat,sensorlon,sensoralt)
satellite = position_of_radec(ra/15,dec,distance.au,t=t,center=sensor)
It appears that the sensor is represented as a Geocentric vector, while the satellite position is represented by a Geometric vector, and I can't figure out how to combine the two into a geocentric vector for the satellite.
satellite_icrf = sensor.at(t) + satellite
# Gives the following exception:
# Exception has occurred: TypeError
# unsupported operand type(s) for +: 'Geocentric' and 'Geometric'
I also tried simply changing the centre of the satellite Geometric vector, but that didn't appear to change the vector in any way.
print(satellite.position.km) # Prints [2.25697530e-12 3.68592038e+04 3.22476248e+03]
satellite.center = earth
print(satellite.position.km) # Still prints [2.25697530e-12 3.68592038e+04 3.22476248e+03]
Can someone explain how to convert from this Geometric vector into a GCRF vector?
I think I've found a workaround.
sensor.at(t) + satellite
does not work, but it looks like:
-(-sensor.at(t)-satellite) does work, and gives the required GCRF vector for the satellite.
This seems a bit hacky though, surely there's a more 'correct' method. I won't mark this as the accepted answer just yet, but I will if no one else posts a better method.