Fixed Az computation.

[mlandriau]

The (Ra,Dec) to (Alt, Az) computation had a bug when the Az was in some quadrant. Now fixed and also corrected the standard convention with Az = [0, 360[ .

[sbailey]

Please write some unit tests to check this; I'm getting some NaNs in output, e.g. desisurvey.utils.radec2altaz(ra=0, dec=32.96, lst=0.0) -> (89.003972222222345, nan)

[dkirkby]

I am in favor of eventually using the routines in astropy.coordinates for these types of low-level operations, since those are well vetted and documented, and then we can focus on the DESI-specific stuff in this package.

However, if this fixes a bug then please add a unit test to demonstrate correct behavior and lets merge and push ahead.

[mlandriau]

The last commit was to fix the bug I was looking for when I uncovered the one in the ra,dec to alt,az conversion. @sbailey does not understand why the code as it was was giving me a rubbish answer instead of throwing an error.

@dkirkby Not sure when I wrote this function, but I did remove some time ago all calls to astropy functions that involve LST.

[sbailey]

I suspect that the avoidance of LST-related calculations in astropy was motivated by an early version of astropy that didn't know how to do LST and AltAz transforms for future dates because the IERS earth axis wobble tables only covered the past. i.e. formerly useless for DESI survey sims. Recent versions now print a warning but at least proceed:

WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar
motion from the 50-yr mean for those. This may affect precision at the 10s of arcsec level
[astropy.coordinates.builtin_frames.utils]

For the purposes of figuring out the LST or moon avoidance of a tile, having a simplified function like this is fine, but we definitely do need tests that it is correct.

[mlandriau]

Instead of the 50-year mean, I implemented a formula that is in the preamble of the IERS bulletin, which is not as accurate as their predictions, but can be used to extend the range of the table (according to the Naval Observatory).

[dkirkby]

I am ok to merge this now (with a new unit test) but I would rather let the astropy folks worry about IERS, etc, in future. (That warning message is quite informative, and can easily be muted.)

[mlandriau]

I added a unit test. It tests 3 values against the ones from astropy. The interface is clunky and not very transparent, so I'm not sure how to get the results I want, e.g. the way I ran it, it gets corner cases ra=lst wrong (well correct within a couple of degrees), while my code gets them right.

[dkirkby]

A couple of degrees discrepancy sounds scary. Can you give an example of an astropy calculation that is giving the wrong answer?

[mlandriau]

I was exaggerating a bit when I said a couple. One of the unit tests has an azimuth at transit of 0.4deg and it should be 0. Another example, which shows the largest discrepancy I've seen: ra=lst=15*6.116476882227936 , dec=20 gives alt, az = 78.0288521, 178.84983768 whereas my code gives 78.036027777777804, 180.0. So the astropy result is 1.15deg of for the azimuth (the differences are always largest for the azimuth), which is too large to be accounted for by the precession of the equinoxes.

[dkirkby]

How exactly did you calculate the astropy values? (I just see hard-coded numbers in your unit test).

[mlandriau]

I fiddled a lot on the command line and copied and pasted the values I used for the unit test. Going through my terminal, I see that this is what I ended up using:

---

from datetime import datetime from astropy.coordinates import EarthLocation, SkyCoord, AltAz from astropy import units as u from astropy.time import Time

x = EarthLocation.of_site('Kitt Peak National Observatory') t = datetime(2017,4,15,0,0,0) tt = Time(t, location=x) tt.sidereal_time('mean') # Need this info to input into my function

y = SkyCoord(ra=0.0u.degree, dec=45.0u.degree) y.transform_to(AltAz(obstime=t, location=x))

Note that the conversion does not allow to input the time as LST, which is what we need for the afternoon planning.

On Thu, Apr 6, 2017 at 1:59 PM, David Kirkby notifications@github.com wrote:

How exactly did you calculate the astropy values? (I just see hard-coded numbers in your unit test).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/desihub/desisurvey/pull/31#issuecomment-292318480, or mute the thread https://github.com/notifications/unsubscribe-auth/ASznwtmiRSvFpAJpG77cUFYmSCw7-xyyks5rtVJIgaJpZM4Mzg6K .

-- Martin Landriau Lawrence Berkeley National Laboratory 1 Cyclotron Road Mailstop 50R5008 Berkeley, CA 94720 USA

Tel: +1 510 486 7934 Fax: +1 510 486 6738

[weaverba137]

If no timezone information is supplied in the datetime object, are you quite certain that the time that astropy.time.Time and astropy.coordinates.AltAz uses is the time you think it is?

[mlandriau]

The differences I'm seeing are of the order 0.1-1deg in Az and 0.01-0.1 in Alt; an error of just one timezone would produce errors of ~15deg. The difference between mean and apparent sidereal time is of the order 0.001deg.

On Thu, Apr 6, 2017 at 2:47 PM, Benjamin Alan Weaver < notifications@github.com> wrote:

If no timezone information is supplied in the datetime object, are you quite certain that the time that astropy.time.Time and astropy.coordinates.AltAz uses is the time you think it is?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/desihub/desisurvey/pull/31#issuecomment-292333869, or mute the thread https://github.com/notifications/unsubscribe-auth/ASznwghHtQl2Hv2yXGLlsaCMJ3dfvXuxks5rtV18gaJpZM4Mzg6K .

-- Martin Landriau Lawrence Berkeley National Laboratory 1 Cyclotron Road Mailstop 50R5008 Berkeley, CA 94720 USA

Tel: +1 510 486 7934 Fax: +1 510 486 6738

[dkirkby]

Note that the conversion does not allow to input the time as LST, which is what we need for the afternoon planning.

Isn't the LST currently calculated from an MJD, so we could use the MJD directly?

tt.sidereal_time('mean') # Need this info to input into my function

I think this should be 'apparent', not 'mean', since that's what desisurvey.utils.mjd2lst() uses.

I suspect the discrepancies between the routines in desisurvey.utils and astropy.coordinates are coming from your earthOrientation() implementation, which seems to introduce an offset of dLST ~ 0.57 deg:

[dkirkby]

For reference, the plot above was created with:

def plot_mjd2lst(mjd1=157300, mjd2=157800, n=500):
    mjd = np.linspace(mjd1, mjd2, n)
    tvec = Time(mjd, format='mjd', location=kpno)
    lst_mean = tvec.sidereal_time('mean').to(u.deg).value
    lst_app = tvec.sidereal_time('apparent').to(u.deg).value
    lst = np.empty(n, float)
    for i, t in enumerate(tvec):
        lst[i] = mjd2lst(t.mjd)
    dlst_mean = np.fmod(lst_mean - lst + 360, 360)
    dlst_app = np.fmod(lst_app - lst + 360, 360)
    plt.plot(mjd - mjd1, dlst_app, 'b-', label='astropy(app) - mjd2lst')
    plt.plot(mjd - mjd1, dlst_mean, 'r-', label='astropy(mean) - mjd2lst')
    plt.xlabel('MJD - {0}'.format(mjd1))
    plt.ylabel('dLST [deg]')
    plt.legend(loc='upper left')

[mlandriau]

From what Stephen said in a previous emails, astropy defaults to a constant value, while the function uses an approximate function (from the Naval Obs). I remember checking the values were consistent with the ser7 tables (according to the Naval Obs, it can't reproduce exactly the tables, but can be used to extend it), but I can re-verify.

But that doesn't matter here because I don't use this function to compute LSTs in the examples I gave above. I used the one that astropy produces so that my equatorial to horizontal function and theirs have exactly the same input. The input to their function is: ra, dec, datetime object, earth location object corresponding to KPNO; the input to mine are the same ra and dec, the lst computed form astropy using the same datetime and earth location object and - implicitely - the latitude of KPNO (which I've just checked is different from what astropy has by 2.4", which corresponds to 74m, smaller than the extent of the site).

What I was trying to figure out when obtaining numbers with which to write the unit test is the following: is there a series of input that will make the astropy function give the correct answer in the case lst=ra? If yes, what is it? If not, why?

Martin

On Thu, Apr 6, 2017 at 4:37 PM, David Kirkby notifications@github.com wrote:

For reference, the plot above was created with:

def plot_mjd2lst(mjd1=157300, mjd2=157800, n=500): mjd = np.linspace(mjd1, mjd2, n) tvec = Time(mjd, format='mjd', location=kpno) lst_mean = tvec.sidereal_time('mean').to(u.deg).value lst_app = tvec.sidereal_time('apparent').to(u.deg).value lst = np.empty(n, float) for i, t in enumerate(tvec): lst[i] = mjd2lst(t.mjd) dlst_mean = np.fmod(lst_mean - lst + 360, 360) dlst_app = np.fmod(lst_app - lst + 360, 360) plt.plot(mjd - mjd1, dlst_app, 'b-', label='astropy(app) - mjd2lst') plt.plot(mjd - mjd1, dlst_mean, 'r-', label='astropy(mean) - mjd2lst') plt.xlabel('MJD - {0}'.format(mjd1)) plt.ylabel('dLST [deg]') plt.legend(loc='upper left')

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/desihub/desisurvey/pull/31#issuecomment-292362010, or mute the thread https://github.com/notifications/unsubscribe-auth/ASznwlJ1u6F_hIYT0eP_V7ZUgfeyI4wKks5rtXcigaJpZM4Mzg6K .

-- Martin Landriau Lawrence Berkeley National Laboratory 1 Cyclotron Road Mailstop 50R5008 Berkeley, CA 94720 USA

Tel: +1 510 486 7934 Fax: +1 510 486 6738

[dkirkby]

I just noticed that the MJD range in my plot is crazy (mjd1=157300, mjd2=157800). With more reasonable values (mjd1=57300, mjd2=57800) I get:

The corresponding date range is 2015-10-05 to 2017-02-16, so the discontinuity is probably where the IERS table ends.

[mlandriau]

After some more investigating, I found that the discrepancy is due, in large part, to the precession of the equinoxes.

The code below illustrates the problem I've encountered with Astropy when writing the unit test for desisurvey.utils.radec2altaz(ra, dec, lst). The difference in Alt for this sample is smaller than 0.01deg, except for one point nearest to zenith where it is 0.018deg; this is also the point where the difference in azimuth is largest (6.9deg!!!).

However if I use noon on 1 January 2000 in stead of the date in the code below, the errors in both Alt and Az are about 0.005deg, except near zenith when the error in Az is 0.03deg. Still not exactly 180.0 and 0, but better. It does seem to point that the astropy solution is unstable near zenith.

from datetime import datetime
from astropy.coordinates import EarthLocation, SkyCoord, AltAz
from astropy import units as u
from astropy.time import Time
from desisurvey.utils import radec2altaz

x = EarthLocation.of_site('Kitt Peak National Observatory')
t = datetime(2017,4,15,0,0,0)
tt = Time(t, location=x)
tsid = tt.sidereal_time('mean')
lst = 15.0 * float(tsid/u.hourangle) * u.degree

print("Alt, Az at transit for various declinations.")
dec0 = -90.0
while dec0 <= 90.0:
    y = SkyCoord(ra=lst, dec=dec0*u.degree)
    z = y.transform_to(AltAz(obstime=t, location=x))
    az = (z.az)/u.degree
    alt = (z.alt)/u.degree
    print("Dec =", dec0, "degrees: Alt =", alt, "degrees, Az =", az, "degrees.")
    alt2, az2 = radec2altaz(float(lst/u.degree), dec0, float(lst/u.degree))
    print("                          ", alt2-alt, "        ", az2-az)
    dec0 += 10.0

[dkirkby]

the errors in both Alt and Az are about 0.005deg, except near zenith when the error in Az is 0.03deg. Still not exactly 180.0 and 0, but better. It does seem to point that the astropy solution is unstable near zenith.

The Az solution is inherently ill defined near zenith so dAz is not a good metric. Instead, I would compare the opening angle between (Alt, Az) vectors from different methods.

The effect of polar motion on (Alt,Az) is real, so needs to be included in the transforms used by desisurvey, but cannot be calculated from the (Ra,Dec,LST) passed to desisurvey.utils.radec2altaz().

I still think we need to move towards astropy functions (and I haven't seen evidence yet that they are insufficiently accurate), but lets defer that to a separate issue. Since this PR fixes a known bug, I am in favor of merging now and pushing ahead on other fronts.

[sbailey]

OK, this is dragging on too long; we can merge now while deferring other updates to #32. I still do want some sort of sanity check unit tests though, e.g. this function originally was returning NaNs for some cases that unit tests would have caught...

[mlandriau]

I committed a change to the unit test after this was merged (only just saw the message). The NaNs were a result of an out of range argument for arccos which was corrected as soon as the problem was reported.

The rationale behind this is to have a simple interface which allows the computation for arbitrary LSTs for the HA assignment stage of the survey planning. The airmass for each tile is computed for each LST bin in advance.

On Sat, Apr 8, 2017 at 4:13 PM, sbailey notifications@github.com wrote:

Merged #31 https://github.com/desihub/desisurvey/pull/31.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/desihub/desisurvey/pull/31#event-1035306158, or mute the thread https://github.com/notifications/unsubscribe-auth/ASznwqXJP15BdBkqmsNUhi4po6WXINHxks5ruBSOgaJpZM4Mzg6K .

-- Martin Landriau Lawrence Berkeley National Laboratory 1 Cyclotron Road Mailstop 50R5008 Berkeley, CA 94720 USA

Tel: +1 510 486 7934 Fax: +1 510 486 6738