refactor of barycentric correction code from #185

[sean-lockwood]

Based on Josh Lothringer's #185 PR. Validation products are available there.

[stscirij]

It's a pity there isn't a working version of the old odelaytime code that would enable a comparison of results. Is there any way of doing an independent test?

[sean-lockwood]

@jlothringer pointed out that the Earth ephemeris used in the IRAF code was very old. He might have more to say on the effectiveness of that comparison...

[jlothringer]

@jlothringer pointed out that the Earth ephemeris used in the IRAF code was very old. He might have more to say on the effectiveness of that comparison...

Yeah, we spent a lot of time (~5 years?) trying to get IRAF results to make sense with a Python port, but ultimately chose to validate the methodology of the new tool using a number of other independent 3rd party barycentric correction tools, namely pintbary, barycorr, and barycorrpy. A lot of those tests are shown in an Odelay_Verification.ipynb notebook linked here: #185 (comment)

We also further validated the calculation of Earth's position via Astropy+JPL Horizons by directly using the SPICE kernels.

[sean-lockwood]

@jlothringer -
Give this version a look and let me know if you see any issues.

The docs are also available for proofreading: https://stistools--189.org.readthedocs.build/en/189/

[sean-lockwood]

@jlothringer -

What happens if a dataset crosses the boundary between two hst_orbfiles?
E.g., od9m72040 over P350000R + P370000R

Is the user expected to concatenate the files to provide coverage? We might want to document this procedure if so.

Update:
One possible answer might be that at least some of the orb files overlap one another:

p350000r.fit  57817 - 57820
p370000r.fit  57819 - 57822

[sean-lockwood]

@jlothringer -

I've also noticed that COS x1dsum files contain the keywords EXPSTRTJ & EXPENDJ. Should we be updating these too?

Additionally, would you expect the various human-readable date/times (DATE-OBS + TIME-OBS) to be updated?

[jlothringer]

@jlothringer -

What happens if a dataset crosses the boundary between two hst_orbfiles? E.g., od9m72040 over P350000R + P370000R

Is the user expected to concatenate the files to provide coverage? We might want to document this procedure if so.

This has always been in the back of my mind - I don't actually know. I figure concatenating them would be the way to go. I didn't know some had overlap - at least some do not with as P2L0000R' and P2O0000R'.

[jlothringer]

@jlothringer -

I've also noticed that COS x1dsum files contain the keywords EXPSTRTJ & EXPENDJ. Should we be updating these too?

Additionally, would you expect the various human-readable date/times (DATE-OBS + TIME-OBS) to be updated?

The COS functionality was something kind of last-minute and I'm unfamiliar with their files. It sounds like it would make sense that those keywords (and any keywords the pipeline or user would use to build a timestamp) should be updated too.

The human readable times, I'm not so sure. I mostly took the file handling functionality from the odelaytime Python port, so it doesn't seem those times were updated in odelaytime? In the comments to the TDATEOBS and TTIMEOBS keywords in the primary header, it does mention that they are UT times, so we would at least have to make sure we didn't change them to TDB. It might make sense to keep these as times at Earth too, but we maybe could be explicit that they are not modified by the script.

[sean-lockwood]

To truncate an orb file for testing purposes:

from astropy.io import fits

# Make a shorter test table:
hdu = fits.HDUList()
with fits.open('p350000r.fit') as f:
    h0 = f[0].header
    h1 = f[1].header
    h0['DEFINEND'] = '2017.066:00:00:00.000'
    h1['LASTMJD'] = 57819.
    d = f[1].data
    hdu.append(fits.PrimaryHDU(header=h0))
    hdu.append(fits.BinTableHDU(data=d[d['Time'] <= 172800], header=h1))
    hdu.writeto('p350000r_truncated.fit', checksum=False, overwrite=True)

To combine two neighboring tables together:

from astropy.table import Table, vstack
from astropy.io import fits

# Stack arrays together, excluding the repeated point in the middle:
t1 = Table.read('p350000r_truncated.fit', hdu=1)
t2 = Table.read('p370000r.fit', hdu=1)
t2['Time'] += (t2.meta['FIRSTMJD'] - t1.meta['FIRSTMJD']) * 24 * 3600.
t = vstack([t1[:-1], t2])
assert len(t) == len(set(t['Time'])), 'Non-unique times found!'

# Combine into a single output FITS file:
hdu = fits.HDUList()
with fits.open('p350000r_truncated.fit') as f1, fits.open('p370000r.fit') as f2:
    h0 = f1[0].header
    h1 = f1[1].header
    h0['DEFINEND'] = f2[0].header['DEFINEND']
    h1['LASTMJD'] = f2[1].header['LASTMJD']
    hdu.append(fits.PrimaryHDU(header=h0))
    hdu.append(fits.BinTableHDU(data=t, header=h1))
    hdu.writeto('combined_orb.fits', overwrite=True)

[sean-lockwood]

@jlothringer -
I added a new function to combine orbfiles, stistools.barycentric_correction.combine_hst_orb_files().

Let me know if you have any feedback on it.

[sean-lockwood]

Ok, I've updated the code to apply to the 'EXPSTRTJ' and 'EXPENDJ' keywords in COS files.

Code should be ready now.

[jlothringer]

I've tested the code against the previous version (pre-re-factor) and all time corrections appear to agree. In those tests, I calculated the corrections for different file types (tag, raw, flt, x1d) and using JPL Horizons versus the HST orbital files for HST's position determination. All agreed as expected with disagreements below 1 ms.

I ran into no errors using COS files. And the orbital file combination function also worked as expected (even when switching the order of the files).

[sean-lockwood]

Thanks, @jlothringer for the thorough check!

@stscirij -
I think this PR is now just pending on your review.

[stscirij]

I pulled the branch and ran the tests and they all passed. Everything looks good to my level of expertise. I did make 1 suggestion, but wouldn't insist on holding everything up if you're OK with keeping the code as it is.
Good job @sean-lockwood and @jlothringer !