Skip to content

Add Barycentric Correction Tool#185

Closed
jlothringer wants to merge 6 commits into
spacetelescope:masterfrom
jlothringer:master
Closed

Add Barycentric Correction Tool#185
jlothringer wants to merge 6 commits into
spacetelescope:masterfrom
jlothringer:master

Conversation

@jlothringer

Copy link
Copy Markdown
Contributor

We've developed a new script, barycentric_correction.py, to replace the old functionality of the odelaytime IRAF routine. This new tool uses modern tools (astropy, JPL Horizons) for coordinate transformations and positional queries. Users can still use older HST orbital files, which will agree with the JPL Horizons calculations to within about 1 ms.

The script works with any STIS file (tag, raw, flt, x1d) and inherits its file modification logic from the original odelaytime port. The script has been tested against the old Python port of odelaytime and has been found to agree within about 1 ms. Standalone functions are also available to calculate barycentric corrections without file modifications.

A couple tests have also been added to tests/ with the test files available at this box link. A file in docs/ has been added and init.py has also been updated.

@jlothringer

Copy link
Copy Markdown
Contributor Author

I also have this testing script and notebook that may be helpful. I also attach the odelaytime.py port from before, which was tested against.
odelaytime.py
Odelay_Verification.ipynb
final_odelay_testing.py

@jkcarlberg

Copy link
Copy Markdown

I assume the decision to define unit transformations were adopted from the original "odelaytime" python conversion, but it seems that these should really rely on astropy units and constants to avoid limiting the precision unnecessarily. (E.g., LYPERPC = (1 * u.pc).to('lightyear').value)

Of note, LYPERPC (defined as 3.261633 in the code) disagrees with the astropy calculation beginning in the 4th decimal place
In [21]: (1. *u.pc).to('lightyear')
Out[21]: <Quantity 3.26156378 lyr>

@sean-lockwood

Copy link
Copy Markdown
Member

Will need to add it to the sphinx index:
https://github.com/spacetelescope/stistools/blob/master/doc/source/index.rst

Comment thread tests/test_barycentric_correction.py Outdated

def test_raw_JPL(self):
"""Compare output for a x1d file."""
self.get_data("input", "od9m97020_raw.fits".fits")

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like some typos here in the filename. (And below on L26)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch- done

Comment thread tests/test_barycentric_correction.py Outdated

bary_corr("od9m97020_raw.fits", hst_orb="p2o0000r_TIME.fit")

outputs = [(output, "od9m97020_raw_orbfile_out.fits")]

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jlothringer -
Let me know where this file lives and I can add it along with the input data to Artifactory.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Linked from the PR description here

Comment thread stistools/barycentric_correction.py Outdated

#should be done on the tag file, right? cuz changes will propogate to, like, inttag
#nah, make it generic, may not TIME-TAG data!
def bary_corr(table_names, verbose=True, distance=1e9, hst_orb = None,

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For consistency with other functions in stistools, I think we should either rename your main function (and associated references) to barycentric_correction or rename the module (file name) to bary_corr/bary_corr.py.

I'd vote for the former.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with former- there is already another tool called barycorr: https://github.com/tronsgaard/barycorr

Comment thread tests/test_barycentric_correction.py Outdated
def test_raw_JPL(self):
"""Compare output for a x1d file."""
self.get_data("input", "od9m97020_raw.fits".fits")
#output = "od9m97020_barycorr_raw_out.fits"

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The output variable is not defined but used in the comparison. Fix in both tests.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Comment thread tests/test_barycentric_correction.py Outdated

bary_corr("od9m97020_raw.fits")

outputs = [(output, "od9m97020_raw_JPL_out.fits")]

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The second filename should be the reference we upload to Artifactory. You should probably rename out to ref in both tests.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed and renamed.

Comment thread stistools/barycentric_correction.py Outdated
including barycorr, astroutils, and pintbary.

HST's changing location around the Earth can lead to time-delay
differences of up to ~23 seconds. HST's location can be determined

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

23 s or ms?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ms, though that's relative to the geocenter, so true maximum error is more like 42.5 ms: https://www.wolframalpha.com/input?i=%28Earth+diameter%2B2*550+km%29+%2F+speed+of+light. Fixed.

@jlothringer

Copy link
Copy Markdown
Contributor Author

I assume the decision to define unit transformations were adopted from the original "odelaytime" python conversion, but it seems that these should really rely on astropy units and constants to avoid limiting the precision unnecessarily. (E.g., LYPERPC = (1 * u.pc).to('lightyear').value)

Of note, LYPERPC (defined as 3.261633 in the code) disagrees with the astropy calculation beginning in the 4th decimal place In [21]: (1. *u.pc).to('lightyear') Out[21]: <Quantity 3.26156378 lyr>

I've updated the script to use values derived from the Astropy units where appropriate.

sean-lockwood added a commit that referenced this pull request May 28, 2026
Authored-by: jlothringer <jlothringer@stsci.edu>
Contributions from:  Sean Lockwood, Chris Hayes, Robert Jedrzejewski
@sean-lockwood

Copy link
Copy Markdown
Member

Code merged in #189. Closing this PR.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants