0

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?

Michael
  • 31
  • 4
  • Welcome to [Stack Overflow.](https://stackoverflow.com/ "Stack Overflow") This is not a code-writing or tutoring service. We help solve problems with code you are writing, not open-ended requests for advice. Please edit your question to show what you have tried so far, and what specific problem you need help with. See the [How To Ask a Good Question](https://stackoverflow.com/help/how-to-ask "How To Ask a Good Question") page for details on how to best help us help you. **DO NOT** post images of code, links to code, data, error messages, etc. - copy or type the text into the question. – itprorh66 Dec 01 '22 at 19:13
  • 1
    I'm sorry, but I'm not sure what I've done wrong. This isn't an open-ended request for advice, I'm trying to solve a specific problem (How do I convert from topocentric ra/dec to gcrf using Skyfield?) I included two examples of things that I've already tried (adding the `sensor` and the `satellite` vectors directly, and changing the value of `satellite.center`) I didn't post any "images of code, links to code, data, error messages, etc." I read through the link you posted and didn't see anything obviously wrong with my question. If something's wrong, could you please be more specific? – Michael Dec 01 '22 at 20:40
  • 1
    And it's even trickier than you imply, because GCRF is independent of rotation, so you also have to take into account the TIME of your observation. – Tim Roberts Dec 01 '22 at 21:11
  • Is this a geosynchronous object? If so, then the GCRF position changes with time. Much of skyfield was designed for distant objects. – Tim Roberts Dec 01 '22 at 21:19
  • Yes, this is a geosynchronous object, however I have included the time of the observation as `t`. When creating the satellite object, I included the observation time in the `position_of_radec()` function (the `t=t` argument) so I would expect it to give the GCRF vector at the instant of the observation. – Michael Dec 01 '22 at 21:23
  • Skyfield does appear to handle near earth objects such as satellites, but all of the examples on the website seem to assume that you're starting with a TLE, rather than an observed RA/Dec. – Michael Dec 01 '22 at 21:24

1 Answers1

2

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.

Michael
  • 31
  • 4