Question: Return phase from an input time array and ephemerides file

[TarekHC]

Hi!

I just discovered this library, and I was wondering if it has the functionality generally used within optical/gamma astronomy... Is it currently possible to provide a time array and ephemerides file, and return the phase of each time value according to the ephemerides?

Maybe its a super dumb question... Apologies if that is the case, as I'm definitely not an expert. If it is possible, would it be possible to get a small snippet so I can try to work it out?

Thanks!

[mattpitkin]

I think you can do something similar to what you want using the fakepulsar function within toasim.py (see this notebook), e.g.:

Note: don't trust me that this works correctly - I could well have misinterpreted something!

from libstempo.toasim import fakepulsar
import numpy as np

# we'll use the par file given with libstempo
parfile = "B1953+29_NANOGrav_dfg+12.par"

# create a set of times at which you want to get the phase
obstimes = np.arange(53000,54800,30) + np.random.randn(60)

toaerr = 1e-12  # set some very small observation errors so the code doesn't complain that they are zero

# create fake pulsar data
psr = fakepulsar(parfile=parfile, obstimes=obstimes, toaerr=toaerr, iters=0)  # iters is zero so that it doesn't try to alter the TOAs

# get the residuals (given in seconds)
res = psr.residuals()

# convert residuals to phase (in cycles)
phase = res * psr["F0"].val

# add on the nearest pulse number (if you want to)
fullphase = psr.pulsenumbers(updatebats=False, formresiduals=False) + phase

[mattpitkin]

In fakepulsar you can use the freq argument if you want to set the observation frequency (in MHz) and the observatory argument to set the observatory name (using the tempo2 name codes given in the fourth column of this file).

[TarekHC]

Wow, thank you for answering so quickly! :D

I will test this right away! Although I have another simple question: I understand the times to be used are UTC times (in MJD) taken at the observatory location (so not barycentered or anything)?

[vallis]

Thanks Matt for replying! I can do my bit and say that the times are indeed meant to be the pulse arrival times at the observatory. Time conversions, barycentering, etc. are done internally to compute residuals with respect to the timing model, which includes spin evolution, binary dynamics, etc.

[TarekHC]

Thank you both! @vallis and @mattpitkin !!

Good news and bad news:

• I was able to make it work relatively easily. So thank you!
• I seem to have hit the float64 precision wall of numpy arrays... My times require precision well beyond the ms level (it is a very high speed optical camera, essentially). I'm able to reproduce the Crab pulsar lightcurve with other software (so I know it is not the problem of the data), but I would like to switch to python.

Taking a look to the notebook you linked, I saw I can provide observation times as a ".tim" file. I don't have a lot of experience with that format. Is there any way I can provide data with enough precision in a different way? Or I should create a ".tim" file for sure?

[mattpitkin]

Within fakepulsar it actually constructs a .tim file based on the supplied obstimes values. You can pass it a np.array with float128 precision and it should still work ok, e.g.

import numpy as np

# float128 array (of zeros)
obstimes = np.zeros(3, dtype=np.float128)

# add in some times
obstimes[0] = np.float("56543.969253958453854835495934")
...

[TarekHC]

Hi again @mattpitkin

Ok, I'm actually ashamed I didn't know there is actually a np.float128...

Not sure it is relevant to you guys... But I'm not able to reproduce the Crab lightcurve I was expecting (while I do get it with the other C software I use, same data).

A quick question: is it normal that the phase distribution I get is not as even as I would expect (while my sampling is constant over 5 min of data, very high frequency)?

I'm providing a Crab pulsar .par file:

PSRJ                 J0534+2200
RAJ                  05:34:31.972
DECJ                 22:00:52.07
F0                   29.5986693695511
F1                   -3.67823E-10
F2                   9.01E-21
TRES                 27.028242047
PEPOCH               59260
POSEPOCH             59260
START                59246
FINISH               59274
TZRMJD               59260.000000046
TZRFRQ               0.0
TZRSITE              COE
EPHEM                DE200
EPHVER               2
CLK                  TT(TAI)
UNITS                TDB
PLANET_SHAPIRO       Y
CORRECT_TROPOSPHERE  N
T2CMETHOD            IAU200B
TIMEEPH              FB90

While times are in MJD. I'm just surprised I'm not reconstructing the Crab lightcurve, given that over such a short observation time I would expect that most of the problems would not really affect much... But I may be totally wrong!