Get time delays using PINT

time_delay.py

"""
Get solar system barycentring time delays using PINT and lalsuite.
(see the example at https://github.com/nanograv/PINT/blob/master/examples/TimingModel_composition.ipynb)
"""

import subprocess as sp

# create a fake pulsar (units must be TDB as PINT can't do TCB yet!)
parcontent = """PSRJ   J0123+4501
RAJ    01:23:00.0
DECJ   45:01:00.0
F0     123.456789
F1     -1.23456e-12
PEPOCH 56000
EPHEM  DE405
UNITS  TDB
"""

parfile = 'J0123+4501.par'

with open(parfile, 'w') as fp:
    fp.write(parcontent)


# use TEMPO2 to generate some simlated TOAs
p = sp.Popen('tempo2 -gr fake -f {} -start 55000 -end 55365 -rms 0.00006 -ndobs 1 -nobsd 1 -randha y'.format(parfile),
             shell=True)
out, err = p.communicate()

# read in the par file and TOAs with PINT
from pint.models import get_model
m = get_model(parfile)

print(m.DelayComponents_list)

# remove any component that are not AstrometryEquatorial or SolarSystemShapiro if necessary e.g.
# > component, order, from_list, comp_type = m.map_component('DispersionDM')
# > from_list.remove(component)

# set planet Shapiro delays to zero
component, order, from_list, comp_type = m.map_component('SolarSystemShapiro')
component.PLANET_SHAPIRO.value = False

from pint.toa import get_TOAs

t = get_TOAs('J0123+4501.simulate', ephem='DE405', include_gps=True)

# get the delays
total_delay = m.delay(t)

# functions for using LAL
def detector2ssb(tseries, detector, ra, dec, ephem='DE405', units='TDB'):
    """
    Calculate the barycentring time delay between detector and SSB.

    Args:
        tseries (array_like): a times series of GPS time of arrivals at the SSB
        detector (str): a detector name
        ra (float): a source right ascension (rads)
        dec (float): a source declination (rads)
        ephem (str): the JPL solar system ephemeris type (default is DE405)
        units (str): the coordinate system units (TDB or TCB)

    Returns:
        An array of time delay values.
    """

    if ephem not in ['DE200', 'DE405', 'DE421', 'DE430']:
        raise ValueError("Ephemeris '{}' not available".format(ephem)

    # load ephemeris files (download if necessary)
    earthephem = 'earth00-40-{}.dat.gz'.format(ephem)
    sunephem = 'sun00-40-{}.dat.gz'.format(ephem)

    # check if files are in a known path
    earthfcheck = lalp.PulsarFileResolvePath(earthephem)

    if earthfcheck is None:
        # try downloading and caching the file (using astropy)
        try:
            from astropy.utils.data import download_file
        except ImportError:
            raise ImportError("Could not import astropy")

        earthurl = 'https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/{}'.format(earthephem)

        # download the file (or use the cached file)
        try:
            earthephem = download_file(earthurl, cache=True)
        except Exception as e:
            raise e("Problem downloading Earth ephemeris file")
    
    sunfcheck = lalp.PulsarFileResolvePath(sunephem)

    if sunfcheck is None:
        # try downloading and caching the file (using astropy)
        try:
            from astropy.utils.data import download_file
        except ImportError:
            raise ImportError("Could not import astropy")

        sunurl = 'https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/{}'.format(sunephem)

        # download the file (or use the cached file)
        try:
            sunephem = download_file(sunurl, cache=True)
        except Exception as e:
            raise e("Problem downloading Sun ephemeris file")

    # load ephemeris files
    try:
        edat = lalp.InitBarycenter(earthephem, sunephem)
    except IOError:
        raise IOError("Problem loading ephemeris files")

    # get time units file
    if units == 'TDB':
        timefile = 'tdb_2000-2040.dat.gz'
        ttype = 1
    elif units == 'TCB':
        timefile = 'te405_2000-2040.dat.gz'
        ttype = 2
    else:
        raise ValueError("Units must be 'TDB' or 'TCB'")
    
    timefcheck = lalp.PulsarFileResolvePath(timefile)
    
    if timefcheck is None:
        # try downloading and caching the file (using astropy)
        try:
            from astropy.utils.data import download_file
        except ImportError:
            raise ImportError("Could not import astropy")

        timeurl = 'https://git.ligo.org/lscsoft/lalsuite/raw/master/lalpulsar/src/{}'.format(timefile)

        # download the file (or use the cached file)
        try:
            timefile = download_file(timeurl, cache=True)
        except Exception as e:
            raise e("Problem downloading time file")

    # load ephemeris files
    try:
        tdat = lalp.InitTimeCorrections(timefile)
    except IOError:
        raise IOError("Problem loading time file")
                         
    # set the site information
    site = lalp.GetSiteInfo( detector.upper() )

    tdelay = np.zeros(len(tseries))

    # make sure sky positions are in the allowed range
    rat, dect = lal.NormalizeSkyPosition(ra, dec)
    
    # initialise variables
    bary = lalp.BarycenterInput()
    bary.site.location = site.location/lal.C_SI

    bary.alpha = rat
    bary.delta = dect
    bary.dInv = 0.  # assume no parallax

    earth = lalp.EarthState()
    emit = lalp.EmissionTime()
                         
    for i, ttime in enumerate(tseries):
        # convert time to LIGOTimeGPS structure
        bary.tgps = lal.LIGOTimeGPS(ttime)

        lalp.BarycenterEarthNew( earth, gps, ephem, tephem, ttype[units] )
        lalp.Barycenter( emit, bary, earth )

        # get the time delay between detector and SSB
        tdelay[i] = emit.deltaT

    return tdelay