Discrepancy in phase prediction between PINT and tempo2 polycos.

[theXYZT]

Hi,

I am trying to figure out what's causing a particular discrepancy in pulsar phase prediction between PINT and generating polycos with tempo2.

Here is a par file I am using, b1937.par (it is a modified version of the NanoGrav 12.5 year narrowband par file for B1937+21):

PSR              B1937+21
LAMBDA   301.9732445257902  1     0.0000000054618
BETA      42.2967522555839  1     0.0000000058937
PMLAMBDA           -0.0184  1              0.0041
PMBETA             -0.4007  1              0.0060
PX                  0.2819  1              0.0509
POSEPOCH        55599.0000
F0    641.9282322943167856  1  0.0000000000229207
F1     -4.330889418718D-14  1  1.374778598040D-19
PEPOCH        55599.000000
DM                 71.0201
DMEPOCH          58245.375
SOLARN0               0.00
EPHEM               DE436
ECL                 IERS2010
CLK                 TT(BIPM2017)
UNITS               TDB
TIMEEPH             FB90
T2CMETHOD           TEMPO
CORRECT_TROPOSPHERE N
PLANET_SHAPIRO      N
DILATEFREQ          N
NTOA                 17024
TRES                  1.82
TZRMJD  55616.41780189249631
TZRFRQ            1742.942
TZRSITE                  GB
MODE                     1

When I use tempo2 to generate polycos, it generates the following relevant entry:

B1937+21    7-May-18   90000.00   58245.37500000000            71.020167 -0.706 -6.330
 146774936445.058162  641.928232294317   ao   90   12   327.000                
 -3.00147350705225993e-07  2.72493109523728450e+00 -1.19447953427931179e-04
 -5.36771151182867317e-09  1.79612760415288616e-10  4.84839643606338133e-13
  2.66425801966225905e-15 -5.05089187190312671e-16 -2.43473904299295782e-19
  2.19028473393043386e-19 -9.16295733411736226e-23 -3.35126294498144066e-23

Using PINT, I execute the following code:

import pint.models
import pint.toa

m = pint.models.get_model(folder / "b1937.par")

ts = pint.toa.get_TOAs_list([pint.toa.TOA((58245, 0.375), 0, freq=327 * u.MHz, obs="ao")],
                            ephem=m.EPHEM.value, bipm_version="BIPM2017")

ph = m.phase(ts, abs_phase=False)

Here are the results:

Phase prediction at MJD 58245.375:
tempo2 polyco:  146774936445.0581616998526493
pint-pulsar:    146774936445.0755193680524826

The fractional part of the predicted phase is different between them. What am I doing wrong that's causing this discrepancy?

[paulray]

For absolute phase comparisons, always use abs_phase=True. How does it look then?

[dlakaplan]

I'm trying to investigate using:

import io
import numpy as np
from astropy import units as u
from pint.polycos import Polycos
from pint.models import get_model
import pint.toa as toa
import pint.simulation
import pint.logging

pint.logging.setup(level="WARNING")
freq = 327
obs = "ao"

mjd = 58245.375
mjd_start = mjd
mjd_end = mjd + 1
mjds = np.linspace(mjd_start, mjd_end, 11)
model = get_model("b1937.par")

ts = toa.get_TOAs_list(
    [toa.TOA(mjd, obs=obs, freq=freq) for mjd in mjds], ephem=model.EPHEM.value
)
polyco_text = """B1937+21    7-May-18   90000.00   58245.37500000000            71.020167 -0.706 -6.330
 146774936445.058162  641.928232294317   ao   90   12   327.000                
 -3.00147350705225993e-07  2.72493109523728450e+00 -1.19447953427931179e-04
 -5.36771151182867317e-09  1.79612760415288616e-10  4.84839643606338133e-13
  2.66425801966225905e-15 -5.05089187190312671e-16 -2.43473904299295782e-19
  2.19028473393043386e-19 -9.16295733411736226e-23 -3.35126294498144066e-23"""

q = Polycos()
q.read_polyco_file(io.StringIO(polyco_text))

p = Polycos()
p.generate_polycos(model, mjd_start, mjd_end, obs, 90, 12, freq)

ph_pintpolyco = p.eval_abs_phase(mjds[0])
ph_tempo2polyco = q.eval_abs_phase(mjds[0])
ph_model = model.phase(ts, abs_phase=True)
ph_model_noabs = model.phase(ts, abs_phase=False)

print(
    f"Tempo2 polyco prediction from PINT: {(ph_tempo2polyco.int + ph_tempo2polyco.frac)[0]}"
)
print(
    f"PINT polyco prediction from PINT: {(ph_pintpolyco.int + ph_pintpolyco.frac)[0]}"
)
print(f"PINT model prediction from PINT: {(ph_model.int + ph_model.frac)[0]}")
print(
    f"PINT model prediction from PINT (abs_phase=False) {(ph_model_noabs.int + ph_model_noabs.frac)[0]}"
)

This mimics what is in tests/test_polycos.py, which compares the phase from PINT's model and the polyco. But then it also reads in the polyco above and compares PINT's phase with abs_phase=True and `False.

[dlakaplan]

When I run it with the starting parameters (Arecibo, 327 MHz) I get:

Tempo2 polyco prediction from PINT: 146774936445.05817
PINT polyco prediction from PINT: 145809053984.0521
PINT model prediction from PINT: 145809053984.0521
PINT model prediction from PINT (abs_phase=False) 146774936445.07553

So PINT's model and PINT's polyco agree (which the test tests for). That is good.

Both of those disagree even at the integer level with the tempo2 polyco.

The tempo2 polyco and the result with abs_phase=False disagree at the ~0.02 cycle level, as @theXYZT found.

[dlakaplan]

I tried this using fake TOAs returned by pint.simulation.make_fake_toas. Those actually disagreed even between the PINT model and PINT polyco. The procedure above (toa.get_TOAs_list) is what the test does. This may be because it doesn't apply observatory clock corrections (I think).

[dlakaplan]

The polyco output by PINT is:

B1937+21   07-May-18   90000.00   58245.37500000000            71.020100  0.000  0.000
 145809053984.052083  641.928232294317   ao   90   12   327.000                
 -1.66452517437094157e-07  2.72492679128434689e+00 -1.19429226408798771e-04
 -5.18832680537685495e-09  1.83698955572492937e-10  1.37495035056639063e-13
 -1.27410974958984366e-16 -1.35407839776891287e-16  1.14459885507404128e-19
  6.17368677135918858e-20 -3.46339609727139525e-23 -1.04182120118969413e-23

I notice that it doesn't have the Doppler shift term. I don't know if this is intentional, or if this could be the cause of the difference. It looks like the PINT polyco code doesn't compute anything with the Doppler term. If I compute the difference in DM delay with and without that term (and it's effect on the topocentric RF) I get a difference of 0.4 ms, so a significant fraction of a cycle. Probably too big or too to actually account for the differences above, though.

[theXYZT]

As @dlakaplan found, I used abs_phase = False because setting it to True resulted in completely different integer parts for the predicted phase. When I contributed to the polyco part of PINT 2 years ago, I was using abs_phase=False to generate polycos. This was then changed in #1115, it seems.

@dlakaplan: The dopper shift term is not really used for phase prediction, so the polyco writer just ignores that term for simplicity.

[dlakaplan]

@theXYZT : could you try polyco generation for a few other observatories and/or frequencies? I'd like to understand if this is an issue related to clock corrections or other topocentric issues.

[paulray]

For specificity, could you give the command used to generate the tempo2 polyco file?

[paulray]

I think polycos should always be made with abs_phase=True since you want the phase to be 0.0 at the fiducial TOA defined by the TZR* parameters.

[theXYZT]

For specificity, could you give the command used to generate the tempo2 polyco file?

tempo2 -tempo1 -f b1937.par -polyco "58245.34 58245.42 90 12 12 ao 327"

[dlakaplan]

One thing to note: if you generate TOAs using the command above:

ts = pint.toa.get_TOAs_list([pint.toa.TOA((58245, 0.375), 0, freq=327 * u.MHz, obs="ao")],
                            ephem=m.EPHEM.value, bipm_version="BIPM2017")

It does do the BIPM and GPS corrections. So the resulting TOA is actually at MJD 58245.37500000032. This difference (3.2e-10d = 2.7438950005900858e-05s) I think comes from the BIPM correction.

But the difference between MJD 58245.37500000032 and 58245.375 is 2.7438950005900858e-05s = 0.01761533 cycles, which is almost (but not exactly) the original difference above. If I instead turn off those corrections and get a MJD of 58245.375 for the comparison, I get a phase (using abs_phase=False) of 146774936445.05777. Not the same as the original, but now differs by -0.00039635 cycles instead of 0.017357668199833302. [Note that PINT still lists a clock correction of ~ -2e-7s, but the MJD now seems to be exactly 58245.375]

And I'm still not sure why this is only close when using abs_phase=False

[theXYZT]

I assumed I was supposed to include BIPM since it is mentioned in the par file and I assumed tempo2 would probably use this information, so PINT should too. I have little to no understanding of how tempo2 / PINT actually do phase prediction so am mostly poking in the dark here.

[paulray]

The BIPM correction is 27 µs, as David noted above. There could be some discrepancy in how the 2 codes handle clock corrections within the polyco computation, or what time scale is assumed when you ask for a polyco to be evaluated. So 27 µs differences wouldn't surprise me too much.

What I don't get is why PINT and Tempo2 agree better when you have abs_phase=False. I had thought that Tempo2 computed polycos giving absolute phases relative to TZRMJD.

[theXYZT]

Is there a much more minimal parfile somewhere that might help narrow down where the discrepancy comes from?

[paulray]

Try making TZRSITE @ and TZRFRQ 0.0 and see if that changes things

[dlakaplan]

For specificity, could you give the command used to generate the tempo2 polyco file?

tempo2 -tempo1 -f b1937.par -polyco "58245.34 58245.42 90 12 12 ao 327"

I don't do this regularly, but I tried that and I got different results. It was in a different format: is there some plugin I need? And if I try to get the right block I find:

TEMPO2: POLYCO TEMPO1 emulation
B1937+21  
 7-May-18
90000.00   
58245.37500000000000000000
71.0201674144            
-0.7062657
1.389  
146774936445.058227583765984
641.92823229431678561108
ao   
90   
12   
327.000   
0.000000 0.00000
-3.357534738702538500837446733271e-07
2.724931187635792411002796931108e+00
...

so the absolute phase changed once again to 146774936445.058227583765984.

[theXYZT]

@dlakaplan It should generate two files. You want the one with extension .dat

[dlakaplan]

got it (there were 2, although both .dat)

[theXYZT]

got it (there were 2, although both .dat)

My bad! I was thinking of the .tim that also gets generated.

[dlakaplan]

If I look at the right file I still get small differences compared to what @theXYZT had posted (assuming this is the right entry):

B1937+21    7-May-18   90000.00   58245.37500000000            71.020167 -0.706  1.389
 146774936445.058228  641.928232294317   ao   90   12   327.000                
 -3.35753473870253850e-07  2.72493118763579245e+00 -1.19448419458989229e-04
 -5.35038518677386951e-09  1.79234446086681942e-10  4.54332746231087971e-13
  3.27638352816407964e-15 -4.82500512582389774e-16 -6.32884393364234230e-19
  2.11944854967265608e-19 -1.13728262846463262e-23 -3.27590308951086241e-23

So the phase is different by ~1e-4, and the first coefficient is also different. Others are too. I don't know if this is due to a different clock file or something else: I don't use tempo2 standalone often, so something could be out of date.

[theXYZT]

For reference, I used tempo2==2022.05.1 via conda-forge.

[dlakaplan]

Still looking at this. I tried to simplify the par file:

• Set F1 to 0
• Set TZRSITE to @
• Set TZRFREQ to 0
• Set DM to 0 This was trying to remove questions about barycentering, DMs, Doppler shifts, etc. I then computed polycos and phases for site=@ (which are frequency-independent, as expected).

Comparing the PINT phases with abs_phase=True vs False, they differ by a constant amount which is equal to the difference between TZRMJD and PEPOCH, with an extra 66s (for some time standard difference? this seems familiar but I don't remember explicitly).

Comparing the PINT abs_phase=False to the polyco phase, I find they differ by about 0.55 cycles (almost but not exactly constant). I don't know why that is. I do see a number of hardcoded constants in the tempo2 polyco code (for date/time conversions): I don't know if this are sufficient.

More important, I don't know why the tempo2 polyco predictions are only close to those with abs_phase=False. It looks like tempo2 is using the TZR data.

[kerrm]

66s is the difference between UTC and TT: 32.184s for TT to TAI, and an additional 34 (I guess) for leap seconds

On Fri, 5 Aug 2022 at 15:18, David Kaplan @.***> wrote:

Still looking at this. I tried to simplify the par file:

• Set F1 to 0
• Set TZRSITE to @
• Set TZRFREQ to 0
• Set DM to 0 This was trying to remove questions about barycentering, DMs, Doppler shifts, etc. I then computed polycos and phases for site=@ (which are frequency-independent, as expected).

Comparing the PINT phases with abs_phase=True vs False, they differ by a constant amount which is equal to the difference between TZRMJD and PEPOCH, with an extra 66s (for some time standard difference? this seems familiar but I don't remember explicitly).

Comparing the PINT abs_phase=False to the polyco phase, I find they differ by about 0.55 cycles (almost but not exactly constant). I don't know why that is. I do see a number of hardcoded constants in the tempo2 polyco code (for date/time conversions): I don't know if this are sufficient.

More important, I don't know why the tempo2 polyco predictions are only close to those with abs_phase=False. It looks like tempo2 is using the TZR data.

— Reply to this email directly, view it on GitHub https://github.com/nanograv/PINT/issues/1369#issuecomment-1206783554, or unsubscribe https://github.com/notifications/unsubscribe-auth/AF3LAG72XIFNXSRKVN4MXWDVXVSIBANCNFSM55SXYJPQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[dlakaplan]

OK, yes, thanks. I knew it looked familiar.

[dlakaplan]

A good test to try: if you generate a polyco for time==TZRMJD, what is the phase?