Pulsation Analysis

Modules to analyze events with pulsation signal.

Crab Timing

Several tools developed only for Crab’s analysis

tatpulsar.pulse.Crab.retrive_eph.get_par(mjd, eph_df)

get the best JBL Crab ephemeris for a given time

Parameters

mjd : float

the time (MJD) to get the corresponding ephemeris parameters

eph_df : DataFrame

the pandas DataFrame (return from function retrieve_ephemeris())

Returns

par : list

The best parameter list

tatpulsar.pulse.Crab.retrive_eph.retrieve_ephemeris(write_to_file=True, ephfile='Crab.gro')

retrieve the Jordrell Bank Ephemeris of Crab pulsar return the Pandas DataFrame.

Parameters

write_to_file : boolean

Whether write the retrived ephemeris table into a TXT file

ephfile : str

If write the Table to a TXT file, then assign the TXT file name to write

Returns

df : Pandas DataFrame

Return the Pandas DataFrame

Pulse Analysis

Barycentric Correction

This is a Python based barycentric correction tool. To run the code one needs to have access to one of the ephemeride files provided by jpl: https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/ The default JPL-eph is DE421.

Thanks to Victor, the initial version of this program is composed by V. Doreshenko. Youli made some optimizations. Apply barycentric correction to vector date (in MJD) assuming coordinates ra, dec given in degrees optionally use orbit of a satellinte in heasarc format This includes geometric and shapiro corrections, but einstein correction is not implemented for now. Accurate to about 3e-4s. If return_correction=True (default), correction in seconds is returned. approx einstein is used to define number of points to calculate interpolated einstein correction (0 means do for each data point). see also https://asd.gsfc.nasa.gov/Craig.Markwardt/bary/ for detail and https://pages.physics.wisc.edu/~craigm/idl/ephem.html for the codebase which at least inspired the code below.

tatpulsar.pulse.barycor.barycor.barycor(date, ra, dec, orbit=None, return_correction=False, approx_einstein=10, jplephem=None, accelerate=False)

The core function to compute the barycentric correction.

Parameters

date : array-like

The time series observed by observatory or satellite, (must in MJD)

ra : float

The right ascension of the observed source (in degree)

dec : float

The declination of the observed source (in degree)

orbit : str, default is None

The corresponding orbit file for given time. The corresponding orbit file for the given time. In case of data observed by satellites, the position correction of the satellite orbit needs to be considered. If this parameter is None, the geometric center of the Earth is considered as the position of the periodic signal detection.

Todo

the data from Earth observatory are not supported

return_correction : boolean, default is False

Whether to returns the deviation from the pre-correction data, or the corrected time series. Default is False, which means return the time series.

approx_einstein : int, default is 10

The rate of einstein correction, default is 10

jplephem : File

The JPL solar ephemeris file, default is DE421 stored in the packages. You can assign your required ephemeris, check https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/ for details.

accelerate : boolean, default is False

An accelerate method to boost the computating consumption. The algorithm divided the time series into 60 intervals. And only one time in each interval are computed for barycentric correction. And interpolate the corrected time based on those 60 corrected time.

Binary System Correction

Module for binary related analysis and correction

tatpulsar.pulse.binary.doppler_cor(time, f0, f1, f2, axsini, Porb, omega, e, T_halfpi)

correct the observed freqency of neutron star, convert the frequency moduled by the binary orbital Doppler effect to the intrinsic frequency of NS

Parameters

time : array-like

The epoch of observed time series

f0 : float

the observed frequency at reference epoch t0

f1 : float

the observed frequency derivative at reference epoch t0

f2 : float

the observed frequency second derivative at reference epoch t0

axsini : float

Projected semi-major orbital axis (in units of light-sec)

Porb : float

The period of binary motion (in units of second)

omega : float

Longitude of periastron

e : float

The orbital eccentricity

T_halfpi : float

The mean longitude, with T_halfpi the epoch at which the mean longitude is pi/2 (in units of second)

Returns

f_intri : array-like

The corrected intrinsic frequecy of neutron star

Notes

Reference: Galloway et. al, 2005, ApJ, 635, 1217

tatpulsar.pulse.binary.freq_doppler(time, f0, axsini, Porb, omega, e, T_halfpi)

calculate the frequency modulated by Doppler effect

tatpulsar.pulse.binary.orbit_cor_bt(t, Porb, axsini, e, omega, Tw, gamma)

use numerical method to solve Kepler equation and calculate delay BT model (Blandford & Teukolsky, 1976)

Parameters

t : array-like

The time series without binary correction

Porb : float

The period of binary motion (in units of second)

axsini : float

Projected semi-major orbital axis (in units of light-sec)

e : float

The orbital eccentricity

omega : float

Longitude of periastron

Tw : float

The epoch of periastron passage (in units of seconds, same time system with parameter t)

gamma ; float

the coefficient measures the combined effect of gravitational redshift and time dilation

Returns

new_t : array-like

The time series after the binary correction

Notes

Reference: Blandford & Teukolsky, 1976, ApJ, 205:580-591

tatpulsar.pulse.binary.orbit_cor_deeter(time, Porb, axsini, e, omega, Tnod)

Correct the photon arrival times to the photon emission time Deeter model.

Parameters

time : array-like

The time series before binary correction

Porb : float

The period of binary motion (in units of second)

axsini : float

Projected semi-major orbital axis (in units of light-sec)

e : float

The orbital eccentricity

omega : float

Longitude of periastron

Tnod : float

The epoch of ascending node passage (in units of seconds, same time system with parameter t)

Returns

t_em : array-like

The photon emission time

Notes

Reference: Deeter et al. 1981, ApJ, 247:1003-1012

tatpulsar.pulse.binary.orbit_cor_kepler(time, Tw, ecc, Porb, omega, axsini, PdotOrb=0, omegadot=0, gamma=0)

Corrects observed times of photons to their emission times based on the Kepler function for a binary system.

Parameters

time: float

The time of observed photon in MJD

Tw: float

The periastron time in MJD

ecc: float

Eccentricity

Porb: float

Orbital period in second

PdotOrb: float (optional)

Second derivative of Orbital period in sec/sec default is 0

omega: float

Long. of periastron in radians.

omegadot: float (optional)

second derivative of Long. of periastron in rad/sec

axsini: float

projection of semi major axis in light-sec

gamma: float

0

Returns

t_em: float

The emitting time of photon in binary system in MJD

Examples

>>> observed_time = np.array([...]) # Time series in MJD
>>> Tw = 54424.71               #Barycentric time (in MJD(TDB)) of periastron
>>> e = 0.68                    #orbital eccentricity
>>> Porb = 249.48 * 86400.      #Orbital period (s)
>>> omega = -26. *np.pi/180.    #Longitude of periastron (radians)
>>> axsin = 530.                 #Projected semi-major axis (light seconds)
>>> Porb_dot = 0                #First derivative of Orbital period
>>> binary_corrected_time = orbit_cor_kepler(observed_time
                                            Tw=Tw,
                                            ecc=e,
                                            Porb=Porb,
                                            Porb_dot,
                                            omega,
                                            0,
                                            asin,
                                            0)

Search Periodic Signal

tatpulsar.pulse.search.search(data, **kwargs)

search the best frequency

Parameters

data : array-like

the time series of TDB data

parfile : str, optional

read the parameters from parfile

f0 : float

the init f0 value to search

pepoch : float, optional

time for input frequecy values.

Note

the output frequecy is the frequency at which the middle of each time interval

f0step : float

The step length for frequencies to search

f1 : float, optional

The frequency derivative used to calculate the phase of each photon

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

nbins : int, optional, default is 20

the bin number of profile. default value is 20

f1step : float, optional

The step length to search the frequency derivate as well. Only works if f1step and f1range are both given

f1range : float, optional

the fdot search step range

telescope : str, optional

The name of the mission, support mission are {‘fermi’, ‘hxmt’, ‘nicer’, ‘gecam’, ‘nustar’, ‘ixpe’} parameter determine the MET refenrence time. default is “fermi”.

pepochformat : str, optional

the format of pepoch, “mjd” or “met”. The default if “mjd”

Returns

chi_square : dictionary

The Chi Square distribution of Epoch folding. The valid keys are

• T0: the reference time in MJD

• ChiSquare: the calculated \(\chi^2\) array (return a multi-dimensinal array if both f0 and f1 are searched)

• Profile: the counts of profile folded by best search frequency

• Pars: The dictionary of input and output frequency and frequency derivative {F0, F1, F0_init, F1_init, F2_init, F3_init, F4_init}

Epoch-folding

tatpulsar.pulse.fold.align_profile(profile_list, template)

use ccf function to align each profile in the list to the given template

Parameters

profile_list: list of array

A list of profile. each element in that list is an array of profile

template: array-like

the array of template profile.

Returns

new_list: list

The list of aligned profiles

tatpulsar.pulse.fold.cal_phase(time, pepoch, f0, f1=0, f2=0, f3=0, f4=0, format='met', phi0=0, to_1=True)

calculate the phase for given time or time series.

Parameters

time : array-like

the time series to calculate the phase. The default format of time is “MET”, if your input time is in “MJD” format, set format to ‘mjd’

pepoch : float

time for input frequecy values. The default format of pepoch is “MET”, if your input pepoch is in “MJD” format, set format to ‘mjd’. NOTE: the output frequecy is the frequency at which the time of middle of the time interval

f0 : float

frequency

f1 : float, optional

fdot.

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

format : str, optional, default is ‘met’

The time system of given time and pepoch. Optional input are {‘met’, ‘mjd’}

to_1 : boolean, optional, default is False

normalize phase from 0 to 1

Returns

phase : array-like

The phase value for given time and timing parameters

tatpulsar.pulse.fold.fold(time, parfile=None, pepoch=None, f0=None, f1=0, f2=0, f3=0, f4=0, nbins=20, phi0=0, gti=None, use_data_gti=False, format='met', **kwargs)

Epoch folding the photon array and return the folded profile.

Parameters

time : array-like

the time series of TDB data

parfile : str ,optional

read the parameters from parfile

pepoch : float, optional

time for input frequecy values. NOTE: the output frequecy is the frequency at which the time of middle of the time interval

f0 : float, optional

frequency

f1 : float, optional

fdot.

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

nbins : int, optional, default is 20

the bin number of profile. default value is 20

phi0 : float, optional, default is 0

the reference phase of the profile, if given, the phase is calculated by:

\[\phi = \phi - \phi_{0}\]

gti : N-darray, optional

The good time intervals of given data. default is None, if given, the exposure correction for each phase bin will be considered. You can also use use_data_gti parameter to calculate the GTI for given time array.

use_data_gti: bool, optional, default is False

whether calculate the GTI for given time (if True the time must be an array instead of a value).

format : str, optional

the format of time and pepoch, “mjd” or “met”. The default if “met”.

Warning

The format of event array and the reference time should be the same time format (MJD or MET).

Returns

profile : tatpulsar.data.profile.Profile object

return the profile object define in tatpulsar.data.profile.Profile

tatpulsar.pulse.fold.fold2d(time, y, nseg, parfile=None, pepoch=None, f0=None, f1=0, f2=0, f3=0, f4=0, nbins=20, phi0=0, gti=None, use_data_gti=False, format='met')

Epoch folding the two array into 2 dimensional histogram

Note

histogram2d may be implemented in Profile object.

Parameters

time: array-like

time array

y: array-like

another dimensional information of each photon, such as energy channel, or time

nseg: int

the segment number to uniform split y

parfile : str ,optional

read the parameters from parfile

pepoch : float, optional

time for input frequecy values. NOTE: the output frequecy is the frequency at which the time of middle of the time interval

f0 : float, optional

frequency

f1 : float, optional

fdot.

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

nbins : int, optional, default is 20

the bin number of profile. default value is 20

phi0 : float, optional, default is 0

the reference phase of the profile, if given, the phase is calculated by:

\[\phi = \phi - \phi_{0}\]

gti : N-darray, optional

The good time intervals of given data. default is None, if given, the exposure correction for each phase bin will be considered. You can also use use_data_gti parameter to calculate the GTI for given time array. In fold2d case, gti will apply on each profile extended in y dimension.

use_data_gti: bool, optional, default is False

whether calculate the GTI for given time (if True the time must be an array instead of a value).

format : str, optional

the format of time and pepoch, “mjd” or “met”. The default if “met”.

Warning

The format of event array and the reference time should be the same time format (MJD or MET).

Returns

profiles : list

return a list of Profile object (tatpulsar.data.profile.Profile)

Raises

IOError

If pepoch, f0, and parfile not given, or if data is empty.

tatpulsar.pulse.fold.fold_lightcurve(time, counts, pepoch, f0, nbins=16, dt=1, counts_err=None, f1=0, f2=0, f3=0, f4=0, format='met', **kwargs)

Fold the pulse profiles from the net light curves.

Parameters

time : array-like

the time series of TDB data

counts : array-like

the counts of the light curve

pepoch : float, optional

time for input frequecy values. NOTE: the output frequecy is the frequency at which the time of middle of the time interval

f0 : float, optional

frequency

nbins : int

the bin number of profile. default value is 16

dt : float

the bin size of the light curve, default value is 1 s

counts_err : array-like

the counts error of the light curve

f1 : float, optional

fdot.

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

format : str, optional

the format of time and pepoch, “mjd” or “met”. The default if “met”.

Returns

profile : tatpulsar.data.profile.Profile object

return the profile object define in tatpulsar.data.profile.Profile

tatpulsar.pulse.fold.merge_aligned_profile(profile_list, template)

align each profile in the profile_list and return the profile that sum up the counts of each profile

Parameters

profile_list: list of array

A list of profile. each element in that list is an array of profile

template: array-like

the array of template profile.

Returns

new_profile: array

The array of merged aligned profiles

tatpulsar.pulse.fold.phase_exposure(gti, nbins, f0, f1=0, f2=0, f3=0, f4=0, pepoch=0, phi0=0, format='met')

calculate the exposure correction coefficients for each phase bin.

Parameters

gti: ndarray or list

the list of GTI array, example [[gti0_0, gti0_1], [gti1_0, gti1_1], …]

nbins : int

the bin number of profile

f0 : float

frequency

f1 : float, optional

fdot.

f2 : float, optional

the second derivative of frequency

f3 : float, optional

the third derivative of frequency

f4 : float, optional

the fourth derivative of frequency

pepoch : float, optional

time for input frequecy values. NOTE: the output frequecy is the frequency at which the time of middle of the time interval

format : str, optional

the format of time and pepoch, “mjd” or “met”. The default if “met”.

Returns

expo_corr: array-like

the coeffients for each phase bin to calculate the exposure correction. to correction the count in each phase bin, the count (\(R(\phi)\)) should divide the corresponding expo_corr value (\(c(\phi)\)).

Timing Analysis

tatpulsar.pulse.residuals.cal_residual(toas, toa_errs, f_set_all, PEPOCH_all, start_time, stop_time, inperiod=False, phi_shift=0)

calculate the residuals for toas in each Ephemeride

return the residuals as one

tatpulsar.pulse.residuals.cal_residual_from_parameters(toas, toa_errs, F_set_array, PEPOCH, inperiod=False)

calculate the residuals for toas in one set of Ephemeris parameters return the residuals

Parameters

toas : array-like

ToAs in MJD unit

toa_errs : array-float

ToAs error in microsecond unit

F_set_array : array-float

An array of f0, f1, f2, … etc.

PEPOCH : float

reference time in MJD unit

inperiod : bool

flag to calculate residuals in period unit

Returns

residuals : array-like

The timing residuals in either period or time

residual_err : array-like

The error of timing residuals

residuals_rms : array-like

The root-mean-square of total ToA set from the weighted mean ToA values:

\[\sqrt{\sum \frac{(T - <T>)^2}{N-1}},\]

where \(T\) is the ToAs, and \(<T>\) is the weighted mean value of ToAs

tatpulsar.pulse.residuals.read_toa(timname)

read the standard TEMPO2 format *.tim ToA file and return the ToA and ToA errors

Simulation

functions for the profile sampling and simulation

tatpulsar.simulation.profile_sim.draw_event_from_phase(phase, tstart, tstop, f0, f1=0, f2=0, f3=0, f4=0, pepoch=0)

sampling the photon arrival time given the phase of each photon

Parameters

phase: array-like

The phase of each photon, the phase is normalized to [0, 1]

tstart: array-like

The start time to generate arrival time (MJD)

tstop: array-like

The stop time to generate arrival time (MJD)

f0: float

The frequency of the pulsar

f1: float, optional

The frequency derivative of the pulsar

f2: float, optional

The second derivative of the frequency

f3: float, optional

The third derivative of the frequency

f4: float, optional

The fourth derivative of the frequency

pepoch: float

The reference time of the timing parameters (MJD)

Returns

event_list: array-like

The sampled arrival times of photons

tatpulsar.simulation.profile_sim.poisson_rejection_sampling(x, y, nphot)

rejection sampling the poisson like signal x-y

Parameters

x: array-like

Time or Phase array

y: array-like

counts, Poisson distributed

nphot: int

The output amount of sampled photons

Returns

xval: array-like

The x sample that satisfy the rejection rule