Refined sky subtraction for HPF based on twilight flats

[gully]

A while back @avanderburg and I discussed how to refine the HPF sky subtraction. We discussed some PCA-related ideas among other approaches. @grzeimann recommended that we simply use the twilight flats to estimate the wavelength dependence of the target-to-sky fiber throughput. Here is a pull request to implement the latter approach, as described in the in-progress paper text:

(below reposted and abridged from a Slack conversation)

At the Exoplanet Atmospheres group meeting on Friday @astrocaroline asked about the time-stability and wavelength dependence of the target-fiber-to-sky-fiber ratio (what I define as beta in the paper). Separately @dkrolikowski asked me about the wavelength dependence. Caroline expressed that the twilight flats appear uniquely useful at getting at this calibration term for sky subtraction, and quantifying its drift in time may motivate a more frequent acquisition time. Here I show a figure that quantifies both the wavelength and time dependence of this term. Here are 25 twilight flats from January to August of this year, roughly a few times per month. We see a smooth wavelength dependence across orders, as well as intra order structure that is repeatable from night-to-night. We see a much smaller time dependence with January acquisitions a fraction of a percent higher than, say, August acquisitions. The point-to-point scatter is much larger than the time drift, but it is indeed systematic:

The format of this Pull Request will be to upload a wavelength dependent template---packaged within muler source code---that serves as a lookup table to scale the sky spectrum. The adopted approach should work whether you install the code via pip or the developer version.

[gully]

Ok, I made a smoothed version of the beta value to use as the estimator for beta:

[gully]

A point about uncertainty. The median residual around this beta trend is 0.5%. In some places near the water bands, the residuals are much larger, maybe 30% or greater, so essentially unconstrained. I included the standard deviation of the residuals of the 25 twilight flats as the "uncertainty" around the estimator, so if you want to inspect the uncertainty you may examine the third column of the csv file HPF_sci_to_sky_ratio_beta.csv. In practice there is additional systematic bias in the bias-variance tradeoff for the GP estimator I used.