libstempo
libstempo copied to clipboard
Question: Return phase from an input time array and ephemerides file
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!
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
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).
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)?
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.
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?
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")
...
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!