GammaLib

#! /usr/bin/env python

# ==========================================================================

# Reimplements COMPASS taks PULPHI

#

# Copyright (C) 2022 Juergen Knoedlseder

#

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

#

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

# GNU General Public License for more details.

#

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

# ==========================================================================

import math

import gammalib

# Static stuff for bvceph

jdch0  = 0

earth  = []

sun    = []

tdbtdt = []

tdtut  = []

tdbdot = []

# ======================================================= #

# Extrapolate pulsar frequency, derivative and zero phase #

# ======================================================= #

def pula32(f0, fdot, f2dot, zphas, t0, tm):

    """

    Implement AL-032 to extapolate the pulsar frequency, derivative and

    zero phase to the center of the epoch

    """

    # Set constants

    pow2 = 2147483647.0

    # Calculate time difference in seconds

    tdif = float(tm[0] - t0[0]) * 86400.0 + float(tm[1] - t0[1]) / 8000.0

    # Extrapolate frequency, all times are in seconds

    f0e = f0 + tdif * fdot + (tdif*tdif * f2dot) / 2.0

    # Extrapolate frequency derivative

    fdote = fdot + tdif * f2dot

    # Extrapolate the zero phase

    x = zphas + f0*tdif + ((fdot  * math.pow(tdif,2)) / 2.0) + \

                          ((f2dot * math.pow(tdif,3)) / 6.0)

    # Accuracy checks

    if x > pow2:

        npow = int(x / pow2)

        x    = x - (float(npow)*pow2)

    # Calculate the residual by removing the integer part

    zphasp = x - float(int(x))

    # Fix if < 0

    if zphasp < 0.0:

        zphasp = 1.0 + zphasp

    # Return

    return (f0e, fdote, zphasp)

# =================== #

# Get BVC information #

# =================== #

def pbvcrr(bvc, irow):

    """

    Get BVC information

    """

    # Extract information for row irow

    time = [bvc['TJD'].integer(irow), bvc['TICS'].integer(irow)]

    ssb  = [bvc['SSB_X'].real(irow), bvc['SSB_Y'].real(irow), bvc['SSB_Z'].real(irow)]

    etut = bvc['TDELTA'].real(irow)

    # Return information

    return (time, ssb, etut)

# =================== #

# Get EVP information #

# =================== #

def pevpsr(evp, irow):

    """

    Get EVP information

    """

    # Extract information for row irow

    time = [evp['TJD'].integer(irow), evp['TICS'].integer(irow)]

    # Return information

    return (time)

# ========================== #

# Evaluates the GRO position #

# ========================== #

def pulbvc(evtime, tag1, tag2, y1, y2):

    """

    Evaluates the GRO position at the time of the event by interpolating the

    positions supplied every 16 sec. by BVC. The interpolating function is

    a 1st order curve :  Y = A * X + B

    """

    # Compute times

    xtime1 = float(tag1[0])   * 86400.0 + float(tag1[1])   / 8000.0

    xtime2 = float(tag2[0])   * 86400.0 + float(tag2[1])   / 8000.0

    etime  = float(evtime[0]) * 86400.0 + float(evtime[1]) / 8000.0

    dt     = xtime2 - xtime1

    #

    a = [(y2[0] - y1[0])/dt,

         (y2[1] - y1[1])/dt,

         (y2[2] - y1[2])/dt]

    #

    ssb = [a[0] * (etime - xtime1) + y1[0],

           a[1] * (etime - xtime1) + y1[1],

           a[2] * (etime - xtime1) + y1[2]]

    # Return

    return ssb

# ============================================= #

# Convert event time to solar system Barycentre #

# ============================================= #

def pula04(ftime, ssb, pos1, pos2, pos3, etut):

    """

    This routine implements the algorithm PUL-AL-004 to calculate the arrival times of

    detected photons at the Solar System Barycentre (SSB)

    """

    # Set constants

    radius = 0.49254909e-5

    # Find the light travel time from the satellite to SSB along the pulsar direction

    travt = (pos1*ssb[0] + pos2*ssb[1] + pos3*ssb[2]) * 1.e-6

    # Find the relativistic delay due to the sun

    r     = 1.e-6 * (math.sqrt((ssb[0]*ssb[0] + ssb[1]*ssb[1] + ssb[2]*ssb[2])))

    relat = -2.0 * radius * math.log(1.0 + (travt / r))

    # For debugging

    tdelta = travt - relat + etut

    # Find the time at SSB

    ssbtim = [ftime[0], ftime[1]+int(((travt - relat + etut) * 8000.0) + 0.5)]

    # Handle overflow and underflow

    if ssbtim[1] >= 691200000:

       ssbtim[1] = ssbtim[1] - 691200000

       ssbtim[0] = ssbtim[0] + 1

    elif ssbtim[1] < 0:

       ssbtim[1] = ssbtim[1] + 691200000

       ssbtim[0] = ssbtim[0] - 1

    # Return

    return ssbtim, tdelta

# ====================== #

# Compute residual phase #

# ====================== #

def pula05(ctime, f, fdot, f2dot, dm, fobs, zphas, tm):

    """

    This routine implements the algorithm PUL-AL-005 to calculate the residual phase

    corresponding to each photon arrival time at SSB.

    """

    # Set constants

    pow2 = 2147483647.0

    # Make all in double precision

    dmd   = float(dm)

    fobsd = float(fobs)

    # Apply dispersion measure

    if fobsd > 1.0e-6:

        the0 = (-dmd * 10000.0) / (2.41 * fobsd*fobsd)

    else:

        the0 = 0.0

    #

    dt = float(ctime[0] - tm[0]) * 86400.0 + float(ctime[1] - tm[1]) / 8000.0

    # Theta is the decimal residual of the formula below

    x = the0 + (f * dt) + (fdot * dt*dt / 2.0) + (f2dot * dt*dt*dt / 6.0) + zphas

    # Accuracy checks

    if abs(x) > pow2:

        npow = int(x / pow2)

        x    = x - (float(npow)*pow2)

    # Calculate the residual by removing the integer part

    resid = x - float(int(x))

    # Fix if < 0

    if resid < 0.0:

        resid = 1.0 + resid

    # Return

    return resid

# ===================== #

# Returns JPL ephemeris #

# ===================== #

def bvceph(jdutc, frc, nleap, filename='de200_new.fits'):

    """

    Reads and interpolates ephemeris files created from JPL data by the CONVERT

    program, returning positions of the Earth and Sun with respect to the Solar

    System barycenter (SSBC) and the difference between Ephemeris Time (TDB)

    and Universal Time (UTC) at the same epoch.

    """

    # Declare global variables

    global jdch0

    global earth

    global sun

    global tdbtdt

    global tdtut

    global tdbdot

    # Set constants

    ndaych = 16

    daysec = 1.0 / 86400.0

    xxx    = 10.0343817

    etatc  = 32.184 - 0.0343817

    # Choose set of coefficents to use; if out of range, read new chunk

    nset = 1 + jdutc - jdch0

    if nset > ndaych or nset < 1:

        # Open table

        fits  = gammalib.GFits(filename)

        ephem = fits.table('EPHEM')

        # If this is the first call, open the file for direct access and

        if jdch0 == 0:

            jd0 = int(ephem.real('TSTART'))

            jd1 = int(ephem.real('TSTOP'))

        # Read earth and sun positions and tdb-tdt from the ephemeris

        # Get at-ut from the leap second table in a1utcf

        jdch1  = min([jdutc+(ndaych-1),jd1])

        jdch0  = jdch1 - (ndaych-1)

        earth  = []

        sun    = []

        tdbtst = []

        tdtut  = []

        tdbdot = []

        for nset in range(ndaych):

            # Get record index (starting from 0)

            nrec = nset + (jdch0-jd0)

            #print(nrec)

            # Read information for one row

            earth.append({'earth_x'   : ephem['EARTH'][nrec,0],

                          'earth_y'   : ephem['EARTH'][nrec,1],

                          'earth_z'   : ephem['EARTH'][nrec,2],

                          'earth_dx'  : ephem['EARTH'][nrec,3],

                          'earth_dy'  : ephem['EARTH'][nrec,4],

                          'earth_dz'  : ephem['EARTH'][nrec,5],

                          'earth_d2x' : ephem['EARTH'][nrec,6],

                          'earth_d2y' : ephem['EARTH'][nrec,7],

                          'earth_d2z' : ephem['EARTH'][nrec,8],

                          'earth_d3x' : ephem['EARTH'][nrec,9],

                          'earth_d3y' : ephem['EARTH'][nrec,10],

                          'earth_d3z' : ephem['EARTH'][nrec,11]})

            sun.append({'sun_x' : ephem['SUN'][nrec,0],

                        'sun_y' : ephem['SUN'][nrec,1],

                        'sun_z' : ephem['SUN'][nrec,2]})

            tdbtdt.append(ephem['TIMEDIFF'][nrec])

            # Compute tdtut

            tdtut.append(etatc + nleap + xxx)

        # Close FITS file

        fits.close()

        # Make rough computation of time derivative of tdb-tdt

        for nset in range(ndaych-1):

            tdbdot.append(tdbtdt[nset+1]-tdbtdt[nset])

        tdbtdt[ndaych-1] = tdbtdt[ndaych-2]

        # Set nset (starting from 1)

        nset = 1 + (jdutc - jdch0)

    # DEBUG

    #print(nset-1)

    # Now use taylor series to get the coordinates at the desired time

    dt   = frc + tdtut[nset-1] * daysec

    dt2  = dt  * dt

    dt3  = dt2 * dt

    etut = tdbtdt[nset-1] + dt * tdbdot[nset-1] + tdtut[nset-1]

    rcs  = [sun[nset-1]['sun_x'], sun[nset-1]['sun_y'], sun[nset-1]['sun_z']]

    rce  = [earth[nset-1]['earth_x'] + \

            earth[nset-1]['earth_dx']  * dt + \

            earth[nset-1]['earth_d2x'] * 0.5 * dt2 + \

            earth[nset-1]['earth_d3x'] * 1.0/6.0 * dt3,

            earth[nset-1]['earth_y'] + \

            earth[nset-1]['earth_dy']  * dt + \

            earth[nset-1]['earth_d2y'] * 0.5 * dt2 + \

            earth[nset-1]['earth_d3y'] * 1.0/6.0 * dt3,

            earth[nset-1]['earth_z'] + \

            earth[nset-1]['earth_dz']  * dt + \

            earth[nset-1]['earth_d2z'] * 0.5 * dt2 + \

            earth[nset-1]['earth_d3z'] * 1.0/6.0 * dt3]

    vce  = [(earth[nset-1]['earth_dx'] + earth[nset-1]['earth_d2x']  * dt) * daysec,

            (earth[nset-1]['earth_dy'] + earth[nset-1]['earth_d2y']  * dt) * daysec,

            (earth[nset-1]['earth_dz'] + earth[nset-1]['earth_d2z']  * dt) * daysec]

    # DEBUG

    #print(dt)

    # Return

    return rce, etut

# =========== #

# pulphi task #

# =========== #

def pulphi(pars, evpname, bvcname, maxevents=10000000):

    """

    Generates barycentric times and pulse phases for a given pulsar and

    outputs events lying in the given phase windows.

    """

    # Set constants

    zztic  = 691200000

    pimath = 3.14159265358979

    # Read parameters

    # CALL PPPTRD ( TPAR, PUNAME, PPTUN, EPOCH, PPT0, PPT1, PPT2, RAPU,

    #1              DECPU, FOBS, DM, DISPU, BIPE, BIPED, BIOE, SEMA,

    #2              PERI, PERID, TPERI, ZPHAS, G, RETGP, RETNO )

    # Derive some items for printing and write the contents of the PPT

    # and derived items in log file

    timejd = float(pars['epoch'][0]) + float(pars['epoch'][1]) / (8000.0 * 86400.0) + 2440000.5

    timemj = timejd - 2400000.5

    print('%s%d.%d' % (gammalib.parformat('Pulsar epoch'), pars['epoch'][0], pars['epoch'][1]))

    print('%s%f' % (gammalib.parformat('Pulsar epoch JD'), timejd))

    print('%s%f' % (gammalib.parformat('Pulsar epoch MJD'), timemj))

    # Convert frequency to true units

    f0    = pars['ppt0']

    fdot  = pars['ppt1'] * 1.0e-12

    f2dot = pars['ppt2'] * 1.0e-20

    print('%s%f' % (gammalib.parformat('Frequency'), f0))

    print('%s%e' % (gammalib.parformat('Frequency derivative'), fdot))

    print('%s%e' % (gammalib.parformat('Second freq. derivative'), f2dot))

    # Write rest to log file

    print('%s%f' % (gammalib.parformat('Right Ascension'), pars['rapu']))

    print('%s%f' % (gammalib.parformat('Declination'), pars['decpu']))

    print('%s%f' % (gammalib.parformat('Radio frequency'), pars['fobs']))

    print('%s%f' % (gammalib.parformat('Dispersion measure'), pars['dm']))

    print('%s%f' % (gammalib.parformat('Pulsar distance'), pars['dispu']))

    print('%s%f' % (gammalib.parformat('Radio phase'), pars['zphas']))

    # Open EVP dataset

    evpfile = gammalib.GFits(evpname)

    evp     = evpfile.table('COMPTEL_EVP')

    #print(evp.header())

    # Calculate the centre time of the current observation

    intsta = [evp.integer('VISDAY'), evp.integer('VISTIM')]

    intend = [evp.integer('VIEDAY'), evp.integer('VIETIM')]

    tm     = [int(intsta[0]+intend[0]/2), int(intsta[1]+intend[1]/2)]

    # Check to make sure half a day is not lost

    if 2*tm[0] == (intsta[0] + intend[0] - 1):

        tm[1] = tm[1] + 345600000

        if tm[1] >= zztic:

            tm[1] = tm[1] - zztic

            tm[0] = tm[0] + 1

    # Centre in JD

    tmjd = float(tm[0]) + float(tm[1]) / (8000.0 * 86400.0) + 2440000.5

    # Centre in MJD

    tmmj = tmjd - 2400000.0

    # Calculate the length of the current observation

    dt = (float(intend[0] - intsta[0]) * 86400.0) + (float(intend[1] - intsta[1]) / 8000.0)

    # Use AL-032 to extapolate the pulsar frequency, derivative and zero phase to the

    # center of the epoch

    f0e, fdote, zphasp = pula32(f0, fdot, f2dot, pars['zphas'], pars['epoch'], tm)

    # Calculate the cosinus of the pulsar direction

    ra   = pars['rapu']  * (pimath/180.0)

    d    = pars['decpu'] * (pimath/180.0)

    pos1 = math.cos(d) * math.cos(ra)

    pos2 = math.cos(d) * math.sin(ra)

    pos3 = math.sin(d)

    # Open the BVC dataset for SSB vectors

    bvcfile = gammalib.GFits(bvcname)

    bvc     = bvcfile.table('COMPTEL_BVC')

    #print(bvc.header())

    # Get the first 2 vectors

    start, ssb1, etut = pbvcrr(bvc, 0)

    endt,  ssb2, etut = pbvcrr(bvc, 1)

    # Open ASCII result file

    f = open('pulphi.dat', 'w')

    # Main event processing loop

    ievp   = 0

    ibvc   = 2

    nevin  = 0

    nevout = 0

    scan   = True

    while scan:

        # Read event time

        fftime = pevpsr(evp, ievp)

        ievp += 1

        # Increment number of input events

        nevin = nevin + 1

        # Make time correction if necessary

        if (fftime[0] < 8798) or ((fftime[0] == 8798) and (fftime[1] < .28800000)):

            ftime = [fftime[0], fftime[1] - (2.042144 * 8000.0)]

            if ftime[1] < 0.0:

                ftime[0] = ftime[0] - 1

                ftime[1] = ftime[1] + zztic

        else:

            ftime = [fftime[0], fftime[1]]

        # While the event time is after the end of validity of the current

        # SSB vector get the next one.

        while True:

            # If arrival time is < first time in BVC

            if (ftime[0] < start[0]) or ((ftime[0] == start[0]) and (ftime[1] < start[1])):

                tdif = float(ftime[0] - start[0]) * 86400.0 + float(ftime[1] - start[1]) / 8000.0

                print('<<< %d:%d %e %d:%d' % (ftime[0],ftime[1],tdif,start[0],start[1]))

                ssb = [ssb1[0],ssb1[1],ssb1[2]]

                break

            # If arrival time is within the interval supplied (START < t < ENDT)

            elif ((ftime[0] > start[0]) or ((ftime[0] == start[0]) and (ftime[1] > start[1]))) and \

                 ((ftime[0] <  endt[0]) or ((ftime[0] ==  endt[0]) and (ftime[1] <= endt[1]))) :

                tdif = float(endt[0] - start[0]) * 86400.0 + float(endt[1] - start[1]) / 8000.0

                if tdif > 900.0:

                    print('=== %d:%d %e' % (ftime[0],ftime[1],tdif,start[0],start[1]))

                ssb = pulbvc(ftime, start, endt, ssb1, ssb2)

                break

            # If arrival time is > ENDT

            elif (ftime[0] > endt[0]) or ((ftime[0] == endt[0]) and (ftime[1] > endt[1])):

                # Replace start by endt

                start[0] = endt[0]

                start[1] = endt[1]

                ssb1[0]  = ssb2[0]

                ssb1[1]  = ssb2[1]

                ssb1[2]  = ssb2[2]

                # Was end of BVC file reached?

                if ibvc >= bvc.nrows():

                    print('*** End of BVC file reached. Stop')

                    scan = False

                    break

                # Read new BVC vector

                endt, ssb2, etut = pbvcrr(bvc, ibvc)

                ibvc += 1

        # Was end of BVC file reached?

        if not scan:

            break

        # Correct the photon arrival time using SSB vector 1

        ssbtim, tdelta = pula04(ftime, ssb, pos1, pos2, pos3, etut)

        # Do not consider binary

        ctime = [ssbtim[0],ssbtim[1]]

        # Calculate phase residuals

        resid = pula05(ctime, f0e, fdote, f2dot, pars['dm'], pars['fobs'], zphasp, tm)

        # Write result in ASCII file

        f.write('%8d  %4.4d:%9.9d  %4.4d:%9.9d  %4.4d:%9.9d  %14.8f  %10.8f\n' %

               (ievp-1, fftime[0], fftime[1], ftime[0], ftime[1], ctime[0], ctime[1], tdelta, resid))

        # Was end of EVP file reached?

        if ievp >= evp.nrows():

            print('*** End of EVP file reached. Stop')

            scan = False

            break

        # Was maximum number of events reached?

        if ievp >= maxevents:

            print('*** Maximum number of events reached. Stop')

            scan = False

            break

    # Close ASCII file

    f.close()

    # Return

    return

# ================================ #

# Convert parameters into PPT file #

# ================================ #

def pulcpp(pars):

    """

    Convert parameters into PPT file

    """

    # Set constants

    mj1968 = 2400000

    daysec = 24.0 * 3600.0

    convd  = 4.0 * math.atan(1.0) / 180.0

    rsw    = 0.49254909e-5

    factor = [1.0, 0.5, 1.0/6.0]

    # Convert EPOCH(1) from TJD to MJD

    epoch = [int(pars['t0geo']) - 40000.0, 0]

    #tgeoin = int(pars['epoch'][0] + 40000)

    tgeoin = int(pars['t0geo'])

    print('%s%d' % (gammalib.parformat('Integer part Tgeo (MJD)'), tgeoin))

    #

    #dgeofr = pars['tgeofr']

    dgeofr = pars['t0geo'] - tgeoin

    print('%s%e' % (gammalib.parformat('Fraction part Tgeo (MJD)'), dgeofr))

    #

    tpkcor = pars['tpkcor']

    print('%s%f' % (gammalib.parformat('Correction to Tgeo (ms)'), tpkcor))

    tpkcr8 = float(tpkcor) * 1.0e-3 / daysec

    print('%s%f' % (gammalib.parformat('Correction to Tgeo (days)'), tpkcr8))

    #

    jd = int(tgeoin + mj1968 + 1)

    #

    aid1 = dgeofr

    aid2 = tpkcr8

    aid3 = 0.5

    #

    frc = aid1 + aid2 - aid3

    print('%s%d' % (gammalib.parformat('Calculated JD (int)'), jd))

    print('%s%f' % (gammalib.parformat('Calculated JD (frac)'), frc))

    #

    f     = pars['ppt0']

    fdot  = pars['ppt1'] * 1.0e-12

    fddot = pars['ppt2'] * 1.0e-20

    # Compute number of leap seconds (-10)

    nleap = pars['nleap']

    print('%s%d' % (gammalib.parformat('Number of Leap seconds'), nleap))

    #

    freqt = [f, fdot, fddot]

    print('%s%f' % (gammalib.parformat('Frequency'), freqt[0]))

    print('%s%e' % (gammalib.parformat('Frequency derivative'), freqt[1]))

    print('%s%e' % (gammalib.parformat('Second freq. derivative'), freqt[2]))

    #

    decdp = pars['decpu'] * convd

    radp  = pars['rapu']  * convd

    #

    pos1 = math.cos(decdp) * math.cos(radp)

    pos2 = math.cos(decdp) * math.sin(radp)

    pos3 = math.sin(decdp)

    print('%s%f' % (gammalib.parformat('Pos1'), pos1))

    print('%s%f' % (gammalib.parformat('Pos2'), pos2))

    print('%s%f' % (gammalib.parformat('Pos3'), pos3))

    #

    print('%s%f' % (gammalib.parformat('Pulsar position RA 2000'), radp/convd))

    print('%s%f' % (gammalib.parformat('Pulsar position DEC2000'), decdp/convd))

    #

    rce, etut0 = bvceph(jd, frc, nleap)

    #

    ctest = nleap + 10 + 32.184

    etut  = etut0

    dtest = etut - ctest

    #

    ssb = [rce[0], rce[1], rce[2]]

    # Calculate traveltime from SSB to site (0,0,0) for source at (POSX,POSY,POSZ)

    bclt   = ssb[0]*pos1 + ssb[1]*pos2 + ssb[2]*pos3

    ssbnrm = math.sqrt(ssb[0]*ssb[0] + ssb[1]*ssb[1] + ssb[2]*ssb[2])

    # Relativistic delay

    relat = -2.0 * rsw * math.log(1.0 + bclt / ssbnrm)

    # Time difference (in sec) between Arrival at SSB and nominal epoch

    delay = bclt - relat + etut

    # Log

    print('%s%f sec' % (gammalib.parformat('UTC--> TDT'), ctest))

    print('%s%f sec' % (gammalib.parformat('TDT--> TDB'), dtest))

    print('%s%f sec' % (gammalib.parformat('UTC--> TDB'), etut0))

    print('%s%f sec' % (gammalib.parformat('Traveltime difference'), bclt))

    print('%s%f sec' % (gammalib.parformat('Relativistic delay'), relat))

    print('%s%f sec' % (gammalib.parformat('Total delay'), delay))

    #

    difsec = (dgeofr + tpkcr8) * daysec + delay

    print('%s%f sec' % (gammalib.parformat('Time diff. to nom. epoch'), difsec))

    #

    radiof = 0.0

    dlyaid = 1.0

    for i in range(3):

        dlyaid = dlyaid * difsec

        radiof = radiof + factor[i] * freqt[i] * dlyaid

    #

    radiof = radiof - int(radiof)

    #

    if radiof < 0.0:

        radiof = radiof + 1.0

    #

    zphas = radiof

    print('%s%f' % (gammalib.parformat('Radio-phase 0 at nom.epoch'), zphas))

    #

    zphas = -zphas

    # Set parameters

    pars['epoch'] = epoch

    pars['zphas'] = zphas

    # Return

    return pars

# ======================== #

# Main routine entry point #

# ======================== #

if __name__ == '__main__':

    # Line 53 of psrtime file

    #0531+21  05 34 31.972  22 00 52.07 48377 48407 48392.000000224  29.9485663280898 -3.77641D-10   0.00D+00  1.5 J

    # Set pulsar parameters

    #pars = {'epoch': [48392, 0],

    #        'tgeofr':   0.000000224,

    pars = {'t0geo'  : 48392.000000224,

            'tpkcor' :     0.0,

            'ppt0'   :    29.9485663280898,

            'ppt1'   :  -377.641,

            'ppt2'   :     0.0,

            'rapu'   :    83.6331,

            'decpu'  :    22.0145,

            'fobs'   :     0.0,

            'dm'     :     0.0,

            'dispu'  :     0.0,

            'nleap'  :      16} # valid between 1990-12-31 and 1992-06-30

    # Set filenames

    evpname = '$COMDATA/phase01/vp0001_0/m28511_evp.fits'

    bvcname = '$COMDATA/phase01/vp0001_0/s10150_bvc.fits'

    # Compute zphas using the pulcpp task

    pars = pulcpp(pars)

    # Run task

    pulphi(pars, evpname, bvcname, maxevents=100)