nanograv / enterprise

ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) is a pulsar timing analysis code, aimed at noise analysis, gravitational-wave searches, and timing model analysis.
https://enterprise.readthedocs.io
MIT License
65 stars 67 forks source link

Adding PintPulsar SSB Vectors #244

Closed Hazboun6 closed 3 years ago

Hazboun6 commented 3 years ago

This PR tries to fix Issue #238 and pulls out the needed vectors from the Pint timing objects in order to do solar system ephemeris modeling. The unit tests pass, but the vectors seem to be off on the order of milliseconds. Before merging I'd like to hear from the Pint team about if this sort of difference is expected. Below are two images that show the vectors and the difference (libstempo-Pint) for Earth and Jupiter. These are for PSR J1911+1347 with the NANOGrav 12.5 year data set.

image

image

scottransom commented 3 years ago

Hmmmm... That seems further off than I would expect. And am I reading things correctly in that the errors only seem to be in Y and Z and not X? This will likely take some investigation.

Hazboun6 commented 3 years ago

@scottransom yeah, that's right regarding the coordinate differences. Below is just the x direction difference, which is more like ~50 nanoseconds. image

vallis commented 3 years ago

Indeed, it would be good to understand this, although it won't affect BayesEphem, since the orbits of the planets are multiplied by M_planet/M_sun to yield mass-perturbation partials, and a few milliseconds out of thousands seconds will be a small correction.

Where does PINT get planetary ephemerides? Are they coordinate-transformed before being reported?

In other news, Jeff, could you merge master into the PR so you get the effects of the new continuous integration?

scottransom commented 3 years ago

@vallis So the high-level PINT routine that is called to get the planet positions is compute_posvels() in toa.py. That calls objPosVel_wrt_SSB() which is in solar_system_ephemerides.py. However, the code that actually gets the positions and velocities from the JPL ephemerides is from the astropy coordinates get_body_barycentric() method. So we probably need to dig in there (and perhaps compare with pyephem or some other codes).

Hazboun6 commented 3 years ago

@vallis I merged the recent CI changes. So you think it's okay to merge this PR and start testing BayesEphem runs with Pint?

codecov[bot] commented 3 years ago

Codecov Report

Merging #244 (93ce894) into master (45c523f) will increase coverage by 0.12%. The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #244      +/-   ##
==========================================
+ Coverage   85.56%   85.69%   +0.12%     
==========================================
  Files          12       12              
  Lines        2723     2740      +17     
==========================================
+ Hits         2330     2348      +18     
+ Misses        393      392       -1     
Impacted Files Coverage Δ
enterprise/pulsar.py 90.98% <100.00%> (+0.74%) :arrow_up:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 45c523f...2584d9a. Read the comment docs.

Hazboun6 commented 3 years ago

@paulray, @luojing1211 These tests are failing on py3.8 for some weird error with Pint's urllib download of ephemeris bsp files. Have either of you seen these with Pint before?

paulray commented 3 years ago

Yeah, we get things like that sometimes. There is a good chance that re-running that test will fix it. There is an East Coast-wide internet issue today so that could be related.

aarchiba commented 3 years ago

The code appears to be correct, basically inverting PINT's subtraction to get back the Astropy output:

https://github.com/nanograv/enterprise/pull/244/files#diff-27db36a2fdf65101108346718f2a100d56a8955c4cb7501dc9f4cafa0cf18384R366

versus

https://github.com/nanograv/PINT/blob/daea8fc7a5bc8b8de974907abb4d87d9bd47c613/src/pint/toa.py#L2032

We can (should) have a round-trip test with PINT to confirm that our generated vectors do agree with Astropy, but if you are depending on PINT already you could add a check of PINT and TEMPO2 SSB-planet vectors against Astropy to the test suite here.

Hazboun6 commented 3 years ago

Thanks @paulray rerunning did the trick!

Hazboun6 commented 3 years ago

@aarchiba thank you for the code cross check. That is very helpful.