Geodetic Ground-Station position in ICRS at time x

[concord647]

Hi Brandon, all,

I have a fairly simple problem to solve and I would like to verify my implementation in Skyfield, since my comparison with GMAT showed significant differences - although I am no expert on neither, so it could be very well the case that I made mistakes in my Skyfield-implementations/GMAT implementation or comparing the wrong things here...hence my question.

The task: I want to compute the cartesian coordinates in ICRS/J2000 (the differences between those two should be negligible for my purpose) of a ground station (that is defined via geodetic lat/lon/alt position) at specific times/time vectors.

This is my implementation using Skyfield:

# DEFINE GEODETIC GS COORDINATE SYSTEM
GS_Lat = 34.22475000  # N
GS_Lon = 118.05591667  # W
GS_Alt = 1728  # meter

GS = wgs84.latlon(GS_Lat  * N, GS_Lon  * W, elevation_m=GS_Alt )

ts = load.timescale()
timenow = ts.now()

# define start and end time in utc
time_start = ts.utc(year=2023, month=10, day=24, hour=5, minute=19, second=7.344)
time_end = ts.utc(year=2023, month=10, day=25, hour=5, minute=19, second=7.344)
steps = 24*60*60 # stepsize = 1 second

time_vector = ts.linspace(time_start, time_end, num=steps)

# get position of ground station at time_vector
position_vector_fullday = GS.at(time_vector)

# cartesian coordinates in ICRS in meter:
position_vector_fullday.position.m

Is my understanding correct that I produce cartesian coordinates in ICRS with this code?

I've tried to do the same thing in GMAT, i.e. defined the same ground station, propagated over a full 24 hours of the same epoch with 1 second intervall and reported the ground station cartesian coordinates in EarthICRF. The following plot shows the differences in all three axes of the two approaches:

So I am not sure where my mistake is - I guess my main question is: is this the way you would solve this problem using Skyfield or is that already wrong?

Thank you!

-Fabian

[brandon-rhodes]

I have not used GMAT before, but a first question I would think of asking is whether GMAT might be taking polar motion into account, which would shift the coordinates slightly. Maybe you could try adding code from the example shown here in the Skyfield documentation:

https://rhodesmill.org/skyfield/accuracy-efficiency.html#id1

[concord647]

I had already downloaded and used the latest finals.all unfortunately...

I ran both simulations for a three day period - starting at NOV-08, so this doesn't match the plot from my original post - and that confirms your guess:

[brandon-rhodes]

This almost looks like two sine waves going out of sync with each other — notice what close agreement they have at the start! I wonder if the waves would converge again periodically if carried out for many more days. So, it almost looks like the two time axes are simply using different scales. I would be tempted to double-check and triple-check whether both libraries are really giving you seconds. Maybe you could share how you are forming your x-axis parameter?

Given that you are generating 24*60*60 points for a range of 24*60*60 seconds, note that you are not producing one result per second—you are producing one point every 1.0000116 seconds. If you have a range of 3 seconds, you need 4 points in your linspace:

      second A           second B          second C
|                 |                 |                 |
0.0              1.0               2.0               3.0
^
four points

Maybe try a linspace of 86401 instead of 86400 points?

[concord647]

Oh man, that is embarrissing! You are right of course. I think I only checked the first and last timestamp and didn't realize what's happening in between...anyways, the plot now looks like this:

I guess the difference is now in the range of what is to be expected?

Thank you Brandon for your help, it's much appreciated!

[brandon-rhodes]

I guess the difference is now in the range of what is to be expected?

I'm glad we worked out why the two results were diverging!

The residual difference is a bit more than I would expect. If I'm doing the math right, the position is off by something like an arcsecond:

(30m / earthradius) / (2*pi) * 360 * 3600 = 0.97126581

I wonder if both libraries are using the same ∆T estimate? That would explain the x,y difference, as one position would be slightly trailing the other by the difference in their ∆T value. (If they were using different leap second tables, then an offset would also be created, but hopefully they are both up to date?)

Skyfield can return the exact TT for your starting time, and the exact delta T it is using. Does GMAT return the same values, so you could compare?

The constant difference in z is more concerning. How can they have different z if they are both using WGS84, I wonder?

[concord647]

I couldn't find anything about ∆T for GMAT. However, GMAT can report what they call "TTModJulian", which is a Modified Julian Date with an offset of 2430000.0. Skyfield doesn't seem to use the "Modified" Julian Date, so I adapted my time_vector:

time_vector.mjd = time_vector.tt - 2430000.0

and compared it to the GMAT report:

So it looks like there is a slight difference here, but if I am not mistaken, the maximum difference in MJD of 4.14e-8 should relate to 86400*4.14e-8=0.00357696 seconds?

It could explain the odd looking jump at the very end of the plot, but not the rest?

I also wanted to rule out any rounding errors of the WGS84-Position, as my original input used quite a few decimal places. I set a new position to

GS_Lat = 45 # N GS_Lon = 0 # W GS_Alt = 1000 # meter

and compared again:

The difference in z is actually not constant, but is slightly oscillating too, just on a much smaller scale:

[brandon-rhodes]

Given all the possible differences that could be in play, I am not sure what you should try next, unless you can read the source code of both libraries and add print() calls along the way to see their internal results until you uncover a difference?

Is there a constant time offset that, if applied to one of the two libraries, brings it into perfect agreement with the other and reduces all three lines to zero?

[concord647]

Yeah, I am not sure if I can/want to dive that deep into the differences of both libraries - I think the 20-30 meters discrepancy is fine for my purpose, even though it is a little irritating to not know where it stems from of course.

Anyways, I've had some one with access to STK generate the J2000 positions for the time period of the original post and compared them to both the GMAT and Skyfield results and none of them really match up:

Thanks for you support, it's much appreciated!

[brandon-rhodes]

Did you try adding or subtracting a little from the start time, to see if a constant time offset explains the difference — or at least the difference in x and y? It almost looks like what might be two versions of the same position, one following the other the whole way around the Earth.

Also, let's do a sanity check: using Pythagoras, what are the exact lengths of the three vectors? If they're not the same length from the Earth's center, then a first question to answer is why, which might solve some of your problem at least.

[concord647]

Yes, I did now. I found an "optimum" offset of -0.09201 seconds, i.e. Skyfield running 0.09201 ahead of GMAT. Nonetheless, that didn't get the error down to zero but it plateaued at ~ 7.3 meters. Not sure if I've only found a local optimum, but for larger offsets the error seems to grow significantly. I couldn't get the error < 7 meter.

Sanity check with Pythagoras: I've looked at the vector magnitudes and they seem to match fairly well. Interestingly, the STK vector is far less constant than with GMAT or Skyfield...

[brandon-rhodes]

Wow, that's close agreement on the length of the vector! I'm kind of amazed, actually, because the differences in the separate x,y,z coordinates are so far off, that I didn't imagine that the difference in length would have so many leading zeros. That's a very small number of meters.

I wonder if the two libraries are using the same precession-nutation routine? That's about the only thing I can think of that would knock the two vectors out of alignment by so far. Skyfield uses IAU 2000A, as explained here:

https://rhodesmill.org/skyfield/accuracy-efficiency.html#precession-and-nutation

If you make any progress on figuring out the difference, I'm very interested to hear — 20 meters is a huge difference for something like a known latitude and longitude! But unless someone makes a dive into GMAT's code, I'm probably out of suggestions for what to look at next to explain the difference. In any case, I hope your project proceeds smoothly!

[concord647]

Thanks a lot for your input and help! I am not sure how to narrow this down either, but will definitely let you know if I do find out more somehow.