Repository navigation
CIRS is not at the ultimate accuracy #10356
Description
Activity
Now that GCRS is pretty much all set, it's time to address the other elephant in the room which are the GCRS to CIRS and CIRS to observed transforms. The problem here is that instead of following the
atcochain, you chose to useatciandatioinstead. What SOFA fails to mention, is that the later is completely incompatible with solar system bodies where diurnal parallax comes into play (although, if you look closely, the documentation foratci13()does say that it is specifically for star data) This is especially true of the Moon, so that is what I've concentrated on here. There, I can't help but notice a rather large (to my standards) error in the returned Alt/Az coordinates. Here is the error between your Alt/Az and mine in the latest data file that I have from you.
This is in coordinate space and on this date, the Moon reaches lower meridian transit at 13:05 UT. At that time the elevation is -80.5 degrees, so I also plotted the actual offset on the sky which looks better, but is still pretty bad.
I was able to reproduce this exactly (to 20 uas) by altering my program to follow the steps that you use.- Transform the topocentric GCRS coordinates back to BCRS..
- Drop the EarthLocation and transform to geocentric GCRS.
- Apply the BPN matrix to transform to geocentric CIRS.
- Add the EarthLocation back in to transform to topocentric CIRS.
- Use
apio()to set up the CIRS to observed transform. - Call
atioq()to transform to Alt/Az.
It turns out that the root of the problem is in the transform to geocentric GCRS. Specifically the geocentric aberration. When observing the Moon, there is a contribution to diurnal aberration from the diurnal parallax. This comes in to the aberration calculated by
ab()due to the diurnal offsets in the topocentric direction vector from the observer to the Moon. I was able to confirm this by altering the steps that you use. In the transformation to geocentric GCRS, I retained the topocentric BCRS unit vectors, but calledapcg()instead ofapco()to tell SOFA that I was using geocentric positions even though I wasn't. This eliminated the velocity based component of the diurnal aberration but retained the position based component. Then comparing the resulting CIRS Ra/Dec to that obtained by using the full topocentric GCRS transform yielded a nice velocity based diurnal aberration. Very sinusoidal with the Declination component dropping to zero at upper and lower meridian transit.
Then I usedapio()as before to set up the velocity based diurnal aberration during the transformation to Alt/Az. The resulting Alt/Az positions were in agreement with my normal data (full topocentric transform usingapco()) to within 15 uas.
On the other hand, the CIRS RA/Dec offsets obtained from following your procedure exactly weren't nearly as pretty.
There is a phase shift and an offset in the Declination sinusoid due to the missing position based diurnal aberration.Clearly, the way to solve this problem is to merely apply the BPN matrix directly to yield topocentric CIRS coordinates without any folderol of transforming the GCRS coordinates to geocentric first. Then you can still use
apio13(), if you'd like, to set up the Alt/Az transformation, but you must be sure to setastrom['diurab']to zero before callingatioq(). Since this is rather inefficient, it would better to make a modified version ofapio()dropping the call topvtob()and settingastrom->diurabto zero directly.Note that you can also apply the BPN to the
gcrs_frame.obsgeolocattribute to get the appropriatecirs_frame.obsgeolocattribute to apply to topocentrc CIRS (topocentric intermediate place). Yes, I know that, as currently implemented, CIRS doesn't have anobsgeolocattribute, but it should. Also note that the true equator true equinox (TETE) frame that I requested (issue #10230) can also be realized in the same manner by merely applying the equinox based BPN to either geocentric GCRS (proper place) or topocentric GCRS (local place) to derive geocentric TETE (apparent place) or topocentric TETE (topocentric apparent place)@maxnoe - it may be worth holding your trigger on this for a few days. I'm about to push an implementation of apparent place coordinates (TETE) which, when finished would make a reasonable starting point for how to implement CIRS with obsgeoloc/vel.
Most of the methods I use could be ported to the corrections in CIRS quite simply, but the ICRS<->CIRS transforms will need fixing, or indeed removing, since in principle one can rely on ICRS->GCRS->CIRS without the direct transformation.
I second what @StuartLittlefair has said here. The only difference between CIRS and TETE is the BPN matrix that is used. If topocentric intermediate place inherits its observer location from the GCRS frame, then the observer's geodetic longitude and latitude can easily be extracted for use in conversion to AltAz or HADec.
But to answer your question, the idea here would be to start with a GCRS frame (either geocentric or topocentric), then call
get_cip()andc2ixys()to build the CIO based BPN matrix for transformation to CIRS. Or equivalently, callgcrs_to_cirs_mat()fromintermediate_rotation_transformsto get the BPN matrix directly. Please note that, as currently implemented,gcrs_to_cirs()fromintermediate_rotation_transformswon't work for this as it converts from topocentric to geocentric.On the other hand, just like TETE, one can envision an AltAz or HADec frame with it's unique time and observer location as is currently done. But then the CIRS location would need to be converted to GCRS in order to do the transformation back to ICRS/BCRS to get the source position at the given time and then pass back through GCRS to finally arrive at CIRS ready for conversion to AltAz or HADec. .
Seems the long way around. Of course,
apco()will do that for you. It would only be necessary to extract the GCRS position of the observer thatapco()constructs to derive the ICRS/BCRS position of the source, then use the resultingastromdata structure for conversion to observed coordinates. Seeapco13()for an example of how to do this.I'm getting somewhere with this now. I've added a topocentric position to CIRS (in GCRS coordinates for consistency with other frames). I think I have the ICRS<->CIRS and GCRS<->CIRS figure out. Just have to try and fix CIRS<->AltAz following @mkbrewer's suggestions
Reacted by Marten van KerkwijkWell, partial success; I have managed to reduce the discrepancies by following the two routes @mkbrewer suggested - either using
apio13and settingastrom['diurab']=0, or usingapcoto get theastromcontext. However, the discrepancies still show some residual structure for @mkbrewer's moon test case.
Note: I get identical results using the direct ICRS->AltAz transform sketched out in #10887
For some reason I am unable to get a good comparison with Horizons, or SkyField's output. One thing that makes this difficult is that I am not sure I have an up-to-date file of the output from @mkbrewer's script. I'm pretty sure the version I have is out-of-date. Do you think you could share this here?
I'd like the correct Az and El for the Moon on Apr 6 2017 from an observatory at 22.959748 deg S, 67.787260 deg W and 5186m...
What version of the output data that I gave you are you using here? If it's
moon_topo_4_6_2017_mkb_v5_ap_polar_motion.csv, then I'd like to see the Ra/Dec comparison also. Perhaps the best thing to do is for you to post your output so that I can look into the apparent discrepancy.It’s a lot older than that - the file name is just
moon_topo_4_6_2017_mkb_ap.csv. If I have a copy of the latest data then I know if the things I am doing are working.My file from my efforts so far is attached. The fact that I get the same answer from the snippet in #10887 must be telling me something, but I don't know what yet!
Here you go.
5 remaining items
Note that you can also apply the BPN to the gcrs_frame.obsgeoloc attribute to get the appropriate cirs_frame.obsgeoloc attribute to apply to topocentrc CIRS (topocentric intermediate place).
#10356 (comment). The reason that I suggested using the observer's CIRS location in the CIRS frame rather than the GCRS
obsgeolocis that it avoids the necessity of recomputing the GCRS->CIRS BPN matrix as you already have that available during the transform from GCRS to CIRS.In the above case, this BPN matrix is actually computed twice. Once in
EarthLocation.get_gcrs_posvel()to get the observer'sobsgeolocin the GCRS and a second time to convert this to the observer's location in the CIRS. This could be obviated by simply creating a functionEarthLocation.get_cirs_pos().OK - so having now convinced myself that astropy is doing the right thing when converting geocentric CIRS to topocentric CIRS I am almost ready to add some tests to my branch and submit a PR. I have one outstanding question for @mkbrewer. I find excellent agreement with your data (
moon_topo_4_6_2017_mkb_ap_v5_10356.txt), but poor agreement (up to 1" out) with JPL Horizons and Skyfield.
The agreement gets much better if I turn polar motion off:
Although interestingly there I recover the patterns that I saw compared to an older output from your code (
moon_topo_4_6_2017_mkb_ap.csv).Comparing HORIZONS to Skyfield reveals a similar pattern so it seems almost everyone disagrees at this level apart from us!
I am surprised that HORIZONS seems to ignore polar motion by default, but my question is whether you understand the pattern that remains after I turn polar motion off?
@StuartLittlefair Sorry to be long in responding. I didn't notice this post for awhile and then I had issues with my installation of Skyfield that turned out to be due to it no longer being able to download EOP data from the USNO. After upgrading to the latest version, it's working again. If I disable polar motion in my program and Skyfield, I get very good agreement.
Here is the script that I used:
Regarding JPL Horizons, that's a whole can of worms there. I have never been able to get good agreement in AltAz with them. They do use polar motion, but the EOPs that they use are significantly different than the IERS published data. You can request the EOPs that they use if you wish. I did that once and I still couldn't get good agreement. This led me to compare my output to Skyfield, but unfortunately, there are rather significant issues that Skyfield has with polar motion. I've been trying to to get Brandon Rhodes to address this since May, so far to little avail. In fact, he recently downgraded the issue I opened to a feature request rather than a bug, so I have no idea when he's going to get around to addressing it.
Thanks @mkbrewer - don't worry about the time it takes you to reply. I'm asleep whilst you are working! Upgrading Skyfield to the latest version made a big difference for me, but I'm not quite as close as you still.
The output I have from your code is
moon_topo_4_6_2017_mkb_ap_v5_10356.csv. Is the only difference between the output you compare to Skyfield and the output I have polar motion? In which case I probably haven't disabled it properly in astropy...
Nevertheless we are looking good to go here - I just need to write some tests. This is why I was trying to compare to HORIZONS - it would be good to have an independent third party!
Yes, I simply set xp and yp to zero when I call
apco().btw, upgrading Skyfield made a difference for me also. I am no longer seeing the oscillations that were previously present. This is most easily seen in the residual for the distance to the Moon.
No, using that code I get the same as Skyfield to better than 0.6 mas.
OK. I don't know how you are disabling polar motion in the standard transform. Is there any possibility that you are disabling it in some places, but not others? One that comes to mind is
get_gcrs_posvel().That was my thought as well. I’m not too concerned about chasing it up because I don’t think we’d be interested in turning it on/off.
Well, If I can ever get Brandon Rhodes to fix polar motion in Skyfield, this won't be an issue anyway. On that, once this is finished, it will be quite helpful as he can compare his results to yours without having to rely on data that I give him.
@StuartLittlefair I am pleased to announce that polar motion is now fixed in Skyfield. Here are the new residual plots between my program and Skyfield.
This fix should be available at the next release of Skyfield.
Copying @mkbrewer's #10230 (comment):
It sounds like we could improve things fairly simply since GCRS seems OK: