remkos
commented
7 years ago Dear Eric,
Your comments are all very valid. Let me give some responses.
(0) getorb works fine when an arc crosses the leap second. The arc itself will be time-stepped in physical steps of 60 or 30 seconds, and the software detects if that margin is crossed.
(1, 2) Indeed, there is no way to indicate on input a time during a leap second, so during that one second the interpolation is wrong. This has never been a big deal, since the altimeter production code generally fails as well, so we lost nothing.
(3, 4) The psuedo-cartesian was a way to store the orbits in a very compact manor. This format was created in 1992 at the time that there was no standard orbit format for interchange, so we made something that was small and easy to transfer across the internet. I have tested the interpolation accuracy extensively at the time, and it was always better than the digitisation of the data (1mm vertically and 1 cm horizontally). And that was even at the time that orbit errors were of the order of tens of centimeters.
Having said that, I long have wanted to replace this code with something more official and up to date. This would include:
- Reading the SP3 format, or something similar, preferably in netCDF format. Because there are still numerous formats around, it would be worthwhile to keep making a transformation into a common format. Indeed the best would be ECEF, and in double (float) precision, and including the velocity vector.
- Using a Hermitian interpolator (the one used in getorb has similar performance, but is not standard).
- The ability to use multiple time scales (UTC, TAI, GPS).
- The ability to use leap seconds on input, when using UTC.
I certainly invite volunteers to write such code.
rads_add_orbit calls getorb, which is not incorporated into RADS4.
I'm forwarding comments about getorb.f90 from Walter H.F. Smith: getorb.f90 has a block (lines 306 to 316) that checks whether the time step between two given ephemeris points in the file is not equal to the expected time step. If it is not, then a leap second is detected and from there on it looks like the logic will add one second.
I believe the getorb.f90 strategy therefore expects that positions in the file are time-tagged in UTC time and so if they are given every 60 [30] physical seconds with sampling initially on the UTC even [half] minute then they start out given as HH:MM:00 and after the leap will be time-tagged on the HH:MM:59 [HH:MM:30]. The strategy would then move the time tags after the leap so that time_b - time_a would compute correctly IF time_a (a local origin?) is before the leap.
A few things are not clear from that snippet of code:
(1) the time when the interpolation is desired would also have to be modified in the same way. How can this be done, since utc_time_a = 23:59:60.xxxx (a bit before the leap) and utc_time_b = 00:00:00.xxxx (next day, a bit after the leap) will both compute as the same time?
(2) what happens exactly on the leap? Presumably 23:59:60 is computed in UTC as equal to 24:00:00, i.e., midnight, and so it is the next 60 [30] second sample after that where one needs to begin adding 1 sec to each UTC time tag.
(3) pseudo-cartesian and (4) GRS80 ellipsoid These do not sound like they are precise-enough ways to compute lat,lon,height to 1 mm vertically and 1 cm horizontally from scratch. But if this routine is not used to compute the satellite position in its entirety, and instead is used to compute only a small perturbation, then maybe hard-wiring the GRS80 ellipsoid is OK, and maybe the pseudo-cartesian is OK. I would prefer to see ECEF in mm in quad precision, which is the way I do this in my CryoSat code.