astropy / astroplan

Observation planning package for astronomers – maintainer @bmorris3
https://astroplan.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
198 stars 109 forks source link

Potential Sunrise/Sunset Discrepancy #575

Closed prappleizer closed 10 months ago

prappleizer commented 10 months ago

Astroplan documentation claims that 150 grid points should be accurate to within ~1 minute. I have been playing around with it, and at least for this week at Keck Observatory, the sunset and sunrise times appear to be discrepant with the sunset time provided by Keck and other sources by ~10 minutes. I have tried running with n_grid_points = 1000 (should be highly accurate) and the discrepancy still exists. There are different ways of defining the sunset moments --- but both a horizon choice of 0 deg and -0.8333 deg are discrepant (0 more discrepant). Here is the code I've run

k = Observer.at_site('Keck', pressure=0.6/1.03*u.bar)
time = Time('2023-10-14 00:00:00')+10*u.hr # in utc
astroplan_sunset = k.sun_set_time(time=time,horizon=-0.8333*u.deg, which='previous',n_grid_points=1000)

(I've tried different choices of atmospheric pressure but none make a major difference).

In this example, (with a reference time at midnight hawaii time) the quoted sunset time is

2023-10-13 18:00:25

On the other hand, several other sources quote discrepant times:

I am concerned because these three times are consistent with one another to within ~1-2 minutes, which is quite different than the astroplan time. (Note that the discrepancy is larger if a 0 deg horizon is used). Airmass.org uses the altitude of the observatory in question, timeanddate does not. I do not know what source Keck's own page is using.

Any advice on this? As a note, I have tried using the underlying astropy code directly to compute it myself, and I get a similar 10 minute offset:

delta_midnight = np.linspace(-12,12,300)*u.hr
obs_times = midnight+delta_midnight
frame = AltAz(obstime=obs_times,location=site)
sun = get_body('sun',obs_times).transform_to(frame)
sundown_ind, = np.where( sun.alt.to(u.deg) < -0.83*u.deg)
sundown_right = obs_times[sundown_ind[0]]
sundown_left = obs_times[sundown_ind[0]-1]
sunset = sundown_left+((sundown_right-sundown_left)/2.0) 
sunrise_left = obs_times[sundown_ind[-1]]
sunrise_right = obs_times[sundown_ind[-1]+1]
sunrise = sunrise_left + ((sunrise_right-sunrise_left)/2.0)

Where above I attempt my own grid (here of 300, but I've tried any number of grid points).

So yeah, not sure what is going on --- would appreciate any advice!

wtgee commented 10 months ago

timeanddate.com shows 18:00 for tonight (the 12th) and 17:59 for the night of the 13th for Maunakea, which is accurate.

Where do you see 18:08 for timeanddate and what Keck page are you looking at to get 18:09?

prappleizer commented 10 months ago

Keck page is here: http://mkwc.ifa.hawaii.edu/forecast/mko/

For dateandtime --- good catch in that I had it set to honolulu by accident; using mauna kea as the location, it reports 18:02 pm for tonight, so that one is actually reasonably consistent.

For airmass.org, here is a page with settings tuned to Keck on 10/13: https://airmass.org/chart/obsid:keck/date:2023-10-13/object:ngc%201052/ra:40.269994/dec:-8.255764

wtgee commented 10 months ago

airmass is doing their own thing with notes here that I unfortunately don't have the time to go through: https://airmass.org/notes

I suspect MKWC is showing the times for Honolulu (UH-Manoa) as well but we'd have to ask to find more.

wtgee commented 10 months ago

FYI, @joshwalawender also did some accurate calculations for Keck a number of years back. Hopefully things haven't regressed since then :)

https://github.com/astropy/astroplan/issues/242

prappleizer commented 10 months ago

@wtgee Thanks for this link, I hadn't seen it!

Based on perusing that issue, it looks like my discrepancy is most likely the same thing Josh found, which is that while astroplan is accounting for refraction, the apparent horizon itself is not taken into account (which Josh does explicitly in the example). I will try that out too and see how it goes!

wtgee commented 10 months ago

Closing issue, feel free to re-open if you still think there is a problem or want to post some results. Cheers.