Polar Motion

[mkbrewer]

I am trying to track down what I think may be a bug in SOFA regarding polar motion. In their function Pvtob(), they do the standard rotation of the Earth by xp around the y axis and then yp around the x axis. But then, when calculating the observed azimuth and elevation, they do another rotation of the topocentric CIRS apparent position of the source. This second rotation looks very much like the alternate approximation for polar motion described in equation 3.27-3 on page 140 of the 1992 edition of the Explanatory Supplement to the Astronomical Almanac, just missing a factor of tangent latitude. I have never heard of anyone applying both of these before.

I wrote my own solar system ephemeris program using SOFA. When I compare my output without polar motion to JPL Horizons, I see an offset of about 5 mas in declination and 9 meters in distance over the course of a day for the Moon. The azimuth and elevation, however, are good to the 4 decimal places printed by Horizons. When I compare my output with polar motion to JPL Horizons, the declination offset is reduced to 0.5 mas and the distance offset is reduced to 0.8 meters. The azimuth and elevation, however, are now off by about +/-1 arcsecond. If I disable the second rotation this azimuth and elevation offset disappears.

As a second check, I compared my data to Skyfield. Without polar motion in either, I am seeing an offset of about 0.5 mas in Dec and 1.5 mas in RA, a +/-0.8 meter offset in distance, about +/-17 mas in azimuth and about -3 to +9 mas in elevation. When I add polar motion to my data, but not to Skyfield, the Dec offset increases to about 5 mas, the RA offset ranges from -3.7 mas to +1.8 mas and the distance offset ranges from -8 to +2 meters. The azimuth and elevation offsets are greater than +/-1 arcseconds again consistent with Horizons. Then I tried adding polar motion to Skyfield also. This resulted in a Dec offset ranging from -11 to +0.f mas, a RA offset ranging from -5.5 to -0.3 mas, a distance offset of -4 to +5 meters and the same +/-1 arcsecond offset in azimuth and elevation.

I find it odd that adding polar motion to Skyfield resulted in no improvement in the RA and distance offsets and actually made the Dec offset worse. So my question is: is there a known problem with polar motion offsets in Skyfield?

btw, I used the same DE403 ephemeris that Horizion uses in both my program and Skyfield.

[brandon-rhodes]

I am not aware that anyone has ever used the polar motion terms in Skyfield — so this is a good first test to see if they've been implemented correctly. Could you provide a small script that prints an altitude and azimuth, and that shows the difference you've observed between Skyfield's numbers and yours? That would give me a starting point for seeing how the polar motion terms affect your numbers and how they affect Skyfield's.

[ghost]

Not sure how helpful this is but https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/ explains how CSPICE deals with polar motion (and other high precision Earth orientation parameters).

[mkbrewer]

skyfield_tests.gz

This file is actually a compressed tarball. It contains the script I used to create the Skyfield data which consists of two files: one with no polar motion and the other with polar motion. I also included three files output by my program. One with no polar motion, one with standard SOFA polar motion and one where I zeroed out the two SOFA star independent parameters xpl and ypl to disable the additional polar motion correction that occurs when calculating observed AltAz in atioq(). This third file is the one that gives the best agreement with JPL Horizons. There is a mysterious 53 mas offset in apparent right ascension. I pointed this out to Jon Giorgini at JPL last June. His position was, of course, that it must be an issue with SOFA. I also pointed out to him that Skyfield agreed with SOFA on this and he just presumed that Skyfield must have been based on SOFA then. He did say that he would look into it, but nothing has changed with Horizons other than a note:

R.A.__(airless-appar)__DEC. =
  Airless apparent right ascension and declination of the target center with
respect to an instantaneous reference frame defined by the Earth equator
of-date (z-axis) and meridian containing the Earth equinox of-date (x-axis,
IAU76/80). Compensated for down-leg light-time delay, gravitational deflection
of light, stellar aberration, precession & nutation. Note: equinox (RA origin)
is offset -53 mas from the of-date frame defined by the IAU06/00a P & N system.

So I guess that either he hasn't yet found an issue with Horizons and/or he is still blaming SOFA and Skyfield.

There is also an Excel spreadsheet with plots of the difference between SOFA and Skyfield for 4 different tests:

1. No polar motion in either SOFA or Skyfield.
2. Polar motion in SOFA only.
3. Polar motion in both SOFA and Skyfield.
4. Polar motion in both where I modified SOFA as I explained above.

The main things to look at are differences in equatorial coordinates and distance as the offsets in azimuth and elevation are a side effect of these except in the case of the unmodified SOFA where I think the large offsets are a bug. If you don't have Excel, I could (with a little work) provide a PDF file with the plots only. I just thought a spreadsheet might be helpful for comparing my data with your own tests.

[mkbrewer]

@barrycarter Could you be a bit more specific? All I see there is a directory.

[brandon-rhodes]

@mkbrewer — Thank you for the script and the data! I don’t have Excel installed, but opened a quick IPython Notebook to look at the differences in the CSV files.

I should warn you that dives into small differences like this tend to take me a few weeks, as I always have other projects that need to be making headway as well. But each time I sit down to the problem I try to hang in there until I unearth a new item for investigation.

In this case, I found a first snag just as I was getting started on the simplest case: the two no-polar files. I was unhappy to see jumps in the difference between our equinox-of-date right ascensions:

I mean, I was unhappy enough with the fact that the graph is not centered on zero. And that it's different through the day. But the two places where there is a big jump made me extra unhappy.

It arises from a line of code I've never been quite happy with:

            limb_angle, nadir_angle = compute_limb_angle(
                target_au, gcrs_position)
            include_earth_deflection = nadir_angle >= 0.8

Here, following the USNO and their Astronomical Almanac, the decision about whether to include the effect of the Earth's gravitational deflection is made as a hard cut-off rather than a gradual one. Which produces the discontinuities we can see in the graph.

Since Skyfield only promises 0.5 mas agreement with the USNO's results, and the jump is roughly that size, I suppose I should just consider it noise. But I think this evening I'll go take a look at the Supplement and see their reasoning behind the hard cut-off. Maybe, since their business is to predict visible phenomena, they're not interested in spending time computing deflection for objects that are below the horizon and thus by definition won't be creating events in the Almanac at those moments?

In any case: as I have time, I'll continue posting comments here as I dig into each coordinate and circumstance and start understanding the behaviors I see. Hopefully it will only take me a few steps in before I reach the items you're actually interested in: the xp and yp values!

[mkbrewer]

I'll be watching. Particularly if you keep coming up with interesting tidbits like this. I guess it makes sense that there is no deflection when there is no light to deflect. Just so that you know, I'm using SOFA's Ld() which only takes the Sun into account. (I did account for the fact that the Sun to source vector is different from the observer to source vector which is neglected in SOFA's Ldsun()). Thus, it makes sense that there would be increasing divergence as the Moon approaches the horizon peaking at moonrise and moonset.

[ghost]

@barrycarter Could you be a bit more specific? All I see there is a directory.

The aareadme.txt file should explain more.

@brandon-rhodes If I remember correctly, you are parsing bpc kernels in other parts of Skyfield?

[brandon-rhodes]

I looked at the problem for only a few minutes today, so I'll merely report that the inclusion of Earth's gravitational deflection is one cause — rather than a remedy — of the difference in ICRS RA. Yesterday's graph didn't make it clear (at least to my eyes) whether the numbers diverged when the deflection was turning on or off. So here is the different between no-polar ICRS RAs if I turn Earth deflection "always on":

The effect is so large that it swamps the mean -1.0 mas difference, swinging the difference all the way from +1.0 mas to nearly -3.0 mas. I take it, therefore, that neither your library nor SOFA include Earth deflection.

And the effect, when turned on all the time, certainly does not looks physically realistic. Given the steep transient in the middle of the graph, is it treating the Earth as a point source — as a little atom at the Earth's center that weighs the same as the Earth? That would explain why light passing close to it was appreciably bent.

Contrariwise, if I turn Earth deflection all the way off:

~Our libraries become essentially and admirably identical, while both continuing to disagree with SOFA by -1.0 mas in RA in this example.~ Edit: our libraries agree more closely, and the wandering diurnal difference between them now lacks sharp edges.

For the record, turning off Earth deflection also ~makes our declination agree~ removes the sharp edges from the difference:

~While we both disagree with SOFA at around a 0.3 mas level in declination for this particular computation.~ In any case this leaves me with two action items for tomorrow:

1. Review whether Earth refraction is really an important effect — and if so, what can I do for objects behind the limb to hopefully stop having a sharp transition moment where the effect goes from "on" to "off".
2. Now that I know how to edit Skyfield to make our ICRS RA and dec agree perfectly, it's time to move on to apparent RA and dec. We'll see what I find.

[mkbrewer]

Actually, that disagreement due to gravitational light deflection only affects apparent RA. The blue in your plot is the disagreement between my implementation of ICRS RA and Skyfield's. I don't understand what you mean by "we both disagree with SOFA" as you are only comparing my computation to yours.

I would like to know what is causing the high frequency oscillation though. This is unique to Skyfield and is also quite apparent in the difference between our distances. It looks to me like Skyfield is quantizing some calculation without the interpolation necessary to eliminate the drift between calculations. This is, of course, a minor issue compared to the difference when polar motion is introduced.

[mkbrewer]

Btw, the remaining offsets are also seen in the comparison of my data with that of Horizons. I originally ascribed this to a difference in EOPs, but since this difference is also seen in the comparison of my data to Skyfield, I now believe that it's an issue with my data.

[brandon-rhodes]

I don't understand what you mean by "we both disagree with SOFA" as you are only comparing my computation to yours.

I was, alas, simply confused.

Do you happen to have a CSV for SOFA’s result? It would be most natural for me to compare all three libraries, and returning to the problem today I forgot — despite, alas, my careful axis labels — that all I have been comparing was Skyfield and your library.

So, to state things as I should have: it looks like all the transforms we do in going from astrometric to apparent are identical, as the difference between our libraries in astrometric RA is exactly the difference between them in apparent RA, except that Skyfield follows the USNO in including the Earth's gravitational deflection. But otherwise we must be taking the same steps — so, for example, our aberration routines must be equivalent.

The discussion of Earth deflection in the 2012 Explanatory Supplement to the Astronomical Almanac, surprisingly, indicates that they didn’t include it in their calculations: “Only the Sun’s gravitational field has been included … The gravitational field of the Earth, [which is] ignored here, can deflect light by a few tenths of a milliarcsecond for ground-based observers. It may also have to be allowed for in precise astronomy with Earth satellites.”

I can only suppose that “a few tenths of a milliarcsecond” wound up important for some of the applications for which NOVAS was designed.

I would like to know what is causing the high frequency oscillation though.

Good suggestion. I'll take a look at that next!

Btw, the remaining offsets are also seen in the comparison of my data with that of Horizons.

Interesting! Thanks for the additional graph. If you happen to have a CSV of the HORIZONS data, feel free to include it in this discussion.

[mkbrewer]

There is still some misunderstanding here. There is no solar system ephemeris functionality in the official SOFA distribution. I used a rather old version of jpl_eph.c that I've been using for about 20 years now (the one with functions like get_moon_posn_vel() if you are familiar with that version) and interfaced it to SOFA to create my program, so what you see is what you get in that respect. I like that old version since there is only one source file to deal with, no external dependencies and no makefile necessary either. I also use that version in the VxWorks OS to control several telescopes around the world. I do hope that there isn't some small inaccuracy in that version that is causing this offset. I can't imagine what it would be as it just extracts the data from the DE430 ephemeris.

I noticed that Horizons is now printing their equatorial coordinates to 9 decimal places instead of the 7 decimal places that they were using when I originally downloaded the data (I wish they would do that for az/el also as they currently print only 4 decimal places which makes it rather useless for this type of thing), so I downloaded it again and got a much better plot:

Here I have subtracted the 53 mas offset in apparent Right Ascension so that it doesn't mess up the scale (although it looks more like 52 mas to me). I see that they also include gravitational deflection by the Earth which cuts off when the Moon is below the horizon just as you do.

Here is the file that I downloaded:

moon_topo_4_6_2017_jpl_horizons.txt

I had to rename it to .txt due to github's limitations, but it's actually a .cvs file.

[mkbrewer]

Never mind on the 4 decimal places for az/el. I just noticed that it's now 6 so that's much better.

[brandon-rhodes]

There is still some misunderstanding here. There is no solar system ephemeris functionality in the official SOFA distribution.

Oh! So your library is a combination of both some JPL legacy code (or at least with jpl in the filename?) and the SOFA library. And so the systems under comparison are:

• Skyfield
• HORIZONS
• “mkb” = jpl_eph.c + SOFA + code of your own devising

Now that I’m clear on what we are comparing: my first big question is where is the ~1 mas difference is coming from between your library and Skyfield in the ICRF coordinates. That should be a much smaller difference, I think? The only ingredients are correctly computing a Chebyshev polynomial to generate x,y,z coordinates, and then converting them into spherical coordinates, both of which should be very high precision.

Skyfield agrees with NOVAS on those two tasks to within 0.2 arcseconds in RA and 0.1 arcseconds in declination.

I am trying to imagine where any difference could be cropping up between two libraries that compute the polynomials then convert to spherical coordinates. Could we be using somewhat different formulae for the difference TT – TDB, enough to make a full mas of difference?

I would be interested in comparing raw x,y,z coordinates. That would tell us whether the difference is creeping in before or after the switch to spherical coordinates.

In the meantime, I’m looking into the jaggedness that you noted in Skyfield’s coordinate curves, and might also take a look at its polar motion routine (your actual question!) and see why your polar figures didn’t have the effect you expected.

Never mind on the 4 decimal places for az/el. I just noticed that it's now 6 so that's much better.

Nice, I’m glad more digits are available! Is this with the “extra precision” checkbox checked in the HORIZONS “observer-table settings”?

[brandon-rhodes]

A quick update: the jaggedness might arise at least partly from Skyfield's use of a single 64-bit float for time rather than using two. I have been planning to upgrade that. Maybe it should be a project for later this weekend or next week, and then I could re-compare and see if the jaggedness persists.

[mkbrewer]

Oh! So your library is a combination of both some JPL legacy code (or at least with jpl in the filename?) and the SOFA library. And so the systems under comparison are:

Skyfield
HORIZONS
“mkb” = jpl_eph.c + SOFA + code of your own devising

Yes. jpl_eph just reads the ephemeris, so I need to calculate the down leg light time and iterate. I tested my iteration and adding one more iteration only made a difference of up to 2 mm in distance, so that part is fine. I'm using SOFA for the calculation of the observer's geocentric position. It uses the WGS84 ellipsoid, rotates by polar motion and then back rotates into the GCRS frame using the IAU 2006/2000A BPN matrix.

Skyfield agrees with NOVAS on those two tasks to within 0.2 arcseconds in RA and 0.1 arcseconds in declination.

I hope you meant to say mas here.

Could we be using somewhat different formulae for the difference TT – TDB, enough to make a full mas of difference?

This was the problem. I was using TT rather than TDB. I got to thinking about this last night and realized that since the Moon moves at a rate of ~0.5 arcsec / sec an error of 2 ms would indicate an error of 1 mas. I got up this morning and modified my program to use TDB by including SOFA's calculation of TDB-TT which contains 787 sinusoidal terms in time, 10 terms in orbital elements of the Sun, the Moon, Jupiter and Saturn and 5 terms of adjustment for planetary masses to be compatible with JPL rather than the IAU . This nearly eliminated the error in Ra/Dec and Az/El, but rather surprisingly had almost no effect on the error in distance.

There is still a bit of diurnal variation there that I can't account for.

I also compared my new data to Horizons. There it cut the error in both Ra/Dec and distance by about half. I looked into that and realized that Horizons actually uses the DE431 ephemeris which has a less accurate model of the Moon's orbit. I didn't want to download that as it's huge being divided in 1000 year segments rather than the 50 year DE430 ephemeris that I use. Here is a github compatible zip file containing the new data. This is labeled v5 as the former version of my program was v4.

skyfield_test.zip

Nice, I’m glad more digits are available! Is this with the “extra precision” checkbox checked in the HORIZONS “observer-table settings”?

Yes.

[mkbrewer]

A little further clarification on my "library". This is a single source file written in C++ that implements a class containing a few methods to stitch the JPL ephemeris into SOFA. The code of my own devising consists of a method for converting an array of unix timestamps supplied by the caller into UT1, TT and now TDB two part Julian dates, a method for calculating astrometric positions given a solor system body, a two part TDB Julian date and the observer's geocentric position vector supplied by SOFA and a method for calculating the gravitational deflection of light due to the Sun. This last method was necessary as SOFA's function for going from ICRS->CIRS unfortunately sets the target body's heliocentric position vector to be identical to the target body's topocentric vector. This is fine from the point of view of a distant star where the separation between the Earth and the Sun is negligible, but obviously horribly wrong for a solar system body. All I do there is calculate the Sun's light delayed barycentric position, subtract it from the target body's light delayed barycentric position and pass this vector along with the target body's topocentric position vector to a lower level SOFA function that does the calculation. There is some question as to what light delay to use for the Sun. My 1992 edition of the Explanatory Supplement uses the light delay from the target body to the observer, but SOFA suggests using the light delay from when the light makes its closest approach to the Sun. This would vary from no light delay for a target body in opposition to the Sun to the light delay associated the heliocentric position of the Earth for a target body in conjunction, but I figure that, since the Sun moves very slowly, it doesn't matter much at all, so I just went with the ES.

Speaking of light deflection, I had a thought on this:

Review whether Earth refraction is really an important effect — and if so, what can I do for objects behind the limb to hopefully stop having a sharp transition moment where the effect goes from "on" to "off".

Couldn't you use the cube of the sine of the zenith angle for zenith angles greater than 90 degrees to effectively shrink the volume of the Earth to zero at a zenith angle of 180 degrees?

Also, here's another plot of the az/el difference showing the true azimuth offset on the sky.

[brandon-rhodes]

I hope you meant to say mas here.

Indeed that’s what I meant!

This was the problem. I was using TT rather than TDB.

We found the problem! I’m glad the discrepancy could be tracked down to something solid and palpable.

This nearly eliminated the error in Ra/Dec and Az/El, but rather surprisingly had almost no effect on the error in distance.

Could that simply reflect the fact that the Moon’s motion is mostly transverse with respect to the Earth’s surface, moving it many meters around the Earth for each meter closer or farther?

There is still a bit of diurnal variation there that I can't account for.

Yes, a bit of sway, mostly in right ascension. It looks to be 90° out of phase from the difference in meters. Could it be that the ±0.4 meter difference in distance comes from a 0.4 meter displacement of where we think that geographic location is placed relative to the Earth’s center? A slightly different observation point would have no effect in RA when the Moon is in the direction of displacement, but (if I'm visualizing it correctly) would advance RA or retard it when the Moon’s direction was orthogonal to the difference vector.

A little further clarification on my "library". This is a single source file written in C++ … a two part TDB Julian date …

Thanks for the outline of how your single-file program is laid out! I have myself this weekend upgraded Skyfield to use two-part TT and TDB dates as well, and while it solved a problem with whether Skyfield preserves exact microseconds in round-trips to Python datetimes, it did not make my final output as smooth as yours! I might look at the problem a bit more today. It would be fun to have a smoother second derivative.

SOFA suggests using the light delay from when the light makes its closest approach to the Sun. … I figure that, since the Sun moves very slowly, it doesn't matter much at all, so I just went with the ES.

The effect is probably small, yes! Skyfield uses several other deflectors (like Jupiter), and does a single light-time iteration to use a position for the deflector when it was closest to the light ray.

Couldn't you use the cube of the sine of the zenith angle for zenith angles greater than 90 degrees to effectively shrink the volume of the Earth to zero at a zenith angle of 180 degrees?

That’s a very elegant idea, and I’m going to try it out! Hopefully it won’t knock all my tests against the NOVAS library out of whack — or maybe I can have a way to switch the effect off for running those tests.

[mkbrewer]

Yes, a bit of sway, mostly in right ascension. It looks to be 90° out of phase from the difference in meters. Could it be that the ±0.4 meter difference in distance comes from a 0.4 meter displacement of where we think that geographic location is placed relative to the Earth’s center? A slightly different observation point would have no effect in RA when the Moon is in the direction of displacement, but (if I'm visualizing it correctly) would advance RA or retard it when the Moon’s direction was orthogonal to the difference vector.

Likewise, a slightly different observation point would have no effect on distance when the Moon is orthogonal to the direction of displacement, but would advance distance or retard it when the.Moon is in the direction of displacement. Perhaps we have slightly different ellipsoids. Here are the parameters that SOFA uses:

case WGS84: a = 6378137.0; f = 1.0 / 298.257223563;

Regarding the noise in the position and distance, I've been also working with the Astropy people to get the many bugs out of their solar system code and now it's finally getting accurate enough that I see noise there also:

In their case, I'm pretty sure what's causing this is the way they go about calculating the light time in Python. Rather than doing what I do which is to loop on the given time array using the converged value of light_time(t[n-1]) as the initial value for light_time(t[n]), they just calculate the light time for the entire array and look to see that all of the values have converged to 1.0e-8 seconds (3 meters at the speed of light). My method, on the other hand, creates a smooth curve as the difference in light time from one time to the next is highly correlated and only dependent on the change in distance between the target body and the observer in the interval t[n] - t[n-1].

It would appear that Horizons does likewise as I get a smooth curve in the distance offset between my data and theirs:

Here, the diurnal variation is most certainly due to a slight difference in EOPs between the ones I used and those that Horizons uses.

[brandon-rhodes]

*f = 1.0 / 298.257223563;

The corresponding constant I see in Skyfield is:

IERS_2010_INVERSE_EARTH_FLATTENING = 298.25642

But swapping in your older constant, I don't see the diurnal curve disappearing. Alas.

I'm pretty sure what's causing this is the way they go about calculating the light time in Python.

Light time can certainly cause ugly discontinuities. I had thought that I'd fixed that in Skyfield by choosing to always use the same number of iterations, but looking at the code this morning I see:

        if -1e-12 < min(delta) and max(delta) < 1e-12:
            break

which looks adaptive rather than fixed. I should look into whether that causes sharp discontinuities that I could avoid. Your own approach sounds interesting, and I haven't seen code that does that before!

Your agreement with HORIZONS is looking very impressive in that last graph.

[brandon-rhodes]

Here, the diurnal variation is most certainly due to a slight difference in EOPs between the ones I used and those that Horizons uses.

@mkbrewer — So, where do you find EOPs? I think that one reason I've never seriously exercised Skyfield’s meager support for them is that I’ve never know where to find them. A search for “earth orientation parameters” takes me to:

https://www.iers.org/IERS/EN/Science/EarthRotation/EOP.html

But that site does not appear to have been updated much since the 1990s, except maybe for some editing of the "Coordinates of the pole" page in maybe the early 2000s.

The second link, though, looks more interesting and has several pages that I do not remember seeing before. A bit of clicking brings me to this:

http://hpiers.obspm.fr/eop-pc/index.php?index=PNUT&lang=en

Which has mostly graphs, though some other pages on that site seem to have tables of numbers behind them.

Do you know the most official or reliable source of orientation data in machine-readable form?

[ghost]

Do you know the most official or reliable source of orientation data in machine-readable form?

Brandon, I'm guessing my suggestion of https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/ (which is what CSPICE uses) wasn't helpful?

Also https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html (different URL on IERS site) seems to have daily updates?

[brandon-rhodes]

Brandon, I'm guessing my suggestion … wasn't helpful?

Oh — since you described it as documentation I did not go look at the directory myself. There are data files, it turns out! I’ll take a look at them when I get a chance.

[mkbrewer]

For the test script skyfield_topo_test.py included in: skyfield_test.zip The values there are from Astropy at 12:00 UT on April 6th, 2017. The output of this test script can be directly compared to moon_topo_4_6_2017_mkb_v5_polar_motion.csv as the same EOPs were used there.

[mkbrewer]

Indeed, this link:

Also https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html (different URL on IERS site) seems to have daily updates?

Has the daily updates. Bulletin B there is definitive, but it doesn't include future projections, of course. This is the data that is included in the yearly files at the same location.

[mkbrewer]

Please note that JPL Horizons does not use the IERS values. The name of the EOP file that is used for JPL Horizons output is listed in the output. You may contact JPL to have them send you the values that they used.

[brandon-rhodes]

@mkbrewer — I am interested in returning to this issue now that I’ve had a week of pursuing other loose ends. Could you remind me where things stand with respect to Skyfield’s support for xp and yp? Are we still at the point where they seem to bump the result in the wrong direction?

[mkbrewer]

Sorry for the delay. Yes, nothing has changed in that regard.

[brandon-rhodes]

Re-reading the NOVAS library's cel_pole() routine description, I’m not sure I had realized that it provides two different ways to correct the nutation model for unexpected polar motion. It allows either a pair of angles delta-delta-psi and delta-delta-epsilon, or else a dx and dy that correct the equivalent unit vector.

It looks like Skyfield has bits of both approaches scattered through the code. I might not have realized, as I copied ideas from NOVAS, that the two approaches accomplish the same correction. For example: there doesn't seem to be any protection in Skyfield from someone trying to apply both corrections, and thereby (presumably?) overcorrecting by accident in the other direction.

[mkbrewer]

Are you making any progress on this? I'd really like to be able to use Skyfield to resolve the differences I see between SOFA and JPL Horizons at some point.

[brandon-rhodes]

I have been making progress on several pieces of the altaz computation, but have not specifically made it all the way down to the small effects of polar motion.

Among other things (detailed in the changelog), I increased the precision of the ITRF → GCRS transform for velocity by several orders of magnitude in the past couple of weeks — previously there were too few digits of meaningful precision to even make polar motion visible in the result, so it seemed worth tackling first.

Also, over the past several weeks I have pushed through to nearly everything in the Time class the convention of using two floats for time instead of one. In particular UT1 is now computed as a pair of floats — which as you know is a crucial prerequisite to close alt/az agreement, as the Earth Rotation Angle is by far the fasted moving angle when an Earth location observes a Solar System object. Improving UT1 was a part of the overall improvement to ITRF → GCRS affecting not only velocity but also plain position itself.

With all of those quantities now tighter, I’ll now need to schedule time to re-compare Skyfield with HORIZONs and the results you shared from your own library. Hopefully time to tackle that will be in my schedule within the next few weeks!

[mkbrewer]

In response to a question I had from the astropy team, I made the mistake of looking at your documentation of apparent place. Now I'm well and truly confused. Your documentation states that your apparent place is in the GCRS. That can't be the case. GCRS defines two definitions of place. Proper place for geocentric position and local place for topocentric position. both wrt to the axes of the ICRF. Apparent place, both geocentric and topocentric, is always wrt the true equator and equinox of date. This is clearly the case for Skyfield as, if it wasn't, my apparent place without taking polar motion into account would not match yours.

[brandon-rhodes]

Your documentation states that your apparent place is in the GCRS.

Could you point me at the specific part of the docs you are reading? Thanks!

[mkbrewer]

https://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.Apparent

[brandon-rhodes]

I don’t see any mention there of apparent place. If I added a short section on the concept of “place” to https://rhodesmill.org/skyfield/positions.html do you think that other folks who are, like you, interested in “place” would be likely to find it? I would be happy to write up a clarifying few paragraphs about the relationship between what Skyfield calls a “position” and what older libraries call a “place”, if we can figure out where it would be most likely for users to find it.

[mkbrewer]

I'm not concerned with the semantics of place versus position. I use place as that is what the Explanatory Supplement uses. What I am concerned about is the statement: "In the specific case of an Earth observer, the output reference frame is the GCRS."

[brandon-rhodes]

In the specific case of an Earth observer, the output reference frame is the GCRS.

Ah, thank you for narrowing it down to that specific sentence! That helps a good deal. I think the issue is the semantics of place versus position. A Skyfield “position” is a vector, and the x,y,z vector returned by .apparent() does indeed have its three axes oriented to the three axes of the GCRS. In generating this vector, Skyfield is following the NOVAS library:

https://github.com/brandon-rhodes/python-novas/raw/master/Cdist/NOVAS_C3.1_Guide.pdf

For the Naval Observatory, as explained starting in the sentence “some caution with the semantics is in order”:

• “Apparent place” is .apparent().radec(epoch='date')
• “Topocentric place” is apparent place when the observer on the Earth’s surface
• “Virtual place” is .apparent().radec()
• “Local place” is virtual place when the observer on the Earth’s surface

The above is the rough draft of the docs I ought to add, with perhaps this (hopefully clarifying) introductory idea: that positions are vectors and (in Skyfield) are for simplicity always aligned to the BCRS / GCRS. Whereas “places” are positions reduced to spherical coordinates, raising the question of J2000/ICRS or equinox-and-equator-of-date.

I will certainly mention in the new section on “places” the fact, which I think may already be mentioned in the documentation, that it’s by far the most common to express Astrometric with .radec() and Apparent with .radec(epoch='date'). In addition I have tried to inculcate that habit in users by always using those vectors with those choices of RA/dec equinox; if you see any exceptions in Skyfield’s examples, please let me know and I’ll see if they can be fixed.

[mkbrewer]

OK. Maybe it is semantics then. There are two ways of looking at this. One is that the BPN matrix rotates the vector to a new position in the original coordinate system (ICRS) by adding the offset of the pole, equinox and frame bias from J2000 to the given epoch. The other is that the BPN leaves the vector constant and rotates the coordinate system itself so that the vector is now aligned wrt the pole and equator of date. This second view changes the axes so that they are no longer aligned with the GCRS. The main problem, in my view however, is that the documenation makes no mention of an epoch at all. This leads one (me) to mistakenly think that your apparent place/position is always wrt to the equator and equinox of J2000. Btw, the ES tut-tuts the concept of "virtual place" with a footnote under proper place that states: "Virtual place has exactly the same meaning as proper place. The term virtual place has been used in some U.S. Navel Observatory publications and products; it is mentioned here for completeness. To avoid confusion, this term will not be used in the future." Having no dog in that fight, that's all I'm going to say about that.

[mkbrewer]

Now I'm probably in trouble with the PETA people for using that expression. Hopefully, they don't read these comments. ;).

[mkbrewer]

Btw, I can confirm that the second view is the one that is used by both SOFA and the USNO. If you look at the definition of the rotation matrices on page xi of USNO circular 179, these are passive or alias transforms where the coordinate system (i.e. the Earth) is rotated wrt to the position vector. (see https://en.wikipedia.org/wiki/Active_and_passive_transformation).

[brandon-rhodes]

@mkbrewer — I have found my source for calling the raw GCRS vector an “apparent position”! (Note, not “place”, but “position”.)

… the term apparent place has traditionally been reserved specifically for the star or planet position we obtain by applying all these effects (except refraction) for a geocentric observer, with the coordinates expressed with respect to the true equator and equinox of date. … If the apparent star or planet position is instead expressed with respect to the mean equator and equinox of J2000.0 …” [emphasis mine]

— https://github.com/brandon-rhodes/python-novas/raw/master/Cdist/NOVAS_C3.1_Guide.pdf

Both in this paragraph and, if I recall, in their code, they distinguish “position”, which is always a vector, with “place”, which is always an RA/dec coordinate. And in the above paragraph, they clearly imagine there being but a single apparent position vector, which is the variously expressed in an of-date or a fixed coordinate system.

But as their idea of an “apparent position” vector might be idiosyncratic, I have updated the docs to try to improve those class descriptions you were looking at:

cd008c4e76760aec1428b7836867e355fa56a27d

In particular, I have both removed a reference to the idea of an “apparent position”, and have also mentioned the “places” that you asked about here.

Let me know if you find the revision any improvement.

[mkbrewer]

I think that is a lot more clear. I did leave a few comments on your commit.

[mkbrewer]

Actually, this is more of a bug than a feature request as Skyfield purportedly supports polar motion already.

[brandon-rhodes]

@mkbrewer — Happily, I found enough time last month to switch Skyfield to finals2000A.all to drive ∆T (instead of the old USNO deltat.data), which moves another step closer to using the same file for its polar motion numbers.

I have just made a thorough review of novas.c.

I think I had never properly internalized that there is not one but two systems of polar motion correction inside NOVAS. I am going to review the literature, which will probably make the difference clear, though my initial reaction is of course one of confusion: if the pole moves, why not have a pair of numbers to describe the motion and leave it at that? Why two pairs?

Here are my notes on the two systems as reflected in the code and documentation.

I am assuming that, in the short term, it's what I've arbitrarily called “System I” that you are most interested in seeing added to Skyfield next?

System I

From their NOVAS_C3.1_Guide.pdf:

“Not all aspects of the Earth’s orientation are predictable. Polar motion represents the small shift of the geodetic north pole (the ITRS z-axis) with respect to the rotational axis (the CIP), the largest part of which must be determined from observations. Typically, the total shift amounts to a few tenths of an arcsecond (1-2 µrad, 10 meters) and is specified by the parameters xp and yp. The observational determinations are designated simply as x and y, and current values are available from IERS Bulletin A.”

This correction occurs in their wobble() function which accepts xp and yp as arguments and corrects an ITRS vector for polar motion. It is called by cel2ter(), ter2cel(), and (transitively) equ2hor(), which themselves all accept xp and yp arguments. There appears to be no mechanism for storing xp and yp between calls; they always have to be resupplied by the caller.

System II

From their NOVAS_C3.1_Guide.pdf:

“Finally, the new IAU precession and nutation models are neither perfect nor complete and, for very high-accuracy applications, observational corrections are sometimes needed. These corrections now amount to less than one milliarcsecond (5 nrad). These corrections are available from the same sources as the polar motion determinations, and are designated as dX and dY (note that the units are milliarcseconds). In the rare cases where they are needed, they are pre-specified to NOVAS for use in subsequent calculations for a specific date. [cel_pole]”

“cel_pole can now accept input corrections to the pole positions as either (dψ, dε) or (dX, dY)”

“Important: For compatibility with the precession and nutation models used in NOVAS 3.0, and later, specify type = 2 and use only IERS dX and dY values with respect to “IAU 2000A” (sometimes labeled “IAU 2000”). These pole offset values will generally not exceed 0.5 milliarcsecond and therefore cel_pole would need to be called only when very high accuracy is required.”

Instead of being provided as an argument, this correction is performed by setting two global variables in the C module.

Setting the globals:

cel_pole() accepts either (delta-delta-psi,delta-delta-epsilon) or (dx,dy) arguments, and in response sets the globals PSI_COR and EPS_COR. In the former case, only a change of units needs to be performed, but if the GCRS (dx,dy) are given then a transformation is performed using the frame tie and precession.

Using the globals:

e_tilt() calls nutation_angles() to learn d_psi and d_eps, then adds PSI_COR and EPS_COR to them before returning them to the caller; they also affect the values it returns for the true obliquity and the equation of the equinoxes. The routines that use e_tilt() are the equatorial↔ecliptic transforms; sidereal_time(); geo_posvel(); and nutation() is basically “build a rotation matrix with the angles returned by e_tilt()”. Thus the values wind up affecting all sorts of calculations involving the Earth, as these routines are called everywhere else.

[mkbrewer]

Yes, system I is what we refer to when we speak of polar motion. This is the offset of the ITRS pole with respect to the CIP. What you are calling system II is the offset of the CIP itself from current theory of precession and nutation. This is why dX and dY are expressed in the GCRS. If you need a good reference, I'd recommend Chapter 5 of USNO Circular 179.

[brandon-rhodes]

My problem might be that I jumped directly to 5.4.4 in Circular 179 when I did my reading early this morning. I’ll return to the chapter and try reading the whole thing.

If only it had diagrams! Which makes me realize: I should also double-check that I’m not missing any good diagrams in the Supplement to the Almanac.

Thanks very much for the hint. From your description, I am imagining that one computes in this direction to figure out how a GCRS vector is oriented with respect to the Earth’s instantaneous position:

GCRS —(current theory’s deps and dspi + dX,dY)→ CIP —(xp,yp)→ ITRS pole, which is the physically realized position of the Earth’s crust per our observation stations.

Does that lay out the steps correctly?

[mkbrewer]

Yep, that looks good.

[mkbrewer]

@brandon-rhodes, here's an update on where we stand as of Skyfield version 1.31.

Skyfield_v1.31_SOFA_comparison.pdf

There are three sets of three plots each here. The first set is a comparison with no polar motion. The oscillatory behavior seen in the older version that I used (1.10) is gone. There is still a very small difference in the ICRS positions and moon distance which I don't understand. The only difference that I know of is that I'm using the version of JPL DE430 that is published in ASCII form. I suppose that there could possibly be some difference in TDB - TT given that we use different algorithms to compute that. The differences in apparent positions and Az/El are mainly due to the fact that I neglect light deflection by the Earth.

The second set of three plots show the effect of inclusion of polar motion in my program alone to give an idea of the magnitude of the effect. The third set is the comparison with polar motion included in both my program and Skyfield. As you can see, the residuals in right ascension and declination here are larger than those seen in the second set of plots indicating that inclusion of polar motion in Skyfield still makes things worse instead of better.

[brandon-rhodes]

… indicating that inclusion of polar motion in Skyfield still makes things worse instead of better.

As always, thanks for the thorough and detailed plots!

Looking over the version history, my sadly broken attempt to add xp and yp support to the Topos class dates from 2015, when I was for the first time trying to get agreement with HORIZONS and realized that the size of the remaining error was around the size of xp and yp. So I attempted to add the rotations, saw the error against HORIZONS drop, and declared victory.

It turns out, alas, to have been coincidence that my xp and yp happened to help rather than hurt that particular example.

Apologies for the italics and exclamation points in the account below, but I think they’ll aid my comprehension if I ever need to re-read this issue!

So:

I have just spent the evening looking for a way to test Skyfield against NOVAS, the library on which I based its astrometry — and I ran into a fundamental problem: NOVAS does not have the ability to apply xp and yp when computing the position of latitude and longitude in space!

So I had no example on which to base my code in toposlib.py.

This is a bit surprising, so I took a break from reading the code this afternoon, then returned later. But the conclusion held: the geo_posvel() routine in novas.c, on which I based my Topos, neither takes xp and yp as inputs nor pulls them from anywhere else. It’s innocent of polar motion corrections!

So I looked around and found another NOVAS routine, ter2cel(), which does indeed take xp and yp as arguments; then set up a test vector and had NOVAS rotate it; then worked on getting Skyfield to produce the same vector. After many failed attempts, I finally realized what is going on.

At the root of the problem was a conceptual mistake.

I had imagined that xp and yp described some residual error in where the Earth’s rotational axis points in space. With my usual image of the Earth in space as a spinning top, I imagined xp and yp as describing residual error in the orientation of its axis. (You can see that I was imagining xp and yp as what are actually deps and dpsi!)

So I had to realize that, no — the corrections xp and yp do not point at constant directions in space relative to the stars — rather, they point in constant directions across the Earth’s crust!

I had never even started to imagine such a thing. The corrections xp and yp are not motions of the axis of the spinning top at all. They are, rather, the motions of the surface of the top itself as that surface slides very slightly askew of the axis, sometimes slightly in the direction of Iceland, sometimes of Siberia, and so forth.

Which means that xp and yp are not corrections in a kinematically stable frame — they are corrections in the rotating frame of the Earth!

Which means that I was not simply applying xp and yp backwards. I was applying them in an entirely wrong reference frame — the equinox-and-equator-of-date frame, as it happens — and thus the corrections were spinning on a daily basis from the positions they really should be pointing. Whether they helped or hurt depended on the time of day of a particular calculation!

Because of the way NOVAS structures its latitude and longitude calculation, it turns out to be impossible to apply xp and yp correctly, so I am not surprised my attempt failed. Take a look at the code just above where I tried to apply xp and yp:

https://github.com/skyfielders/python-skyfield/blob/4cfc4c234e2f8bd65d7e82d5f21fd82ad6b0d229/skyfield/toposlib.py#L90

Do you see the call to the terra() routine? Just like the NOVAS terra() routine it’s based on, it skips ever computing an ITRF vector for the observer’s location by adding GAST to the longitude before doing the trigonometry to turn the angles into a vector. By the time a vector exists, it’s already in the equinox-and-equator-of-date frame, and any chance to do adjustments in the rotating ITRF frame have been entirely bypassed.

(I made a half-hearted attempt to apply xp and yp to the longitude and latitude before calling terra(), just for fun, but never quite got it to work.)

So I think that I’ll have to abandon terra() and instead build a vector when the Topos is created that gives the position in the ITRF of the observer. Then, in _at(), I’ll apply xp and yp, then use the GAST to spin the vector into the equinox-and-equator-of-date, and then apply t.MT as usual. It might be a more expensive calculation (which is maybe why NOVAS directly adds GAST to the longitude?) but unless the ITRF is a waypoint, xp and yp can’t be applied.

From all of the above, am I correct in drawing the following two conclusions?

1. The values xp and yp do not affect the matrix t.M (the product of the matrices NPB), because that matrix by definition moves coordinates into the frame whose pole is the CIO.

2. But xp and yp do need to be applied every time the target of a calculation is not the equinox-and-equator-of-date but is instead a real-life longitude and latitude where an observer (or target) is located.

Here is a script where for the first time I apply xp and yp successfully in a way that’s symmetrical with their treatment in NOVAS — now all I’ve got to do is get the Topos to behave the same way:

import numpy as np
from numpy import sin, cos
from novas import compat as novas
from skyfield import api
from skyfield.constants import ASEC2RAD, T0
from skyfield.functions import mxv, rot_x, rot_y, rot_z, tau

# xp_arcseconds = 0.0
# yp_arcseconds = 0.0

xp_arcseconds = 11.0
yp_arcseconds = 22.0

ts = api.load.timescale()
t = ts.utc(2020, 11, 11, 23, 37)

itrs_vector = [1.1, 1.2, 1.3]
v = itrs_vector

if yp_arcseconds or xp_arcseconds:
    sprime = -47.0e-6 * (t.tdb - T0) / 36525.0
    tiolon = -sprime * ASEC2RAD

    if 0:
        # Try it the NOVAS way.  More or less cut-and-pasted from novas.c:

        xpole = xp_arcseconds * ASEC2RAD
        ypole = yp_arcseconds * ASEC2RAD

        sinx = sin (xpole);
        cosx = cos (xpole);
        siny = sin (ypole);
        cosy = cos (ypole);
        sinl = sin (tiolon);
        cosl = cos (tiolon);

        xx =  cosx * cosl;
        yx =  sinx * siny * cosl + cosy * sinl;
        zx = -sinx * cosy * cosl + siny * sinl;
        xy = -cosx * sinl;
        yy = -sinx * siny * sinl + cosy * cosl;
        zy =  sinx * cosy * sinl + siny * cosl;
        xz =  sinx;
        yz = -cosx * siny;
        zz =  cosx * cosy;

        v = mxv(np.array([
            [xx, yx, zx],
            [xy, yy, zy],
            [xz, yz, zz],
        ]), v)

    else:
        # Try using NumPy matrix multiply.

        v = mxv(rot_x(-yp_arcseconds * ASEC2RAD), v)
        v = mxv(rot_y(-xp_arcseconds * ASEC2RAD), v)
        v = mxv(rot_z(-tiolon), v)

v = mxv(rot_z(t.gast / 24.0 * tau), v)
v = mxv(t.MT, v)
print(tuple(v))

EQUINOX = 1
FULL = 0
GCRS = 0
novas_km = novas.ter2cel(t.whole, t.ut1_fraction, t.delta_t,
                         xp_arcseconds, yp_arcseconds, itrs_vector,
                         method=EQUINOX, accuracy=FULL, option=GCRS)
print(novas_km)
print(v - novas_km)

[brandon-rhodes]

Update: I sat down this morning, rewrote how the Topos works to abandon the USNO library’s approach, and now that I had explicit ITRS vectors I was able to work in xp and yp correctly. If you are able to install Skyfield from source:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

— then I am interested in whether polar motion now helps rather than hurts Skyfield’s accuracy!