sbailey / surveyqa

DESI Survey QA Dashboard
BSD 3-Clause "New" or "Revised" License
0 stars 0 forks source link

Add hour-angle histograms and timeseries #30

Closed sbailey closed 5 years ago

sbailey commented 5 years ago

Add code for calculating the "hour angle" of each exposure, which is basically a measure of how far off you are from observing each tile as it passes through the meridian going directly overhead from south to north. Background:

The "local sidereal time" (LST) is a measure of time relative to the distant stars instead of to the sun, i.e. LST depends upon the longitude of where you are on earth and what is overhead right now. The hour angle of a particular exposure = LST - RA, where RA is the Right Ascension of wherever the telescope is pointing.

Detailed routines are available in astropy, but the approximations at https://aa.usno.navy.mil/faq/docs/GAST.php are good enough for here. Converting units from hours to degrees, using MJD instead of JD, and using the longitude of Kitt Peak, I get:

D = MJD - 51544.5
LST = (168.86072948111115 + 360.98564736628623 * D) % 360

Please double check me on those numbers, cross checking to http://www-kpno.kpno.noao.edu/Info/Mtn_Weather/ while noting that MJD = JD - 2400000.5 .

Calculate the hour angle for every exposure and include histograms of that on the summary page and timeseries and histograms on the nightly pages. Coordinate with @ana-lyons and issue #24 for adding placeholders for those plots while @williyamshoe works on adding the calculations to enable those plots.

Watch out for the RA and LST wraparound around 0/360. The hour angle of a DESI exposure should always be between [-60, 60], not something like 359.

Old-school "hour angles" are measured in hours (thus the name...) but let's do everything in degrees.

williyamshoe commented 5 years ago

I double checked the calculations by comparing the derived LST from the MJD on the website and the real LST found on the website. I get approximately the same answer. I also derived an equation that converted D into GMST then into LST, and got approximately the same equation as yours.

However, when I calculate the hour-angles for each exposure, they seem to range from [-60, 80], not [-60, 60]. However, it is worth noting that only a small fraction of exposures are above 60.

This is the code I have for adding the hour-angles. Is there anything wrong with it?

D = exposures["MJD"] - 51544.5
LST = (168.86072948111115 + 360.98564736628623 * D) % 360
exposures["HOURANGLE"] = LST - exposures["RA"]

def change_range(i):
    if i > 180:
        return i - 360
    if i < -180:
        return 360 + i
    return i

exposures["HOURANGLE"] = [change_range(i) for i in exposures["HOURANGLE"]]
sbailey commented 5 years ago

Your code looks right, though the change_range logic could be implemented more efficiently by taking advantage of array vector logic:

exposures["HOURANGLE"] = LST - exposures["RA"]
ii = exposures["HOURANGLE"] > 180
exposures[ii] -= 360
jj = exposures["HOURANGLE"] < -180
exposures[jj] += 360

or more cryptically (and efficiently) in a single step, while including a comment to define intention:

exposures["HOURANGLE"] = LST - exposures["RA"]
# Map [-360,360] into [-180,180]
exposures["HOURANGLE"] = (exposures["HOURANGLE"]+180) % 360 - 180

Please track down specific exposures that have abs(HA) > 60 and post their NIGHT and EXPID here, and track down by hand (looking at their DATE-OBS and RA and using an independent LST calculator like http://www.csgnetwork.com/siderealjuliantimecalc.html). An HA of 80 degrees is crazy large (90 would be looking at the horizon) so if it is real we have a bug in our field selection. Also note: filter only on science exposures: the calib exposures are pointing at a screen inside the dome and their HA doesn't matter.

williyamshoe commented 5 years ago

I get a total of 39 exposures that are over 60 degrees for the HA (none are below -60 degrees). They are sorted by HA in the table below. These are all non calib exposures.

EXPID NIGHT HOURANGLE
2722 20200218 60.0505073587
27 20191203 60.2752482903
2127 20200130 60.35577338
2031 20200123 60.708569284
7154 20200611 60.8413247381
7155 20200611 61.3426936891
6180 20200511 61.7480204607
6181 20200511 62.0404856841
28 20191203 62.6149700426
6965 20200531 62.6325506459
7182 20200612 63.1596488855
4214 20200322 63.4624798289
7183 20200612 64.6637557279
4215 20200322 64.7159021986
3652 20200312 65.5369523395
6963 20200531 65.5950214309
1531 20200114 65.7818533822
110 20191205 65.890377435
3653 20200312 65.9129790537
60 20191204 66.0071841251
4216 20200322 66.303570535
1532 20200114 66.6174682972
6964 20200531 66.8484438006
111 20191205 68.4390029144
6966 20200531 68.4400742732
4362 20200325 69.1232561523
6049 20200502 69.3403809386
6050 20200502 69.9253113805
4363 20200325 70.1259940493
7116 20200610 70.5764954756
6413 20200516 70.6128681876
61 20191204 70.9373120974
7117 20200610 71.1614259174
6414 20200516 71.4902638482
112 20191205 71.5462594512
6967 20200531 73.2866407547
114 20191205 73.4667992311
113 20191205 73.8024197128
115 20191205 75.0126868191
sbailey commented 5 years ago

Thanks for the examples. I forgot a cos(dec) term in my expectations. Plotting hour angle is still useful, but the constraint is we should stay within |HA|*cos(dec)<60, so northern portions of the footprint with cos(dec)<1 can have HA>60 . Spot checking a few cases, it looks like that is what is going on.

Please proceed with making the HA plots and submitting a PR.

Context for the future: at times of night when the galaxy is directly overhead (LST ~ 90 and 300) we don't have any tiles to observe at HA=0, so part of the survey planning process is to assign every tile a design desired HA (or equivalently desired LST). Our current inputs those have those design HA assignments in them, but if we get that into the input files, we may replace these HA plots with HA-designHA instead.

sbailey commented 5 years ago

HA histogram was added to summary page in PR #37; remaining item for this issue is to also fill in the HA vs. time and histogram on the nightly pages. @williyamshoe can you add that (even though normally @rdoshi99 is the nightly plot person)? Thanks.

williyamshoe commented 5 years ago

Alright, just added the hourangle timeseries and histogram, replacing the current placeholders. Check out #40 (also just realized I meant to name the branch as hourangle_nightly_branch, not hourangle_summary_branch).