ITRS Hour angle and declination

[mkbrewer]

To complement altaz, it would be good if Skyfield could also provide hour angle and declination (hadec) for an equatorially mounted telescope. This is the ITRS hour angle and declination taking into account both polar motion and atmospheric refraction.

[brandon-rhodes]

Am I correct that hour angle is simply right ascension relative to LST?

I would almost be tempted to say “a single subtraction is a simple enough operation that it might not be worth a separate method” — except that my own attempt to perform the simple subtraction is not getting close agreement with HORIZONS. Any idea where my script needs to be tweaked?

from skyfield.api import load, wgs84
from skyfield.data import iers

url = load.build_url('finals2000A.all')
with load.open(url) as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)

ts = load.timescale()
iers.install_polar_motion_table(ts, finals_data)

t = ts.utc(2020, 12, 19)
boston = wgs84.latlon(42.3567, 288.9431)

eph = load('de421.bsp')
earth = eph['Earth']
jupiter = eph['Jupiter Barycenter']

a = (earth + boston).at(t).observe(jupiter).apparent()
ra, dec, distance = a.radec('date')
ha = (boston.lst_hours_at(t) - ra.hours) % 24.0

print('Hour angle:', ha)
print('HORIZONS:  ', 4.988613031)

print('Difference (mas):', (ha - 4.988613031) * 15 * 3600 * 1e3)

The output:

Hour angle: 4.988615004962174
HORIZONS:   4.988613031
Difference (mas): 106.59395740475475

Here’s the HORIZONS output:

*******************************************************************************
Ephemeris / WWW_USER Sat Dec 19 14:59:30 2020 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Jupiter Barycenter (5)          {source: DE431mx}
Center body name: Earth (399)                     {source: DE431mx}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 2020-Dec-19 00:00:00.0000 UT      
Stop  time      : A.D. 2020-Dec-29 00:00:00.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_JUPITER_BARYCENTER          {West-longitude positive}
Target radii    : 71492.0 x 71492.0 x 66854.0 km  {Equator, meridian, pole}    
Center geodetic : 288.943100,42.3567000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 288.943100,4720.39581,4274.9653 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : High-precision EOP model        {East-longitude positive}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Sun
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE431mx}
Rel. light bend : Sun, EARTH                      {source: DE431mx}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : CAL 
EOP file        : eop.201218.p210311                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2020-DEC-18. PREDICTS-> 2021-MAR-10
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
*******************************************************************************
 Date__(UT)__HR:MN:SS     R.A._(a-appar)_DEC. L_Ap_Hour_Ang
***********************************************************
$$SOE
 2020-Dec-19 00:00:00  m  302.16423 -20.64094   4.988613031
 2020-Dec-20 00:00:00  m  302.39216 -20.59545   5.039127714
 2020-Dec-21 00:00:00  m  302.62093 -20.54948   5.089585822
 2020-Dec-22 00:00:00  m  302.85052 -20.50302   5.139989347
 2020-Dec-23 00:00:00  m  303.08089 -20.45608   5.190340223
 2020-Dec-24 00:00:00  m  303.31202 -20.40866   5.240640335
 2020-Dec-25 00:00:00  m  303.54390 -20.36077   5.290891515
 2020-Dec-26 00:00:00  m  303.77649 -20.31241   5.341095554
 2020-Dec-27 00:00:00  m  304.00977 -20.26358   5.391254202
 2020-Dec-28 00:00:00  m  304.24372 -20.21429   5.441369167
 2020-Dec-29 00:00:00  m  304.47831 -20.16454   5.491442121
$$EOE
*******************************************************************************
Column meaning:
 
 ** WARNING **

    An observer table has been selected for a barycenter. A barycenter is a
dynamical point, not a physical object. Any output physical quantities (such
as magnitude, radii, latitude & longitude) are based on the corresponding
planetary body.

    Barycenters are often close to, but not necessarily coincident with planet
centers. Thus, any output physical quantities may be slightly biased. It is
up to the user to determine if this bias is acceptable for the intended use.
 
TIME

  Times PRIOR to 1962 are UT1, a mean-solar time closely related to the
prior but now-deprecated GMT. Times AFTER 1962 are in UTC, the current
civil or "wall-clock" time-scale. UTC is kept within 0.9 seconds of UT1
using integer leap-seconds for 1972 and later years.

  Conversion from the internal Barycentric Dynamical Time (TDB) of solar
system dynamics to the non-uniform civil UT time-scale requested for output
has not been determined for UTC times after the next July or January 1st.
Therefore, the last known leap-second is used as a constant over future
intervals.

  Time tags refer to the UT time-scale conversion from TDB on Earth
regardless of observer location within the solar system, although clock
rates may differ due to the local gravity field and no analog to "UT"
may be defined for that location.

  Any 'b' symbol in the 1st-column denotes a B.C. date. First-column blank
(" ") denotes an A.D. date. Calendar dates prior to 1582-Oct-15 are in the
Julian calendar system. Later calendar dates are in the Gregorian system.

  NOTE: "n.a." in output means quantity "not available" at the print-time.
 
SOLAR PRESENCE (OBSERVING SITE)
  Time tag is followed by a blank, then a solar-presence symbol:

        '*'  Daylight (refracted solar upper-limb on or above apparent horizon)
        'C'  Civil twilight/dawn
        'N'  Nautical twilight/dawn
        'A'  Astronomical twilight/dawn
        ' '  Night OR geocentric ephemeris

LUNAR PRESENCE (OBSERVING SITE)
  The solar-presence symbol is immediately followed by a lunar-presence symbol:

        'm'  Refracted upper-limb of Moon on or above apparent horizon
        ' '  Refracted upper-limb of Moon below apparent horizon OR geocentric
             ephemeris
 
 'R.A._(a-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-dat
(z-axis) and meridian containing the Earth equinox of-date (x-axis, EOP-correct
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.

  Units: RA  in decimal degrees, ddd.fffff{ffff}
         DEC in decimal degrees  sdd.fffff{ffff}

 
 'L_Ap_Hour_Ang' =
   Local apparent HOUR ANGLE of target at observing site. The angle between the
observers' meridian plane, containing Earth's axis of-date and local zenith
direction, and a great circle passing through Earth's axis-of-date and the
targets' direction, measured westward from the zenith meridian to target
meridian along the equator. Negative values are angular times UNTIL transit.
Positive values are angular times SINCE transit. Exactly 24_hrs/360_degrees.
EARTH TOPOCENTRIC ONLY.  Units: sHH.fffffffff  (decimal angular hours)


 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information  : https://ssd.jpl.nasa.gov/
     Documentation: https://ssd.jpl.nasa.gov/?horizons_doc
     Connect      : https://ssd.jpl.nasa.gov/?horizons (browser)
                    telnet ssd.jpl.nasa.gov 6775       (command-line)
                    e-mail command interface available
                    Script and CGI interfaces available
     Author       : Jon.D.Giorgini@jpl.nasa.gov

*******************************************************************************

[Bernmeister]

Am I correct that hour angle is simply right ascension relative to LST?

I found these references in my travels a while back, might help:

• https://tycho.usno.navy.mil/sidereal.html
• http://www.wwu.edu/skywise/skymobile/skywatch.html
• https://www.heavens-above.com/whattime.aspx?lat=-33.8675&lng=151.207&loc=Sydney&alt=19&tz=AEST&cul=en

The magic formula being hourAngle = localSiderealTime - bodyRA (measured in radians).

[mkbrewer]

Hmm, that is unexpected. Please look at the difference in both apparent right ascension and LST. I'd expect them to be the same at about 51 mas. This should yield agreement in hour angle. I'll look into this on my side also. (see below)

However, I'm not asking for apparent hour angle and declination here. I'm asking for ITRS hour angle and declination as derived from azimuth and elevation using geodetic latitude. This is the hour angle and declination for an equatorial mount aligned with the ITRS pole rather than the CIP. In a non-airless model, it also includes refraction.

[mkbrewer]

@Bernmeister That is correct, however there is a slight rotation of the pole involved here due to polar motion. This is the difference between TIRS and ITRS hour angle.

[mkbrewer]

Yes. I'm seeing something similar on my side also. For my moon data, the difference in apparent right ascension between my data and that of Horizons is 51.27 mas, while the difference in LST is 110.66 mas. This yields a difference in apparent hour angle of 59.39 mas. Since LST is derived from UT1, I expect that this due to a slight difference in DUT1 between my data and that of Horizons.

[brandon-rhodes]

I'm not asking for apparent hour angle and declination here. I'm asking for ITRS hour angle and declination as derived from azimuth and elevation using geodetic latitude. This is the hour angle and declination for an equatorial mount aligned with the ITRS pole rather than the CIP.

In that case, am I correct that the angle you are asking for is not the angle that HORIZONS is determining? Your hour angle would be a rotation around the ITRS pole. Their hour angle, by contrast, seems to be a rotation around the “Earth's axis of-date” by which I think they would mean the rotational axis — from their HA definition above:

• “The angle between the observers' meridian plane, containing Earth's axis of-date and local zenith direction, and a great circle passing through Earth's axis-of-date and the targets' direction”

Their HA is the angle between two planes, both of which have in common the “Earth's axis of-date”, which I think means their HA can be computed as the difference between two right ascensions in the equinox-and-equator-of-date?

Your emphasizing “geodetic latitude” certainly points out a flaw in my first attempt above: I asked “where is Boston relative to the geocenter?” instead of “where does Boston’s zenith point?” which would not make a difference if polar motion was zero — since longitude, the operative quantity here, ignores Earth flattening and the difference between geodetic and geocentric coordinates — but does indeed matter in the case of polar motion, which knocks the ITRS slightly askew of the equator-and-equinox-of-date and mixes a little bit of longitude into latitude and vice versa in the rotation to the frame of reference of Earth’s rotation.

Alas, though, when I try factoring that into my script, moving my focus from Boston to its zenith, I get farther from the HORIZONS number rather than closer!

from skyfield.api import load, wgs84
from skyfield.data import iers
from skyfield.framelib import itrs
from skyfield.positionlib import Geocentric, Distance, Velocity, mxv

url = load.build_url('finals2000A.all')
with load.open(url) as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)

ts = load.timescale()
iers.install_polar_motion_table(ts, finals_data)

t = ts.utc(2020, 12, 19)
boston = wgs84.latlon(42.3567, 288.9431)

eph = load('de421.bsp')
earth = eph['Earth']
jupiter = eph['Jupiter Barycenter']

# "The angle between the observers' meridian plane, containing Earth's
#  axis of-date and local zenith direction,"

def zenith_of(boston, t):
    zenith = [0,0,1]
    xyz = mxv(boston._R_latlon.T, zenith)
    d = Distance(au=xyz)
    v = Velocity([0,0,0])
    return Geocentric.from_time_and_frame_vectors(t, itrs, d, v)

ra3, dec3, distance3 = zenith_of(boston, t).radec('date')
print(ra3.hours, dec3.degrees)

# "and a great circle passing through Earth's axis-of-date and the
#  targets' direction"

a = (earth + boston).at(t).observe(jupiter).apparent()
ra2, dec2, distance2 = a.radec('date')

# Compute difference:

ha = (ra3.hours - ra2.hours) % 24.0
ha_horizons = 4.988613031

print('Hour angle:', ha)
print('HORIZONS:  ', ha_horizons)
print('Difference (mas):', (ha - ha_horizons) * 15 * 3600 * 1e3)

Output:

1.132898384637877 42.356784878153405
Hour angle: 4.9886153343651145
HORIZONS:   4.988613031
Difference (mas): 124.38171619066907

Instead of matching HORIZONS, let's instead focus on whether my computations can match yours. Can you share the exact hour angle you're getting for these circumstances, and maybe suggest where my calculation could be adjusted to match yours?

[mkbrewer]

OK. Here's my data from my usual run of the Moon on April 6th, 2017. Rather than azimuth and elevation, I am printing ITRS hour angle and declination (HAO, DECO) here for an airless model.

moon_topo_4_6_2017_mkb_sf_v5_hadec.csv.txt

If we take a Cartesian vector in az/el coordinates such that the X axis points South, the Y axis points West and the Z axis points up, then the difference between an Alt/Az mount and an equatorial mount is a counter-clockwise rotation by 90 - latitude around the Y axis such that the Z axis points to the ITRS pole.

Please note that in a non-airless model, this rotation must be done after atmospheric refraction has been applied.

[mkbrewer]

There will be a slight difference between your data and mine as I am neglecting the deflection of light in the Earth's gravitational field.

[brandon-rhodes]

Thanks for the data file! Unless someone else wants to jump in and try replicating the numbers using Skyfield code, I should have some time the week after Christmas to take a look. (And: happy holidays, everyone!)

[mkbrewer]

OK. This works for me using Skyfield version 1.34.

skyfield_topo_hadec_test.py.txt

While doing this, I noticed that I had forgotten that this is a left handed coordinate system. I amended my previous post accordingly.

[mkbrewer]

I also notice that your DUT1 doesn't quite agree with finals2000A.all. For 0 UT on the 6th, it's off by 4 ns increasing smoothly to 12 ns for 0 UT on the 7th. This is likely to be round off error as there is no interpolation at these times, but it seems a bit large for that.

[reza-ghazi]

Merry Christmas, and happy holidays everyone!

[mkbrewer]

@brandon-rhodes It's been nearly a month now. I don't want to push. I realize that you've been busy, just wondering if this somehow fell off your radar.

[brandon-rhodes]

Thanks for checking in! It's still on my radar. The holidays have been busy, and also one or two bugs appeared that needed attention before a new feature (or at least that's how I tend to prioritize).

[brandon-rhodes]

And sometimes I add features that help me answer questions that I myself have. I just landed new code that when the user runs python -m skyfield prints:

Skyfield version: 1.35
jplephem version: 2.14
sgp4 version: 2.13
Built-in leap seconds table ends with leap second at: 2016-12-31 23:59:60 UTC
Built-in ∆T table from finals2000A.all covers: 1973-01-01 to 2021-10-16

That will prevent me from having to double-check those numbers manually.

[brandon-rhodes]

If only GitHub had a way to attach a to-do list to an Issue. (Besides editing the original issue description to include it, which I avoid if I didn’t write the issue myself.) But here are TODO items from this issue at this point in time:

[Edit: switched from a numbered list to bullets, to not falsely imply that these have a necessary order; in particular, the first item will probably be the last tackled.]

• Further investigate the “it's off by 4 ns increasing smoothly to 12 ns” phenomenon for DUT1 interpolation. I think I opted for speed rather than high precision, since the DUT1 values from finals2000A.all are only specified to 7 digits (±100 ns) and have physical significance at least an order of magnitude lower (±0.0000066 seconds = 6600 ns for the most recent observation from the file). But it does look a little confusing to see extra digits hanging off the end that were not there in the finals2000A.all file, and if they could be eliminated without slowing things up, that would certainly be prettier. If I recall I made a bit of progress on that just before Christmas, but am not yet at a conclusion.
• It should have occurred to me months ago, but: let’s add a mechanism for turning off Earth deflection. But where should it live? Should it be an option to .apparent()? Or be a flag that lives on the timescale, which would be more convenient but for the first time would attach “policy” not directly related to time scales to the timescale object? And whether it’s a new option to .apparent() or a flag on the timescale, should it be part of a more general mechanism for specifying which deflectors (Jupiter, Saturn, etc) are used, or be a separate little flag of its own?
• Once I can control Earth deflection I should double-check whether results all across Skyfield suddenly get closer to HORIZONS results, since I think it might not include Earth deflection? @mkbrewer, is that your impression of the HORIZONS behavior as well?
• In any case, I want then to figure out where HORIZONS is getting its hour angle from so that I can match its results.
• It will be interesting to compare that result to the code you provided (thanks, @mkbrewer!) for an hour angle generated not by subtracting RA and LST, but by rotating altaz (possibly affected by atmospheric refraction) back into a spherical coordinate system anchored at the geographic (not rotational) pole.
• Assuming that the two hour angles aren't the same, an API question comes up, which I suppose I can go ahead and start working to decide even before I have numbers: how should I distinguish “HORIZONS hour angle” and “MKB hour angle”? A single method with a switch opting for one or two calculation methods? Or two separate methods? Once I have both approaches coded and sitting next to each other, it may become apparent which API approach is best.

[mkbrewer]

• I realize that this level of error is insignificant. I was just pointing out that there is a slight problem with the math there is all. The fact that the error increases smoothly indicates that the problem is with the starting and ending values rather than the interpolation itself.
• I'm not sure why you'd want to do this. The only reason that I haven't added it to my code is that I don't care much about sub mas errors so long as I know what is causing them.
• JPL adds deflection by the Earth also. There is a discrepancy in the cutoff times though. They cut off the deflection when the Moon is at -5 degrees elevation at moonset, whereas you cut it off when the Moon is at -19 degrees elevation. At moonrise, they cut it off when the Moon is at -27 degrees elevation, whereas you cut it off when the Moon is at -19 degrees elevation again.

• That's easy. JPL's hour angle is just LST - RA; the TIRS hour angle. Here's a plot of the difference between Skyfield's TIRS hour angle and JPL's compensated for the difference in LST (averaging 7.378 ms).

• Two points here: It's not MKB's hour angle, it's SOFA's "observed" hour angle (the proper ITRS hour angle that should be used with an equatorial mount aligned with the geographic pole) and I'm not asking for the TIRS hour angle at all. That is trivial to calculate given the apparent RA and LST that Skyfield already provides.

[mkbrewer]

Oh, and here's the JPL Horizons data if you want to try this out yourself.

moon_topo_4_6_2017_jpl_horizons.csv.txt

[brandon-rhodes]

I can vary the implementation later (for now I directly use the known longitude of the observer), but hadec() is now available and should work if anyone wants to install Skyfield from the development branch:

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

If the numbers look good, we can close this issue, and I can separately track the possible additional to-do items mentioned above.

Oh, and, to respond to the comments about them:

• Yes, the nanoseconds are not important, but each time I see the dangling decimals again after a few weeks away from them, they make me double-check whether I used the wrong timescale for interpolating. Skyfield goes to a bit of trouble in several places to make sure that if you put a number in and ask for it later, you get the same number back out, and I won't mind seeing if that could be accomplished for polar motion numbers.
• (a) Because those jagged graphs annoy me. (b) Because maybe the second most common question I get about the library is "why doesn't it match" some other library or service, and adding a switch for boutique adjustments like Earth deflection would let me write some docs in the future that show "here are all the switches to try throwing if you're trying to get Skyfield to agree with X" and maybe it would help users work out the differences themselves without me in the loop.
• Gah! Asymmetrical cutoffs at –5 and –27!? I wonder if that's deliberate, or the result of a HORIZONS cutoff that pays attention to something other than zenith distance.
• Thanks for the additional notes about the HORIZONS values; they'll be helpful if I ever try to match their HA in Skyfield. But for now I'm going with the SOFA concept, which matches the Explanatory Supplement’s assertion that polar motion affects HA at high precision. (That is: if I'm understanding the Supplement, they agree that polar motion not only affects slightly the xyz coordinate of an observer specified by lat/lon, but also affects the HA angle directly by tilting the equator along which HA is measured.)

[mkbrewer]

I've been hesitant to try installing from the development branch because I'm unsure what that will do. Will it write over my current installation of Skyfield without changing the version number? Or will it install a separate version alongside my currently installed version?

A couple of comments on your implementation:

• First of all, rotating by GAST will rotate from the equinox to the TIO, so you need to add the TIO locator to the geographic longitude there. Another way to think about this is that LST contains the TIO locator, so HA must also so that the ITRS "observed" RA = LST - HA is independent of the TIO locator.
• Second, this implementation is only valid in an airless model as it doesn't include refraction. An equatorial mount is just an AltAz mount that has had its azimuth axis rotated North-South by 90 degrees minus latitude (an equatorial mount located at the pole would be completely equivalent to an AltAz mount). This is why the proper pointing for an equatorial mount is just a rotation of Alt-Az.

[brandon-rhodes]

I've been hesitant to try installing from the development branch because I'm unsure what that will do. Will it write over my current installation of Skyfield without changing the version number?

Yes, that's what it would do. It removes and replaces the current install but the version number will be the current one from PyPI, as I don't increment it til it's time to release (lest I accidentally release early). If you're wary I can close the issue on the strength of the passing test, then you can re-open if you see problems after the next version comes out.

First of all, rotating by GAST will rotate from the equinox to the TIO, so you need to add the TIO locator to the geographic longitude there. Another way to think about this is that LST contains the TIO locator, so HA must also so that the ITRS "observed" RA = LST - HA is independent of the TIO locator.

The raw position vector is topocentric GCRS. The geographic longitude is an ITRS angle, and thus is over on the other side of the transform that includes the TIO locator. So when I apply the ITRS transform to the topocentric GCRS, which includes the TIO locator, aren’t I generating an angle that’s directly comparable to the geographic longitude? Your suggestion sounds like it would apply the TIO locator twice, but maybe I don't see where exactly it would go in the code.

Second, this implementation is only valid in an airless model as it doesn't include refraction.

Understood: at the moment my target is the Hour Angle described in the Supplement, which accounts for polar motion but not for atmosphere. The code will indeed have to become more complicated, as shown in your script, if in the future we go beyond the almanac concept of HA.

[mkbrewer]

Sorry to be long in responding. I saw my error immediately, but I got tangled up in signs until I realized that you use alibi transformations in your rotations so that your wobble matrix subtracts the TIO locator from the sublongtiude rather than adding it. So, yes, you're fine there as is.

[brandon-rhodes]

Sorry to be long in responding.

No problem! It’s the weekend; and, as you are my only correspondent who regularly checks my work at this level of precision, your comments are valuable whether they take a day or a week.

I'm going to close this issue out pending a release soon, and it can be re-opened if the new HA routine winds up breaking for someone. Thanks for the input, everyone!

[mkbrewer]

I tried installing the development branch using your instructions, but it didn't work. After downloading 93.6 MB, I got:

Requirement already satisfied (use --upgrade to upgrade): skyfield==1.35 from https://github.com/skyfielders/python-skyfield/archive/master.zip in /usr/local/lib/python2.7/dist-packages

and it did nothing.

[mkbrewer]

You should probably document the fact that your ITRS hour angle and declination don't include refraction.

[brandon-rhodes]

I tried installing the development branch using your instructions, but it didn't work. After downloading 93.6 MB, I got…

I'll adjust the install-from-a-GitHub-zip template I use to add -U so that it doesn't skip an install of what looks like the same version.

And, yes, I guess the zip includes all the files in the repository, even the ones that aren't needed for the install. I'll think about whether there's a way around that.

You should probably document the fact that your ITRS hour angle and declination don't include refraction.

Done.