GammaLib

Note that BVC files were produced using BVCGEN and are available in the HEASARC archive. The BVC files contain the CGRO spacecraft (X,Y,Z) position in km with respect to the Sun Barycentre as well as the time delay in seconds (TDELTA column) as function of COMPTEL TJD and TICS. Entries are provided every 131072 tics or 16.384 sec, which is the frequency of a superpacket (i.e. the same frequency as the OAD datasets.

For a given VP all BVC information is in a single file. It’s probably sufficient to implement a GCOMBvc class to deal with the information from the file. Probably the only thing that is really needed is a function that returns the time delay (TDELTA column) as function of COMPTEL TJD and TICS. The question is whether an interpolation of the time delays is needed, which would imply that we need to know to which moment within a superpacket the information corresponds.

Here an interesting note on pulsar ephemerides in the CGRO epoch: https://www.jb.man.ac.uk/~pulsar/crab/CGRO_format.html.

After much discussion of the issues, Dick Manchester and I
have settled on a specific format for data files containing
the information you'll need.  Most of the essential information
will be in a file called "psrtime.dat", available by anonymous ftp
from puppsr.princeton.edu (or Internet address 128.112.128.165).
After you login with username "anonymous" and using anything for
a password, type "cd gro" to change to the GRO directory.
Then you can "get psrtime.dat" to download the main data file.
A copy of this message (or suitably updated instructions) will
be available as the file "read.me", and the file containing binary
pulsar data is "psrbin.dat".

A short dummy version of the main file looks like this (you'll have
to list it on a wide or "landscape mode" (120 column) printer to
see the full width):

 PSR B     RA(J2000)   DEC(J2000)   MJD1  MJD2    t0geo(MJD)         
f0(s^-1)      f1(s^-2)     f2(s^-3)  RMS O B
-------------------------------------------------------------------------
----------------------------------------
0531+21  05 34 31.974  22 00 52.05 48282 48316 48299.000000134  
29.9516010684378 -3.77726D-10   0.00D+00  5.5 P 
0655+64  07 00 37.811  64 18 11.26 44546 47829 46187.000002101   
5.1106207921632 -1.76088D-17   0.00D+00  0.7 P *
0833-45  08 35 20.680 -45 10 35.70 48256 48320 48288.000000968  
11.1989616062225 -1.55846D-11   1.08D-21  0.5 P 
1855+09  18 57 36.394  09 43 17.33 46437 48239 47338.000000045 
186.4940816803072 -6.20316D-16   0.00D+00  0.2 P *

The column headings should be reasonably self-explanatory:

PSR    Pulsar name in 1950 coordinates (for consistency with prior usage)
RA    Right Ascension in J2000 coordinates (hh mm ss.sss)
DEC    Declination in J2000 coordinates (sdd mm ss.ss)
MJD1,2    First and last dates for valid parameters (MJD)
t0geo    Infinite-frequency geocentric UTC arrival time of a pulse (MJD)
    Note: the integer part of t0geo is the barycentric (TDB) epoch
    of RA, DEC, f0, f1,and f2
f0    Pulsar rotation frequency (s**(-1))
f1    First derivative of pulsar frequency (s**(-2))
f2    Second derivative of pulsar frequency (s**(-3))
RMS    Root-mean-square radio timing residual, in milliperiods
O    Observer code
B    Blank for single pulsars, "*" for binaries.

The "Observer codes" have the following meanings and contact persons:
A    Australia Telescope National Facility (Dick Manchester)
B    Bologna (Nichi D'Amico)
C    Cornell/Carleton (Jim Cordes, Joel Weisberg)
J    Jodrell Bank (Andrew Lyne)
P    Princeton (Joe Taylor)

For binary pulsars, the orbital parameters will be put in a file
called "psrbin.dat".  A short dummy version of the file looks like this:

 PSR B           pb         a1*sin(i)      e          t0           omega 
   omdot    gamma     pbdot   O
-------------------------------------------------------------------------
-------------------------------
0655+64      88877.061781   4.1255679 0.00000000 46386.41030392 
180.000000 0.00000 0.000000   0.00D+00 P
1855+09    1065067.590793   9.2307810 0.00002157 47529.90001480 
276.493717 0.00000 0.000000   0.00D+00 P

where the parameter definitions are (for details, see Taylor and 
Weisberg,=20
Astrophys. J. 345, 434, 1989):

PSR    Pulsar name
pb    Orbital period (s)
a1sini    Projected semi-major axis (light seconds)
e    orbital eccentricity
t0    Barycentric time (TDB scale) of periastron (MJD)
omega    Longitude of periastron (degrees)
omdot    First derivative of omdot (degrees per Julian year)
gamma    Time-dilation and gravitational redshift parameter (s)
pbdot    First derivative of pb
O    Observer code

Many binary pulsars will not require the last three (relativistic)
parameters.

Finally, a word for programmers.  The Fortran statements used to
write out the file entries are listed below:

Main file:

      write(33,1120) psrname,irh,irm,rsec,decsgn,idd,idm,dsec,
     +      mjd1,mjd2,t0geo,f0,f1,f2,rmsmp,obsflag,binflag
1120      format(a8,2i3.2,f7.3,1x,a1,i2.2,i3.2,f6.2,2i6,f16.9,
     +      f18.13,1p,d13.5,d11.2,0p,f5.1,2(1x,a1))

Binary pulsar file:

      write(34,1100) psrname,pb,a1,e,t0mjd,omz,omdot,gamma,pbdot,
     +      obsflag
1100      format(a8,f17.6,f12.7,f11.8,f15.8,f11.6,f8.5,f9.6,
     +      1p,d11.2,0p,1x,a1)

Note that it appears that the BVC times are the mid-point of the corresponding superpacket. For that reason I modified the GCOMBvc class so that it simply returns the time instead of the start and stop times.

See also this good reference for the computations to do: https://asd.gsfc.nasa.gov/Craig.Markwardt/bary/.

• PULRGL - New task in which radio-phase zero is calculated at nominal epoch of pulsar ephemeris from Princeton catalogue. See COM-RP-ROL-DRG-65 for a description of the applied algorithm (from SCN No: ROL-55, Date: 1994-11-22)

To make progress I investigated how the INTEGRAL OSA is handling pulsars.

• Ephemeris is provided to the program through a FITS file. In the case of the Crab, one set of ephemeris parameters per month is required to follow the Crab phase over periods much longer than one month. The FITS table includes one ephemeris parameter set in each row with a validity time interval (specified through MJDSTART and MJDSTOP). The program can then automatically select a valid ephemeris parameter set (in a given row) for each processed science window. You just have to make sure that the validity intervals encompass the time span of all data included in the analyzed observation group.
• We are aware that using a FITS file may be inconvenient for most X-ray pulsars for which only one set of ephemeris parameters is needed.
• The meaning of the ephemeris parameters are standard, and the required information can easily be found trough the WWW.

• Starting with OSA 6.0, the spi_science_analysis script offers phased-resolved analysis capabilities.

In the spi_science_analysis.cpp there is the following code:

      if (run_phase_analysis == 1) {

      status = execute_spiros_phase(spi_science_analysis,
                    Group,
                    IRFDOL,
                    spiros_catalog,
                    coordinates,
                    nPhaseBins,
                    clobber,
                    chatter);

      if (run_phase_analysis == 1) {

      status = execute_phase_histogram(spi_science_analysis,
                       Group,
                       clobber,
                       chatter);

ISDCTask spi_phase_hist("spi_phase_hist");
 spi_phase_hist["inOG"]=Group;          
 spi_phase_hist["chatter"]=20;            
 spi_phase_hist["clobber"]=clobber;               
 spi_phase_hist["ephemDOL"]="bla";             
 spi_phase_hist["phaseBinNum"]=20;          
 spi_phase_hist["phaseSameWidthBin"]="yes";     
 spi_phase_hist["phaseBounds"]="bla";           
 spi_phase_hist["phaseSubtractOff"]="no";    
 spi_phase_hist["phaseOffNum"]=0;       
 spi_phase_hist["orbit"]="no";                 
 spi_phase_hist["asini"]=0.0;                 
 spi_phase_hist["porb"]=0.0;                   
 spi_phase_hist["T90epoch"]=0.0;              
 spi_phase_hist["ecc"]=0.0;                    
 spi_phase_hist["omega_d"]=0.0;                
 spi_phase_hist["pporb"]=0.0;             

    // load ephemeris data valid for the SWG tstart time
    Ephem * ephem = new Ephem();
    status = ephem->load(tpt->ephemDOL, tstart);
    if (status != ISDC_OK) {
      RILlogMessage(gRILFileRef, Error_2, 
            "Error %d in loading ephemeris file %s", status, tpt->ephemDOL );
      continue;
    } else {
      RILlogMessage(gRILFileRef, Log_2, 
            " Ephemeris loaded: ScW TSTART: %16.11f - ref time: %6.1f = %3.1f days",
            tstart, ephem->tref0, (tstart - ephem->tref0));
      if( tpt->chatter > 50 ) { ephem->dump(); }
    }

    // compute barycentric times for all events
    status = events->baryCalc(swg, *ephem, evtType_max);
    if (status != ISDC_OK) {
      RILlogMessage(gRILFileRef, Error_2, 
            "%d : Error converting evt IJD to barycentric time", 
            status);
      continue;
    }
              // compute orbital correction if applicable
              if ( tpt->orbit ) { 
                anEvtPtr->orbiTime(*tpt);
              } else {    
                anEvtPtr->setOrbi(anEvtPtr->bary());
              }

              // compute the event phase
              phase = anEvtPtr->phase(swg, *ephem);

/*----------------------------------------------------------------------------*/
/*            Evt::baryTime((dal_element * swg, Ephem const & ephem)          */
/* -------------------------------------------------------------------------- */
/* Return the barycentric time of the event                                   */
/*----------------------------------------------------------------------------*/
double Evt::baryTime(dal_element * swg, Ephem const & ephem)
{
  int status = ISDC_OK;
  _bary = _time;
  status = DAL3AUXconvertBarycentIJD(swg, 1, ephem.srcRA, ephem.srcDEC, 
                     DAL3AUX_TDB, DAL3AUX_ANY, TCOR_ANY, 
                     &_bary, status);
  if (status != ISDC_OK) {
    RILlogMessage(gRILFileRef, Error_2, 
          "%d : Error converting evt IJD to barycentric time", 
          status);
    return -1.0;
  }
  return _bary;
}

/*----------------------------------------------------------------------------*/
/*            Evt::orbiTime(Parameters const & tpt)                           */
/* -------------------------------------------------------------------------- */
/* Return the time of the event corrected for orbital motion                  */
/*----------------------------------------------------------------------------*/
double Evt::orbiTime(Parameters const & tpt)
{
  int status = ISDC_OK;

  double time = _bary;  

  double asini = tpt.asini;
  double porb = tpt.porb;
  double T90epoch = tpt.T90epoch;
  double ecc = tpt.ecc;
  double omega_d = tpt.omega_d;
  double pporb = tpt.pporb;

  if (porb==0){
    porb=1e-10;
  }
  if (ecc>=1){
    ecc=0.999999;
  }
  if (ecc<0){
    ecc=0;
  }
  double asini_d = asini/86400.0;
  double twopi = 2.0*3.1415927;
  double t = time;
  double omega=omega_d*3.1415927/180.0;
  double sinw=sin(omega);
  double cosw=cos(omega);
  for (int i=0;i<6;i++){
    double t1=t-T90epoch;
    double m=twopi*(
            t1/porb-0.5 *pporb*t1*t1/porb/porb)
      +3.1415927/2.0-omega;
    double eanom=m;
    for(int j=0;j<5;j++){
      eanom=eanom-(eanom-ecc*sin(eanom)-m)/(1.0-ecc*cos(eanom));
    }
    double sin_e=sin(eanom);
    double cos_e=cos(eanom);
    double z=asini_d*(sinw*(cos_e-ecc)+sqrt(1.0-ecc*ecc)*cosw*sin_e);
    double dz=twopi*asini_d/porb/(1.0-ecc*cos_e)*
      (sqrt(1.0-ecc*ecc)*cosw*cos_e-sinw*sin_e);
    double f=t+z-time;
    double df=1.0+dz;
    if (df==0){
      df=1e-10;
    }
    t=t-f/df;
  }
  _orbi = t;
  return _orbi;
}

/*----------------------------------------------------------------------------*/
/*                            Evt::phase(Ephem const ephem)                   */
/* -------------------------------------------------------------------------- */
/* Return the event detector number (pseudo det number if multiple)           */
/*----------------------------------------------------------------------------*/
double Evt::phase(dal_element * swg, Ephem const & ephem)
{
  // deltaT = time (in seconds) betwwen phase zero next to reference time and 
  //         event time corrected for barycentric and for orbital motions if any
  double deltaT = 86400.0 * ( _orbi - ephem.tref0 ) - ephem.t2peak0;
  // Number of turns during delatT
  double Turns = ephem.freq0 * deltaT + ephem.fdot0 * deltaT * deltaT/2.0 + 
              ephem.fdotdot0 * deltaT * deltaT * deltaT/6.0;
  _phase = (Turns - floor(Turns));
  return _phase;
}

I checked whether there may be some issue with the precision of GTime at typical MJD times of CGRO:

# Check GTime precision at the beginning of CGRO observtion (start of VP0001)
>>> print(gammalib.GTime('MJD 48392.10000000')-gammalib.GTime('MJD 48392.00000000'))
8639.99999988 <=> 8640 (-1.19209289551e-07)
>>> print(gammalib.GTime('MJD 48392.01000000')-gammalib.GTime('MJD 48392.00000000'))
864.000000119 <=> 864 (1.19209289551e-07)
>>> print(gammalib.GTime('MJD 48392.00100000')-gammalib.GTime('MJD 48392.00000000'))
86.3999996185 <=> 86.4 (-3.81469732247e-07)
>>> print(gammalib.GTime('MJD 48392.00010000')-gammalib.GTime('MJD 48392.00000000'))
8.63999974728 <=> 8.64 (-2.52723694416e-07)
>>> print(gammalib.GTime('MJD 48392.00001000')-gammalib.GTime('MJD 48392.00000000'))
0.864000201225 <=> 0.864 (2.01225280771e-07)
>>> print(gammalib.GTime('MJD 48392.00000100')-gammalib.GTime('MJD 48392.00000000'))
0.0864000320435 <=> 0.0864 (3.20434570267e-08)
>>> print(gammalib.GTime('MJD 48392.00000010')-gammalib.GTime('MJD 48392.00000000'))
0.00864005088806 <=> 0.00864 (5.08880615233e-08)
>>> print(gammalib.GTime('MJD 48392.00000001')-gammalib.GTime('MJD 48392.00000000'))
0.000863671302795 <=> 0.000864 (-3.2869720459e-07)
>>> print(gammalib.GTime('MJD 48392.000000001')-gammalib.GTime('MJD 48392.000000000'))
8.60691070557e-05 <=> 0.0000864 (-3.30892944336e-07)

# Check GTime precision for one tic at the beginning of CGRO observations
>>> print(gammalib.com_time(8392,1)-gammalib.com_time(8392,0))
0.000125050544739
>>> print(gammalib.com_time(8392,1)-gammalib.com_time(8392,0)-0.000125)
5.05447387695e-08

# Check GTime precision for one tic at the end of CGRO observations
>>> print(gammalib.com_time(11690,1)-gammalib.com_time(11690,0))
0.000124990940094
>>> print(gammalib.com_time(11690,1)-gammalib.com_time(11690,0)-0.000125)
-9.05990600586e-09

I now compared the Barycentre correction in barycor.py to the one done in GCOMBvcs::tdelta. For line 53 in the psrtime file, which corresponds to VP0001, the

The times for which the evaluation was done is not exactly the same, which may explain why travel is not identical. Yet the agreement seems to be reasonable for about one day of difference. Note also the CGRO is not at the geocentre.

Maybe this code from the ftools/axBary tasks that is found in bary.c has some relevance:

timeparams

// Read line from PSRTIME file
timeline (line, bsource, &ra, &dec, &mjd1, &mjd2, &mjdi, &mjdf, &fd, &dfd, &ddfd, &dddfd, &rms, &fbin, &tdenum, jsource) ;

// Stuff
ptp->f = fd ;                     // f0
ptp->df = dfd ;                   // f1
ptp->ddf = ddfd ;                 // f2
ptp->dddf = dddfd ;               // f3
i = mjdf * 100.0 + 0.5 ;          // Keep the first 3 digits of the fractional part,
mjdftrun = (double) i / 100.0 ;   // ... drop the others (copes with 48499.999999517 style dates)
ptp->t0.MJDint = mjdi ;           // MJD integer part
ptp->t0.MJDfr = mjdftrun ;        // MJD fractional part (should be 0 or 1)
ptp->t0.ts = TDB ;                // Time system (SSB)

// Convert UTC to TT
mtime.MJDint = mjdi ;
mtime.MJDfr = mjdf + roffset / SECDAY ; // roffset: extra offset to be applied to radio zero phase (in s)
mtime.ts = UTC ;
toTT (&mtime) ;  // Converts UTC->TT (i.e. adds leap seconds + 32.184), mtime.ts = TT on output

tr           = barycorr (&mtime, dir, scposn0) ;
dtt          = tr + ( (mtime.MJDint - ptp->t0.MJDint) + (mtime.MJDfr - ptp->t0.MJDfr) ) * SECDAY ;
ptp->radioph = -fd*dtt - dfd*dtt*dtt*0.5 - ddfd*dtt*dtt*dtt/6.0 - dddfd*dtt*dtt*dtt*dtt/24.0 ;
// TDB-TT      = tr
// DT (=Tr-T0) = dtt
// Phase(T0)   = ptp->radioph

The question to answer now is in which time system the COMPTEL times are provided.

The evpbin02.pulbvc.f code takes on input UTC time of the event, which may suggest that COMPTEL times are in fact UTC times.

• it computes the centre time of an observation using the INTSTA and INTEND information in Julian Days and Modified Julian Days and hands this information over to PULA32 to extrapolate the pulsar frequency, derivative and zero phase to the center of the epoch
• it seems that in VP0032, which corresponds to 1992-06-25T01:00:00, a correction of 2.042144 s in the onboard time was applied. There exists in fact a puliex02.timcor.F routine that states GRO Clock too fast by 2.042144 seconds (=16337.152 Tags).
• the pulbvc.f function is called with FTIME, which is the event time read from the EVP dataset. This implies that COMPTEL event times are in UTC.

The following time correction to the EVP data is employed:

      IF ( ( FFTIME(1).LT.8798 ).OR.
     1     ( (FFTIME(1).EQ.8798).AND.(FFTIME(2).LT.28800000)))THEN
         FTIME(1) = FFTIME(1)
         FTIME(2) = FFTIME(2) - ( 2.042144D0 * 8000.D0 )
         IF ( FTIME(2) . LT. 0 )THEN
            FTIME(1) = FTIME(1) - 1
            FTIME(2) = FTIME(2) + ZZTIC
         END IF
      ELSE
         FTIME(1) = FFTIME(1)
         FTIME(2) = FFTIME(2)
      END IF

This time correction needs to be implemented in the gammalib::com_time(), gammalib::com_tjd() and gammalib::com_tics() methods.

Note that I discovered that the methods did already consider that COMPTEL event times are given in UTC!

Here is an selection of PPT datasets that were used at ROL for pulsar analysis in the COMPASS system. If I could get these datasets I may check the validity of the implementation:

R,PPT,98,Crab ephemeris : Obs. 213  from Princeton Catalogue,9047,0,9144,0,KUIPER,X,10/28/93 14:37,,,1,PULCCP,98,,,,,,,
R,PPT,100,Crab ephemeris : Obs. 221  from Princeton Catalogue,9083,0,9173,0,KUIPER,X,10/29/93 16:16,,,1,PULCCP,100,,,,,,,
R,PPT,136,Crab Princeton ephemeris 9296-9379 mid 9337,9296,0,9379,0,KUIPER,X,6/8/94 12:58,,,1,PULCP1,130,,,1,100,,11/27/95 10:13,
R,PPT,138,Crab Princeton ephemeris 9340-9411 mid 9375 valid for Obs. 321.1+321.5,9340,0,9411,0,KUIPER,X,6/20/94 15:02,,,1,PULCP1,132,,,1,100,,11/27/95 10:13,
R,PPT,141,Crab Princeton ephemeris 8331-8412 mid 8371 valid for Obs. 1,8331,0,8412,0,KUIPER,X,8/1/94 10:54,,,1,PULCP1,135,,,1,,,11/27/95 10:13,
R,PPT,146,"Crab (PSR B0531+21) 9047-9144 mid 9095 valid for Obs. 213 , 221",9047,0,9144,0,KUIPER,X,8/9/94 10:45,,,1,PULCP1,140,,,1,100,,11/27/95 10:13,
R,PPT,155,"Crab (PSR B0531+21) 8331-8412 mid 8371 valid for Obs. 0.3,0.4,0.5,1.0",8331,0,8412,0,KUIPER,X,9/8/94 12:36,,,1,PULCP1,149,,,1,100,,11/27/95 10:13,
R,PPT,157,Crab Princeton ephemeris 8406-8493 mid 8449 valid for Obs. 2.5,8406,0,8493,0,KUIPER,X,9/20/94 13:16,,,1,PULCP1,151,,,1,100,,11/27/95 10:13,
R,PPT,158,Crab Princeton ephemeris 8499-8637 mid 8568 valid for Obs. 15,8499,0,8637,0,KUIPER,X,9/20/94 13:18,,,1,PULCP1,152,,,1,100,,11/27/95 10:13,
R,PPT,159,Crab Princeton ephemeris 8637-8801 mid 8719 valid for Obs. 31,8637,0,8801,0,KUIPER,X,9/20/94 13:23,,,1,PULCP1,153,,,1,100,,11/27/95 10:13,
R,PPT,160,"Crab Princeton ephemeris 8801-8913 mid 8857 valid for Obs. 36.0,36.5,39,41",8801,0,8913,0,KUIPER,X,9/20/94 13:26,,,1,PULCP1,154,,,1,100,,11/27/95 10:13,
R,PPT,161,"Crab Princeton ephemeris 8884-8943 mid 8913 valid for Obs. 41,44",8884,0,8943,0,KUIPER,X,9/20/94 13:30,,,1,PULCP1,155,,,1,100,,11/27/95 10:13,
R,PPT,164,Crab ephemeris for VP 337 with MJD 49598,9526,0,9671,0,KUIPER,X,11/23/94 16:37,,,1,PULCP1,158,,,1,100,,11/27/95 10:13,
R,PPT,193,Crab VP 337 49598,9526,0,9671,0,KUIPER,X,2/7/95 16:58,,,1,PULCP2,28,1,258662,1,200,,11/27/95 10:13,
R,PPT,204,Crab VP 0.2+0.3+0.4+0.5+1 Center 8371,8331,0,8412,0,KUIPER,X,5/23/95 14:42,,,1,PULCP2,39,1,274360,1,200,,11/27/95 10:13,
R,PPT,219,PSR 0531+21 (Crab) VP 412.0+413.0 Center 9780,9751,0,9810,0,KUIPER,X,8/3/95 16:11,,,1,PULCP2,54,1,285935,1,200,,11/27/95 10:13,
R,PPT,230,PSR 0531+21 (Crab) VP 419.1+419.5+420 Center 9842,9791,0,9893,0,KUIPER,X,11/1/95 15:17,,,1,PULCP2,65,1,302171,1,200,,11/27/95 10:13,R,PPT,235,PSR 0531+21 (Crab) VP 426.0+427       Center 9961,9931,0,9992,0,KUIPER,X,11/2/95 15:39,,,1,PULCP2,70,1,302415,1,200,,11/27/95 10:13,
R,PPT,238,PSR 0531+21 (Crab) VP 2.5 Center 8449,8406,0,8493,0,KUIPER,X,11/21/95 11:18,,,1,PULCP2,73,1,304244,1,200,,11/27/95 10:13,
R,PPT,239,PSR 0531+21 (Crab) VP 15  Center 8568,8499,0,8637,0,KUIPER,X,11/21/95 11:29,,,1,PULCP2,74,1,304247,1,200,,11/27/95 10:13,
R,PPT,240,PSR 0531+21 (Crab) VP 31  Center 8719,8637,0,8801,0,KUIPER,X,11/21/95 11:33,,,1,PULCP2,75,1,304250,1,200,,11/27/95 10:13,
R,PPT,241,PSR 0531+21 (Crab) VP 36.0+36.5+39+41 Center 8857,8801,0,8913,0,KUIPER,X,11/21/95 11:39,,,1,PULCP2,76,1,304253,1,200,,11/27/95 10:13,
R,PPT,242,PSR 0531+21 (Crab) VP 41+44           Center 8913,8884,0,8943,0,KUIPER,X,11/21/95 11:46,,,1,PULCP2,77,1,304256,1,200,,11/27/95 10:13,
R,PPT,243,PSR 0531+21 (Crab) VP 310             Center 9337,9296,0,9379,0,KUIPER,X,11/21/95 12:06,,,1,PULCP2,78,1,304259,1,200,,11/27/95 10:13,
R,PPT,244,PSR 0531+21 (Crab) VP 321.1+321.5     Center 9375,9340,0,9411,0,KUIPER,X,11/21/95 12:09,,,1,PULCP2,79,1,304262,1,200,,11/27/95 10:13,
R,PPT,248,PSR 0531+21 (Crab) VP 502             Center 10009,9991,0,10028,0,KUIPER,X,1/10/96 11:51,,,1,PULCP2,83,1,307331,1,,,,
R,PPT,273,PSR 0531+21 (Crab) VP 502             Center 10009,9991,0,10028,0,KUIPER,X,3/28/96 11:22,,,1,PULCP2,108,1,317884,1,,,,
R,PPT,274,PSR 0531+21 (Crab) VP 502             Center 10009,9991,0,10028,0,KUIPER,X,3/28/96 11:25,,,1,PULCP2,109,1,317890,1,,,,
R,PPT,292,PSR B0531+21 (Crab) 10204 10235  VP 520 (Correction for Feb 1996),10204,0,10235,0,KUIPER,X,7/18/96 15:54,,,1,PULCP3,17,2,327558,1,,,,
R,PPT,294,PSR B0531+21 (Crab) 10084 10098  VP 510.0/510.5,10084,0,10155,0,KUIPER,X,7/19/96 15:59,,,1,PULCP3,19,2,327696,1,,,,
R,PPT,296,PSR B0531+21 (Crab) 10204 10235  VP 520,10204,0,10235,0,KUIPER,X,8/2/96 10:02,,,1,PULCP3,21,2,328248,1,,,,
R,PPT,300,PSR B0531+21 (Crab) 10142 10259  VP 520,10142,0,10259,0,KUIPER,X,11/4/96 14:57,,,1,PULCP3,25,2,334202,1,,,,
R,PPT,302,PSR B0531+21 (Crab) 10260 10292  VP 523 Crab DEC2000 Modified,10260,0,10292,0,KUIPER,X,11/5/96 17:41,,,1,PULCP3,27,2,334268,1,,,,
R,PPT,305,PSR B0531+21 (Crab) 10267-10310  VP 526 Crab DEC2000 Modified,10267,0,10310,0,KUIPER,X,11/25/96 15:32,,,1,PULCP3,30,2,335891,1,,,,
R,PPT,307,PSR B0531+21 (Crab) 10267-10310  VP 526 Crab new entry,10267,0,10310,0,KUIPER,X,12/9/96 11:31,,,1,PULCP3,32,2,337438,1,,,,
R,PPT,309,PSR B0531+21 (Crab) 10260 10292  VP 523 Crab,10260,0,10292,0,KUIPER,X,12/9/96 11:50,,,1,PULCP3,34,2,337450,1,,,,
R,PPT,311,PSR B0531+21 (Crab) 10297-10325  VP 527/528 (VP 526 Partly) New entry,10297,0,10325,0,KUIPER,X,12/10/96 14:45,,,1,PULCP3,36,2,337792,1,,,,
R,PPT,330,PSR B0531+21 (Crab) 10491-10539  VP 616.1,10491,0,10539,0,KUIPER,X,7/25/97 15:01,,,1,PULCP3,55,2,355079,1,,,,
R,PPT,365,PSR B0531+21  9400- 9551 / mid 9475  Princeton  VP 325,9400,0,9551,0,KUIPER,X,1/5/99 19:47,,,1,PULCP3,90,2,408590,1,,,,
R,PPT,367,PSR B0531+21 10311-10409 / mid 10360 Princeton  VP 520.9,10311,0,10409,0,KUIPER,X,1/6/99 11:59,,,1,PULCP3,92,2,408695,1,,,,
R,PPT,369,"PSR B0531+21 10530-10575 / mid 10552 Princeton  VP (617.1),617.3,617.5,617.7",10530,0,10575,0,KUIPER,X,1/6/99 18:26,,,1,PULCP3,94,2,408743,1,,,,
R,PPT,371,PSR B0531+21 10630-10660 / mid 10645 Princeton  VP 623 Partly,10630,0,10660,0,KUIPER,X,1/8/99 15:12,,,1,PULCP3,96,2,408869,1,,,,
R,PPT,372,PSR B0531+21 10655-10685 / mid 10670 Princeton  VP 623 Partly,10655,0,10685,0,KUIPER,X,1/8/99 15:14,,,1,PULCP3,97,2,408872,1,,,,
R,PPT,388,PSR B0531+21 10957-11057 / mid 11007 Princeton  VP 724.5,10957,0,11057,0,KUIPER,X,1/18/99 13:50,,,1,,,,,1,100,,1/18/99 15:31,
R,PPT,401,PSR B0531+21 11299-11330 / mid 11314 Jodrell Bank VP 815/816,11299,0,11330,0,KUIPER,X,3/21/00 16:16,,,1,PULCP3,126,2,437131,1,,,,
R,PPT,402,PSR B0531+21 11422-11451 / mid 11436 Jodrell Bank VP 829,11422,0,11451,0,KUIPER,X,3/21/00 16:16,,,1,PULCP3,127,2,437134,1,,,,
R,PPT,403,PSR B0531+21 11513-11543 / mid 11528 Jodrell Bank VP 903.1,11513,0,11543,0,KUIPER,X,3/21/00 16:17,,,1,PULCP3,128,2,437137,1,,,,
R,PPT,407,PSR B0531+21 11635-11665 / mid 11650 Jodrell Bank VP 918.5,11635,0,11665,0,KUIPER,X,11/6/00 17:37,,,1,PULCP3,132,2,439030,1,,,,
R,PPT,408,PSR B0531+21 11665-11696 / mid 11680 Jodrell Bank VP 918.5,11665,0,11696,0,KUIPER,X,11/6/00 17:37,,,1,PULCP3,133,2,439033,1,,,,
R,PPT,426,PSR B0531+21 11575-11604 / mid 11590 Jodrell Bank February 2000  Chandra Crab,11575,0,11604,0,KUIPER,X,3/28/01 13:31,,,1,PULCP3,151,2,443930,1,,,,

Note that the phase shift cannot be explained by numerical differences. If identical pulsar positions are used, the pulphi.py software produces identical radio phase shifts. Here the result for the Jodrell bank ephemerides

 Integer part Tgeo (MJD) ...: 48299
 Fraction part Tgeo (MJD) ..: 1.340013e-07
 Correction to Tgeo (ms) ...: 0.000000
 Correction to Tgeo (days) .: 0.000000
 Calculated JD (int) .......: 2448300
 Calculated JD (frac) ......: -0.500000
 Number of Leap seconds ....: 16
 Frequency .................: 29.951601
 Frequency derivative ......: -3.777260e-10
 Second freq. derivative ...: 0.000000e+00
 Pos1 ......................: 0.102808
 Pos2 ......................: 0.921371
 Pos3 ......................: 0.374841
 Pulsar position RA 2000 ...: 83.633217
 Pulsar position DEC2000 ...: 22.014464
 UTC--> TDT ................: 58.184000 sec
 TDT--> TDB ................: 0.001047 sec
 UTC--> TDB ................: 58.185047 sec
 Traveltime difference .....: 255.737957 sec
 Relativistic delay ........: -0.000004 sec
 Total delay ...............: 313.923009 sec
 Time diff. to nom. epoch ..: 313.934586 sec
 Radio-phase 0 at nom.epoch : 0.843473
=== GPulsarEphemeris ===
 Pulsar name ...............: PSR B0531+21
 Pulsar Right Ascension ....: 83.633225 deg
 Pulsar Declination ........: 22.0144583333333 deg
 Time system of ephemeris ..: UTC
 Reference epoch ...........: MJD 48299
 Phase .....................: -0.845097897454252
 Frequency .................: 29.9516010684378 Hz
 Frequency derivative ......: -3.77726e-10 Hz/s
 2nd frequency derivative ..: 0 Hz/s^2
 Period ....................: 33.3871968218011 ms
 Validity MJD range ........: 48281.9993265741 - 48315.9993265741
 Validity UTC range ........: 1991-01-25T23:59:02 - 1991-02-28T23:59:02
 pulphi radio phase ........: -0.210862
 GammaLib radio phase ......: -0.210863
 Difference ................: 0.000001

 Integer part Tgeo (MJD) ...: 48299
 Fraction part Tgeo (MJD) ..: 1.340013e-07
 Correction to Tgeo (ms) ...: 0.000000
 Correction to Tgeo (days) .: 0.000000
 Calculated JD (int) .......: 2448300
 Calculated JD (frac) ......: -0.500000
 Number of Leap seconds ....: 16
 Frequency .................: 29.951601
 Frequency derivative ......: -3.777260e-10
 Second freq. derivative ...: 0.000000e+00
 Pos1 ......................: 0.102807
 Pos2 ......................: 0.921371
 Pos3 ......................: 0.374841
 Pulsar position RA 2000 ...: 83.633225
 Pulsar position DEC2000 ...: 22.014458
 UTC--> TDT ................: 58.184000 sec
 TDT--> TDB ................: 0.001047 sec
 UTC--> TDB ................: 58.185047 sec
 Traveltime difference .....: 255.738011 sec
 Relativistic delay ........: -0.000004 sec
 Total delay ...............: 313.923063 sec
 Time diff. to nom. epoch ..: 313.934641 sec
 Radio-phase 0 at nom.epoch : 0.845098
=== GPulsarEphemeris ===
 Pulsar name ...............: PSR B0531+21
 Pulsar Right Ascension ....: 83.633225 deg
 Pulsar Declination ........: 22.0144583333333 deg
 Time system of ephemeris ..: UTC
 Reference epoch ...........: MJD 48299
 Phase .....................: -0.845097897454252
 Frequency .................: 29.9516010684378 Hz
 Frequency derivative ......: -3.77726e-10 Hz/s
 2nd frequency derivative ..: 0 Hz/s^2
 Period ....................: 33.3871968218011 ms
 Validity MJD range ........: 48281.9993265741 - 48315.9993265741
 Validity UTC range ........: 1991-01-25T23:59:02 - 1991-02-28T23:59:02
 pulphi radio phase ........: -0.845098
 GammaLib radio phase ......: -0.845098
 Difference ................: -0.000000

The pulse shift can however be explained by differences in the assumed position of the Crab pulsar. This is illustrated in the table below where small shifts (in units of degrees) were added to the pulsar position in the pulphi.py computation, while no shifts were applied in the GammaLib computation. A position inaccuracy of 0.54 arcsec in Right Ascension leads to a phase shift of -0.03! In comparison, the difference in Right Ascension between the SIMBAD ICRS position and the position provided in the Princeton file is 0.5355.

Note that according to http://simbad.u-strasbg.fr/simbad/sim-ref?bibcode=2011A%26A...533A..10L the Crab pulsar proper motion in Right Ascension is -13.0±0.2mas/yr and in Declination is +2.9±0.1mas/yr. Since the SIMBAD position is for the epoch 2000 while the ephemerides are for 1991, this translated into a motion of 0.117 arcsec in Right Ascension and 0.0119 arcsec in Declination. Hence the proper motion alone cannot explain the differences in the pulsar position.