Glitches not accounted for when using change_pepoch

[aplangrangian]

Hi, it would seem change_pepoch doesn't account for glitches when forward folding a parameter file. I've tried this twice once with the parameter file below and once by removing all the glitch rows. Changing the epoch results in the same f0, f1 and f2.

Is there a better way to do this? The code I used in python is:

model = models.get_model('/Users/alexlange/Downloads/J0835-4510_tnfit.par') t0 = Time(55197, scale='tdb', format='mjd') model.change_pepoch(t0)

.par file: PSR J0835-4510 EPHEM DE421 UNITS TDB START 54682.71577 FINISH 59500.98129 DILATEFREQ N DMDATA 0.0 NTOA 0.0 CHI2 0.0 RAJ 8:35:20.61149000 0 0.00002000000000000000 DECJ -45:10:34.87510000 0 0.00029999999999999997 PMRA -49.68 0 0.06 PMDEC 29.9 0 0.1 PX 3.5 0 0.2 POSEPOCH 51544.0000000000000000 F0 11.188184713595025169 1 5.0793697580676314082e-08 F1 -1.543185292003074411e-11 1 5.0732618872207964686e-16 F2 1.0455059694879782365e-21 1 9.000466591345223064e-24 PEPOCH 56401.0000000000000000 PLANET_SHAPIRO N TZRMJD 56623.1552660391321759 TZRSITE coe TZRFRQ inf GLEP_1 55408.8000000000000000 GLPH_1 0.004560453727298499 0 0.002775057503383006 GLF0_1 2.1250032625524398e-05 1 2.3106879006399737e-09 GLF1_1 -1.0673117363666227e-13 0 3.4993208846389806e-15 GLF2_1 0.0 GLF0D_1 9.946756155641468e-08 0 8.618210469066635e-09 GLTD_1 10.208572799353567 0 1.7992547326179056 GLEP_2 56555.8082770568799998 GLPH_2 0.8323457251981262 0 0.0032555203771753257 GLF0_2 3.416449630087724e-05 1 2.5190393677039575e-09 GLF1_2 -1.316674947620517e-13 0 3.4676307352359607e-15 GLF2_2 0.0 GLF0D_2 2.5977395910060916e-07 0 9.79572328440877e-09 GLTD_2 9.684181177188016 0 0.6974419246430799 GLEP_3 57734.4848859665900000 GLPH_3 0.003556857729488026 0 0.005022014900437134 GLF0_3 1.5978847458971046e-05 1 2.487134907036957e-09 GLF1_3 -9.088070453917932e-14 0 3.548762554274584e-15 GLF2_3 0.0 GLF0D_3 8.028704134042826e-08 0 1.2122536834256808e-08 GLTD_3 10.288056965477796 0 2.5201798114287834 GLEP_4 58515.5929000000000001 GLPH_4 -0.0006597314072849418 0 0.0025001899198323417 GLF0_4 2.7695690930914995e-05 1 2.2497515502424004e-09 GLF1_4 -9.690934098983648e-14 0 3.4345156804111666e-15 GLF2_4 0.0 GLF0D_4 1.9814862454523047e-07 0 8.468689335602995e-09 GLTD_4 9.657612713042457 0 0.8524673693065955 GLEP_5 59417.6355753565874700 1 0.0013120622178812648036 GLPH_5 0.03 GLF0_5 1.384761391498728e-05 1 2.2361495876293876e-09 GLF1_5 -7.77574990130975e-14 1 2.5608100508571276e-15 GLF2_5 0.0 GLF0D_5 0.0 GLTD_5 0.0 WAVEEPOCH 56401.000000000000000 WAVE_OM 0.001303839121663 0 WAVE1 -1.694743461953648 6.257673555922254 WAVE2 -0.686076945409859 -1.230757244179528 WAVE3 -0.011303342325159 0.168680179535391 WAVE4 -0.034279163921702 0.127362450617002 WAVE5 0.031055694053645 0.045008962565907 WAVE6 0.007814991561899 0.001770432849948 WAVE7 -0.000143057790972 0.005803364468349 WAVE8 0.019195086151892 -0.005821815453603 WAVE9 0.002697263424840 0.008846660164680 WAVE10 0.000500044127211 -0.002373484356075 WAVE11 0.000315522922854 0.000262656616375 WAVE12 0.000695465263366 0.000001269177667 WAVE13 -0.000937497204630 -0.002217897657384 WAVE14 0.000217446636757 -0.000064941464240 WAVE15 -0.000475148850835 -0.000489775035376 WAVE16 0.000081348732408 -0.000834377493197 WAVE17 -0.000434497596120 -0.000298372190613 WAVE18 -0.000416644330520 0.000007474225227 WAVE19 -0.000181727943210 -0.000840901781650 WAVE20 -0.000184700319328 0.000356894613688 WAVE21 -0.000497937041307 0.000206231673858 WAVE22 -0.000184189806883 -0.000349510625280 WAVE23 0.000038889055190 0.000246789551378 WAVE24 -0.000156869770325 0.000077908312360 WAVE25 0.000135550579124 -0.000012612358625 WAVE26 0.000106675399771 0.000123529282148 WAVE27 -0.000199914663131 -0.000044596956841 WAVE28 0.000210962161736 -0.000140626823032 WAVE29 0.000055517112340 0.000173895130347 WAVE30 -0.000096698721675 -0.000035741701158 WAVE31 0.000195625863721 0.000009808046371 WAVE32 0.000020642657974 0.000099171880329 WAVE33 0.000060953078391 -0.000035869650167 WAVE34 0.000110035287220 0.000052512977022 WAVE35 -0.000092684741976 0.000043395711838 WAVE36 0.000066871626283 -0.000071332311983 WAVE37 0.000037358350150 0.000019927924967 WAVE38 -0.000071378728562 0.000001531016208 WAVE39 0.000068128744330 -0.000072118474899 WAVE40 -0.000014276934318 0.000073646177479 WAVE41 -0.000058428736416 -0.000041945624187 WAVE42 0.000104636365959 -0.000013302310891 WAVE43 0.000002519859013 0.000078380660551 WAVE44 -0.000028443949608 -0.000061361393561 WAVE45 0.000111409016616 0.000034675769071 WAVE46 -0.000062527741688 -0.000017381531695 WAVE47 0.000014785310192 -0.000059021924641 WAVE48 0.000065472314363 0.000035374470574 WAVE49 -0.000068541641796 0.000044612862939 WAVE50 0.000039464729172 -0.000014566512658 WAVE51 -0.000000959096499 0.000034757214926 WAVE52 -0.000034216329690 -0.000037029590200 WAVE53 0.000052253537958 -0.000004715266746 WAVE54 -0.000031413664185 0.000019831142146 WAVE55 -0.000009917216846 -0.000024445278488 WAVE56 0.000027073758610 0.000003937647395 WAVE57 -0.000022747392223 0.000009379468023 WAVE58 0.000002455924579 -0.000030466827594 WAVE59 -0.000006202266056 0.000007159236186 WAVE60 -0.000006659021331 0.000018384200242 WAVE61 0.000031558844106 -0.000031284045690 WAVE62 -0.000013782605652 0.000035052265969 WAVE63 0.000002522518717 -0.000000893530619 WAVE64 0.000022833122635 -0.000016329838228 WAVE65 -0.000004614052348 0.000014046179669 WAVE66 -0.000021764354187 -0.000033272898046 WAVE67 0.000003887258811 0.000003628469504 WAVE68 -0.000005350939582 0.000002702094386 WAVE69 0.000006379738064 -0.000022028161116 WAVE70 0.000001382615112 0.000003608553376 WAVE71 -0.000025027928354 0.000000893570457 WAVE72 0.000000227410853 -0.000006404247168 WAVE73 -0.000006668303772 0.000015993917048 WAVE74 0.000003301263375 0.000016503362263 WAVE75 0.000024569571257 -0.000008939529730 WAVE76 -0.000004646437074 0.000022297861715 WAVE77 0.000016973976288 -0.000001526182863 WAVE78 0.000000236676584 0.000001547969098 WAVE79 0.000011434045266 0.000003725920380 WAVE80 0.000013505716856 -0.000020316852110 WAVE81 0.000004866715629 0.000014982484131 WAVE82 0.000011940948987 0.000011583407268 WAVE83 -0.000010031327197 -0.000005108594156 WAVE84 -0.000004724021013 -0.000001019847023 WAVE85 -0.000003820422785 -0.000021109081774 WAVE86 0.000002329194444 -0.000006542073669 WAVE87 -0.000000878288349 -0.000000426430253 WAVE88 0.000001642606918 -0.000008527251279 WAVE89 -0.000001101821063 -0.000005148784844 WAVE90 0.000004772438380 0.000004412842670

LATPHASE 0.9356150486 LATAMP -4.68215 0.06 LATIDX 6.31508 0.12 LATNHARM 90 LATMINWT 0.990000 LATPHOT 483545 LATFN 2 LATLOGL 404552.691 LATHTEST 941287.747 LATWN 2.781e-10

[aplangrangian]

Just adding a comment to give results from both attempts: Par file including glitches: Spindown( floatParameter( F0 11.184356169900447467 (Hz) +/- 5.0793697580676314e-08 Hz frozen=False), MJDParameter( PEPOCH 59297.0000000000000000 (d) frozen=True), floatParameter( F1 -1.517025227117889133e-11 (Hz / s) +/- 5.073261887220796e-16 Hz / s frozen=False), floatParameter( F2 1.0455059694879782365e-21 (Hz / s2) +/- 9.000466591345223e-24 Hz / s2 frozen=False)), Par file not including glitches: Spindown( floatParameter( F0 11.184356169900447467 (Hz) +/- 5.0793697580676314e-08 Hz frozen=False), MJDParameter( PEPOCH 59297.0000000000000000 (d) frozen=True), floatParameter( F1 -1.5170252271178891328e-11 (Hz / s) +/- 5.073261887220796e-16 Hz / s frozen=False), floatParameter( F2 1.0455059694879782365e-21 (Hz / s2) +/- 9.000466591345223e-24 Hz / s2 frozen=False)),

[dlakaplan]

Hi, @aplangrangian I'm not sure what to expect here. Should the changes in the F0/F1/F2 values depend on the glitch? I wouldn't think so, but maybe I'm missing something. Changing the PEPOCH does not change the glitch parameters but I'm not sure if it should.

[paulray]

I guess the proof is in the pudding. If you have a set of TOAs that span the glitch, does model.phase() return the same phases before and after the change_pepoch(). If so, then there is nothing wrong.

[dlakaplan]

I think it works OK, at least with a small test using only a single glitch.

from astropy import units as u, constants as c
from astropy.time import Time
import numpy as np
import pint.logging
import pint.models
import pint.simulation
import io
import copy

pint.logging.setup("WARNING")
par = """
PSR J0835-4510
EPHEM DE421
UNITS TDB
START 54682.71577
FINISH 59500.98129
DILATEFREQ N
DMDATA 0.0
NTOA 0.0
CHI2 0.0
RAJ 8:35:20.61149000 0 0.00002000000000000000
DECJ -45:10:34.87510000 0 0.00029999999999999997
PMRA -49.68 0 0.06
PMDEC 29.9 0 0.1
PX 3.5 0 0.2
POSEPOCH 51544.0000000000000000
F0 11.188184713595025169 1 5.0793697580676314082e-08
F1 -1.543185292003074411e-11 1 5.0732618872207964686e-16
F2 1.0455059694879782365e-21 1 9.000466591345223064e-24
PEPOCH 56401.0000000000000000
PLANET_SHAPIRO N
TZRMJD 56623.1552660391321759
TZRSITE coe
TZRFRQ inf
GLEP_1 55408.8000000000000000
GLPH_1 0.004560453727298499 0 0.002775057503383006
GLF0_1 2.1250032625524398e-05 1 2.3106879006399737e-09
GLF1_1 -1.0673117363666227e-13 0 3.4993208846389806e-15
GLF2_1 0.0"""

m = pint.models.get_model(io.StringIO(par))
t = pint.simulation.make_fake_toas_uniform(54000, 58000, 300, m, add_noise=True)
ph = m.phase(t)

m2 = copy.deepcopy(m)
m2.change_pepoch(Time(57500, format="mjd", scale="tdb"))
ph2 = m2.phase(t)
assert np.allclose(ph, ph2, atol=1e-10)

[aplangrangian]

Changing the PEPOCH does not change the glitch parameters but I'm not sure if it should.

I wouldn't expect it to either. However when trying to determine the F0/F1/F2 spin down values of the pulsar, the change_pepoch() doesn't seem to account for the glitches. Wouldn't the F0/F1/F2 at any given time be dependent on the preceding glitches (ie the F0 would increase).

For example: at 56401 MJD F0/F1 are: 11.188184713595025169 and -1.543185292003074411e-11. Fast forward 5000 days just based on these parameters (no glitches) and we expect: 11.184356169900447467 and -1.5170252271178891328e-11. However, if I do include those glitches and their parameters shown above, I get the same result when we would expect the glitches to spin-up the pulsar and deviate the F0/F1 from the fit.

Or am I misunderstanding something fundamental either with PINT or with pulsar nature?

[aplangrangian]

Did not mean to close. Sorry.

Also thank you both for your quick responses!

[dlakaplan]

I think @paulray 's answer was the right way to think about it. Changing the PEPOCH is used to minimize covariance when fitting, but the phase as a function of time should still be the same. Glitches just represent a delta function (in one or more parameters). The background polynomial represented by PEPOCH/F0/F1/F2/... doesn't change, only its parameterization. The glitch amplitude will also not change, so the same phase(t) is reproduced. At least that's what I understand.

[aplangrangian]

Okay that makes sense. I misunderstood the functionality then.

I suppose if I wanted the spin down parameters at a specific time, I'd have to either fit TOAs of just that time range (resulting in larger errors) or do the calculation manually with all the glitch terms?

[paulray]

I think there is a little misunderstanding here. When a model includes glitches, the glitches are required to model the spin evolution of the pulsar. You can't just change the PEPOCH to be after the glitch and expect to get the right spin period from just reading F0/F1/F2... If you want the local period at some time in a model that has glitches, you can compute that using the d_phase/d_toa. Do you want the barycentric period at that time, or topocentric at some observatory? The latter is traditionally done using polycos, and PINT supports that. In fact, you can generate barycentric polycos as well. Can you give a bit more detail on your use case?

[paulray]

I should have added that the code you probably want is model.d_phase_d_toa(ts), where model is your TimingModel object and ts is a TOAs object. You can look at the test_d_phase_d_toa.py code for an example of using this, where it compares the result to the polycos version.

[aplangrangian]

Can you give a bit more detail on your use case?

In my case, I am looking to find the F0/F1 at a different epoch than given in a .par file computed from Fermi LAT. This will then be compared to a different observation from NuSTAR - so barycentric.

Looking into the test_d_phase_d_toa.py, I would imagine what I'm looking for is the "eval_spin_freq()" function then.

[dlakaplan]

That function uses polycos, which might be OK, but the other function d_phase_d_toa is a bit more direct and will still get you F(t). And you can then compute numerical derivatives to get Fdot(t) etc.

[aplangrangian]

Okay! I'll look into that then! Thank you!

[paulray]

Yes, d_phase_d_toa() is what you want. However, just be careful. If you NuSTAR observation is not between the START and FINISH in the par file, you won't get the right answer. The F0/F1/F2 plus glitches don't fully describe the spin evolution of Vela, because of significant timing noise. This par file parameterizes that using WAVE parameters, which are a Fourier representation of the excess residuals. The d_phase_d_toa() function will correctly evaluate the phase evolution including the WAVE parameters and give you the correct instantaneous frequency at any time between START and FINISH. Extrapolating beyond FINISH will not be correct, as the timing noise is a stochastic process. Also, if you are using an FTOOL like barycorr to barycenter your NuSTAR photons, you may not be correctly accounting for the proper motion. This is likely to be a small effect for you, but just to let you know that you can compute phases for your NuSTAR photons directly using PINT and then the full astrometric model will be applied.

[aplangrangian]

Okay it looks like d_phase_d_toa() worked more or less perfectly. I did find some error loading in the Fermi TOAs:

Traceback (most recent call last):

Cell In[82], line 1 model.d_phase_d_toa(fermi_data)

File /opt/anaconda3/envs/fermipy/lib/python3.9/site-packages/pint/models/timing_model.py:1577 in d_phase_d_toa dt_array = [dt.value] copy_toas.ntoas dt._unit

AttributeError: 'list' object has no attribute 'ntoas'

Not sure if this pops up for anyone else or if it's due to a difference between the load_Fermi_TOAs object and TOAs input for d_phase_d_toa(), but I found a workaround(?) using pint.fermi_toas.toa.get_TOAs_list()

I've attached my small rudimentary code below in case it helps someone else.

from pint import fermi_toas

fermi_data = pint.fermi_toas.load_Fermi_TOAs('/Path/To/Fermi/Ft1.fits') index = whatever_mjd fermi_data[-index].mjd

model = models.get_model('/Path/To/Par/File.par')

test = pint.fermi_toas.toa.get_TOAs_list(fermi_data)

new_f0 = model.d_phase_d_toa(test) print(new_f0[-index])

[dlakaplan]

Yes, this is a known issue. load_Fermi_TOAs is somewhat old and returns a list of TOA objects, not a TOAs object. As you found you need to use get_TOAs_list on the output of load_Fermi_TOAs to turn it into a TOAs object. Alternately, you can use get_Fermi_TOAs and that should work.

We realize this is confusing and are in the process of simplifying the interface (#1541 #1543 )

[aplangrangian]

Awesome. Just thought I should point it out. Otherwise this was very helpful and solves my problems. I suppose I should close?