Repository navigation
Requesting local earth rotation angle and one additional built-in frame #10239
Description
Just as I find it frustrating that the Time package provides local sidereal time, but the coordinate package doesn't currently include a built-in frame providing topocentric position with respect to the true equator and equinox of date (issue #10230), I find it equally frustrating that the Time package doesn't provide the local earth rotation angle to complement the coordinate package's built-in CIRS frame. When working with an equatorial mounted telescope both a set of equatorial coordinates and its complementary time or angle are required in order to calculate the hour angle necessary to point the telescope. The Time package should be enhanced to provide the local earth rotation angle.
In this regard, it would also be nice to provide a built-in frame to transform CIRS coordinates to observed equatorial coordinates to complement the AltAz built-in frame. This should provide both declination and hour angle (or right ascension if the local earth rotation angle is available) adjusted for refraction and polar motion. This would be very simple as the same function that is used by AltAz (atioq()) provides these also.
Activity
mkbrewer commentedon Apr 30, 2020
Observed declination and hour angle would probably be best here. These are the direct analog to azimuth and altitude and are independent of the origin of right ascension in use. It is quite easy to derive the observed right ascension in either frame from these. Subtracting the hour angle from the local sidereal time yields the observed right ascension referenced to the true equinox of date, while subtracting the hour angle from the local earth rotation angle yields the observed right ascension referenced to the CIO.
mhvk commentedon May 20, 2020
@mkbrewer - rephrasing your request just so that I'm sure I understand them correctly.
Is the request to have a method, say,
Time.earth_rotation_angle(longitude)that uses SOFA'sera00and adds the longitude? I guess, likesidereal_time, it would return the result as aLongitudewith units of hour-angle. This would be easy to do (and aftersidereal_time, which I implemented, there will hopefully be less discussion on whether this belongs onEarthLocationor onTime!). Should there be any arguments beyondlongitude(which would default to thelocationstored on the time)?And for the coordinate frame, you are thinking of an "HADec" exactly analogous to "AltAz", but for polar mount? As you note, all that would seem to take is replicating
cirs_to_altazbut catching different outputs inastropy/astropy/coordinates/builtin_frames/cirs_observed_transforms.py
Line 66 in 35085f6
The most difficult part here may be the name...
A more general comment on (2) is that between the various requests, it may be that quite a number of those frames in Fig.1 of the sofa tools docs might be useful on their own...
maxnoe commentedon May 27, 2020
I also just stumbled on this for our simulated gamma-ray events, which don't have a proper observation time and it was wished that the coordinates be not only expressed in AltAz but also in the local equatorial frame.
Having the earth_rotation_angle at
astropy.Timeseems indeed necessary, especially considering that this is "the way to go" now instead of usingsidereal_time.mkbrewer commentedon May 27, 2020
@mhvk Sorry, I somehow missed the notification of your comment. There is one additional element to this, but it is very small. That's the TIO locator given by SOFA's
Sp00(). The local earth rotation angle is then:For reference see
Apco13(),Apco(),Aper13()andAper().mhvk commentedon May 27, 2020
@mkbrewer - thanks! Just to be absolutely sure, there are no different "kinds" of earth rotation angle (like
kind='mean'orkind='apparent'for sidereal time), nor is there a need for a possibly different model? I ask since if there are such differences, I'd rather mimic the argument structure ofsidereal_time(kind, longitude=None, model=None)(though fine if there is a default for all arguments).I realize this all mostly comes from the fact that I have seen (and used) sidereal time a lot but have not encountered earth rotation angle before, and it is leaving me a bit confused that
apcotakes an earth rotation angle argument as well as ansp, i.e., with your equation, the earth rotation angle we would produce here forlongitude=0would not be what one would put as input inapco.p.s. Checking wikipedia (what else...), https://en.wikipedia.org/wiki/Sidereal_time#Earth_Rotation_Angle gives a definition just like
era00, without thesp00term, so if nothing else there will be others who will be confused as well... It also shows that my education well predates 2003...maxnoe commentedon May 27, 2020
I'd have to check with the explanatory supplements again, which are currently in my office and I am in the home office.
But ERA00 should be Earth Rotation Angle at longitude 0, so I'm also a little confused that TIO comes into play here.
mkbrewer commentedon May 27, 2020
Since the CIO is a non-rotating origin, it isn't subject to precession or nutation so there is no "mean" or "apparent" earth rotation angle. The models you refer to regarding sidereal time are just different precession/nutation models, so that isn't applicable. (see how much simpler things are when you use the ERA). The local earth rotation angle that SOFA uses is ITRS, so it is referenced to the TIO. That's why the TIO locator is included. The SOFA "observed" coordinates are also ITRS, so this is the local earth rotation angle can be used to derive the "observed" right ascension from the "observed" hour angle or vice versa, see
Atioq().See
Apco13()for the arguments that SOFA uses forApco(). For longitude=0 or'greenwich'if you will, the ITRS local earth rotation angle is just:Era00() + Sp00(). Seeastrom->along.Of course, one could define a TIRS local earth rotation angle for use with CIRS coordinates to derive the intermediate place hour angle also. This would need to add the offset described in equation (6.20) of USNO Circular 179. Or more appropriately:
which works all the way to the pole. But SOFA doesn't deal with this definition at all.
mkbrewer commentedon May 28, 2020
By the way, in case you are interested, the way that SOFA handles polar motion is rather non-transparent. What they are doing is:
Where
r_obsis the local ITRS observed minus hour angle and declination,Wis the wobble matrix andr_cirsis the CIRS right ascension and declination. This makes sense as the final rotation by geodetic longitude is applied after using the wobble matrix to transform to ITRS.However, with their definition of local earth rotation angle, what they actually do is:
Where
s'is the TIO locator. The expressionR_z(lambda_itrs) W R_z(-lambda_itrs - s')is evaluated byApco()combined withAtioq()as "Polar motion, rotated onto the local meridian".It took me awhile to figure that one out.
mhvk commentedon May 28, 2020
@mkbrewer - I think I am starting to understand (having had to wrap my head around a few other problems hasn't helped). The confusion I had certainly stems from the SOFA implementation, where in
apco13, an "earth rotation angle" is defined asastropy/cextern/erfa/apco13.c
Lines 204 to 205 in 399f7e1
and this angle is passed straight on to
apco, where it is also called "earth rotation angle" and used in a call toaperto seteralor "local earth rotation angle" inastropy/cextern/erfa/aper.c
Line 96 in 399f7e1
So, yes, I can see that that takes figuring out! But that gets to the practical side, would it be useful to be able get the local earth rotation angle also in this form (which perhaps is nameless?). This especially since you mention that "TIRS" version which we could have to complement the "ITRS" one your request is really about. In any case, since our earth locations are all ITRS, that seems definitely the sensible default, so even if expand the options, I think
longitudeshould always be the first argument, i.e., one would get a methodearth_rotation_angle(longitute=None, frame='itrs').Finally, any preference for default units?
era00itself is radians, but degrees may make more sense - at least the wikipedia article suggests that is used (rather than hours, a bit sadly, but probably wisely to avoid confusion).mkbrewer commentedon May 28, 2020
I'm afraid that I don't understand the question here. What form are you talking about?
If your intention is to also support the TIRS local earth rotation angle, it will be necessary to have an option to pass the latitude also. Of course, this is the ITRS geocentric latitude rather than the geodetic latitude which adds further complication. The reason that you can use the geodetic longitude is that in the WGS84 model the geodetic longitude is the same as the geocentric longitude. This is not the case for latitude however.
This should definitely be the same as sidereal time (u.hourangle).
mkbrewer commentedon May 28, 2020
One option is to only support TIRS local earth rotation angle if an
EarthLocationis present, which I suspect would be true 99% of the time anyway. The SOFA functionPvtob()returns the TIRS local earth rotation angle directly, albeit as a Cartesian vector giving the observer's postion wrt the CIO. The longitudinal component of that vector is, in fact, the TIRS local earth rotation angle.mkbrewer commentedon Aug 27, 2020
I thought we were on the way here, but we seem to have stalled. In case that is over confusion between ITRS (observed) and TIRS (topocentic apparent or intermediate place), here is a function that returns the TIRS longitude and an example of one way to use it in the sidereal time function. It will work exactly the same way with the local earth rotation angle. I tried to do this in Python rather than my usual C. Please forgive any errors as I am not an expert at Python.
tirs_longitude.py.txt
12 remaining items
mkbrewer commentedon Nov 22, 2020
This is why I said it would be best to use
longitude=Nonefor returning ERA, GAST or GMST. The local apparent and mean sidereal times at'greenwich'are actuallyGAST + 0 + s'andGMST + 0 + s'.mkbrewer commentedon Nov 22, 2020
I understand that this is a bit confusing. As a matter of tradition, the IAU still uses the terms Greenwich Apparent and Mean Sidereal Time even though they are no longer referenced to the Greenwich meridian.
mhvk commentedon Nov 22, 2020
It is that last point that indeed is very confusing!! But with your explanations I at least understand a bit better why.... It would have been nice if this were pointed out in the comments of the ERFA functions.
mkbrewer commentedon Nov 22, 2020
First line in
erfa.apco():mkbrewer commentedon Nov 22, 2020
If course, if the intended purpose of
sidereal_time('apparent', 'greenwich')is to return GAST, then it should stay that way. One can always usesidereal_time('apparent', 0)to get the local apparent sidereal time at the Greenwich meridian.mhvk commentedon Nov 23, 2020
No, definitely not. It is just that I took the name "Greenwich Apparent Sidereal Time" literally - it never occurred to me that "Greenwich" would not refer simply to the meridian... Now that you explained it is otherwise, I think your suggested solution of having
Nonereturn the output ofgst06directly (without offset) makes a lot of sense.mkbrewer commentedon Jan 28, 2021
@mhvk , @StuartLittlefair In version 1.36, Skyfield has added an
apparent.hadec()method that can be directly compared to the Ha/Dec output oferfa.atioq(). This comparison is good to less than 1 mas:and so can be used to test Astropy's
hadecbuilt-in frame (if it ever gets implemented).mhvk commentedon Jan 28, 2021
Nice! (And, yes, we should have the
HADecbuilt-in frame... It is always a matter of someone getting around to it...)mkbrewer commentedon Feb 27, 2021
@mhvk The new polar motion algorithm introduced in pyerfa 1.7.2 will affect the calculation of the local earth rotation angle and local sidereal time here. The change is minuscule as it is second order in the polar motion quantities.
mkbrewer commentedon Apr 30, 2021
I can implement the local earth rotation angle also, but I'd need some guidance on how to retrieve TT and UT1 from the Time object in order to do so.
maxnoe commentedon Apr 30, 2021
astropy/astropy/coordinates/builtin_frames/utils.py
Line 96 in 339bb16
mhvk commentedon Apr 30, 2021
Fantastic! In principle, as simple as
time.ttandtime.ut1, and the implementation ofsidereal_timeshould give some ideas, but perhaps best to deal with that after #11676? I'll request my own review on that one - hopefully some time next week...mkbrewer commentedon Apr 30, 2021
I'm guessing that this would be something along the lines of:
Is that correct?
mkbrewer commentedon Apr 30, 2021
My plan is to do this as a separate PR from #11676 as that will be cleaner.