BARY
$Id: 00README,v 3.7 2001/09/10 19:20:21 arots Exp arots $

Arnold Rots
USRA, CfA
25 November 1997
Updated 4 May 1998
Updated 2 July 1998
Updated 25 August 1998
Updated 3 September 1998
Updated 15 January 1999
Updated 1 June 1999
Updated 1 November 1999
Updated 31 January 2000
Updated 7 December 2000
Updated 10 September 2001
Updated 6 October 2004
Updated 27 October 2004
Updated 8 June 2006

Update Notes
2 July 1998
Added simple time conversions and FITS formats.
25 August 1998
Version 3.0
3 September 1998
Added axBary.
Consolidated clock corrections and orbit positions
in separate source files (clock.c, scorbit.c).
15 January 1999
Corrected lunar ephemeris calculation (thanks to Jeremy Bailey)
1 June 1999
Changed toTT and toUTC to return change in seconds;
barycorr does not change the MJDTime struct anymore.
1 November 1999
Put checks in for ASC_DATA and/or TIMING_DIR being defined.
Allow for TIMEDEL and DELTAT in scorbit.
31 January 2000
More gracious handling of faulty GTI tables in Chandra products
7 December 2000
Chandra is supported.
Various bugfixes.
10 September 2001
bary.c: make calls to dpleph with TDB time, rather than TT
axBary.c: make sure output file is writable; add TIMZERO to TSTART/TSTOP
6 October 2004
bary.c: correct handling of binary pulsars zero-phase.
fBadd.c: Proper handling of offsets.
fBssum.c: Increased decimals in print statements.
axBary.c: Increased precision of RA, Dec, TimeZero.
baryFase.c: Added Chandra mission.
clock.c: Removed HRC instrumental delay (corrected for in other software).
         Bill Davis value for Chandra clock correction.
27 October 2004
faseBin.c: Errors for zero counts set to zero, in order to prevent
         errors for sums to skyrocket due to empty bin.
8 June 2006
clock.c: Chandra clock correction is taken care of elsewhere; make
         correction contingent of value of TIERABSO

This package is Y2k-compliant.

Note:
Until the Chandra GTI tables are fixed, axBary will still produce a
warning for HDU 2 in Chandra files.  The minimal fix is to add
the keyword MJDREF=50814 to the second extension's header.

A package containing:
  - the JPL DE200 and DE405 ephemerides (full Chebyshev polynomials)
    in FITS format;
  - ANSI standard C code to access those ephemerides;
  - new barycenter correction code that can use either ephemeris and
    starts from any of these time systems: UTC, ET, TT, TCG, TAI;
  - a phase binner that bins photon events in absolute phase, using
    the Princeton radio pulsar database (format), but including
    quartic and quitic term support;
  - a utility, baryFase, that calculates barycenter corrections and
    (optionally) absolute phase for any time and celestial position.
  - a utility to FITSify JPL ephemeris files;
  - Some simple time conversion functions (see recipe section at end);
  - a Makefile.

One should be aware that the JPL ephemeris, as distributed in "Unix"
format, are really IEEE; hence, those files will work just fine on,
say, a Sparc station, but will cause trouble on Alphas and Intel
machines.  The FITS format circumvents this problem.

Although the application is for barycenter corrections, the dpleph
version distributed as part of this package will allow retrieval
of positions for all members of the solar system contained in the
JPL DE ephemerides.

The barycenter package is geared toward space missions and detects
the identity of the mission through the TELESCOP or MISSION keywords
in the FITS headers.  However, it will be easy to include ground-based
observatories by inserting an appropriate geographic position function
for the "orbit ephemeris".  Indeed, the currently supported "missions"
are XTE, AXAF, and the Geocenter.  Missions may be added by
incorporating support in the files bary.h, clock.c, scorbit.c, and,
optionally, phaseHist.c.


This directory contains the following files:

00README       This file
JPLEPH.200     FITS version of the JPL DE200 lunar and planetary
               ephemeris for 1950-2050
JPLEPH.405     FITS version of the JPL DE405 lunar and planetary
               ephemeris for 1950-2050
JPLEPH         Symbolic link to the default ephemeris
tai-utc.dat    Leap second file
psrtime.dat    Pulsar timing ephemeris
psrbin.dat     Binary pulsar orbit ephemeris
tdc.dat        XTE clock correction file
dpleph.c       Functions to access the JPL ephemeris
bary.c         General utilities and functions to calculate barycenter
               corrections
ctatv.c        Function to calculate TDB - TT
clock.c        Function package for mission-specific clock corrections
scorbit.c      Function package for mission-specific orbit positions
axBary.c       Function to effect barycenter corrections in FITS files
axBarym.c      Main to execute axBary from the command line
baryFase.c     Utility to calculate individual barycenter corrections
               and phase
phaseHist.c    Phase binner function: bins events in energy and
               absolute phase according to timing ephemeris
faseBin.c      The phase binner main function, or wrapper, for
               phaseHist
fBREADME       faseBin documentation
bary.h         Header file required by all
xCC.c          XTE clock correction function
bineph2fits.c  Utility to convert ephemerides to FITS

Note that the *.dat files periodically need updating.


Do not ever use DE405 on a FITS file that has:
    TIMESYS='TDB'
  but not:
    RADECSYS='ICRS' or PLEPHEM='JPL-DE405'
TIMESYS='TDB' _with_    RADECSYS='ICRS' and/or PLEPHEM='JPL-DE405'
                           should use denum=405.
TIMESYS='TDB' _without_ RADECSYS='ICRS' or     PLEPHEM='JPL-DE405'
                           should use denum=200.
DE405 may be used in conjunction with FK5 spacecraft orbit ephemeris.
The maximum error is 2 ns for each earth radius that the spacecraft is
removed from the geocenter.
All positions provided by Tempo are on the DE200's reference frame,
however - closer to FK5.  This could introduce errors of up to 0.02 ms.
A function, c200to405, is provided to convert DE200 positions to ICRS.
Even so, it is recommended that DE200 is used for Tempo solutions that
are based on DE200; apparently, efforts are underway to make Tempo
support DE405.
Bary supports a relatively free-format version of psrtime.dat with
quartic terms.


The following data files are expected in either $TIMING_DIR or
$LHEA_DATA (if in both, $TIMING_DIR will take precedence); the URLs
refer to the default files.  You may want to customize LHEA_DIR for
your environment; $LHEA_DATA is used at NASA/GSFC/LHEA (HEASARC),
$ASC_DATA at CfA/SAO (ASC).

JPLEPH, JPLEPH.200, JPLEPH.405  See above and:
         ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/JPLEPH.200
         ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/JPLEPH.405
psrtime.dat    Pulsar timing ephemeris (in the format of jcat from the
               Princeton radio pulsar database):
         ftp://pulsar.princeton.edu/gro/jcat
psrbin.dat     Binary pulsar orbit ephemeris file (in the format of
               psrbin.dat from the Princeton radio pulsar database):
         ftp://pulsar.princeton.edu/gro/psrbin.dat
tdc.dat        XTE clock corrections, from XTE ftp archive:
         ftp://heasarc.gsfc.nasa.gov/xte/calib_data/clock/tdc.dat
tai-utc.dat    leap second file, from USNO:
         ftp://maia.usno.navy.mil/ser7/tai-utc.dat


The following library is required, with its header files:

cfitsio        From HEASARC; fitsio.h, fitsio2.h, longnam.h
               ftp://heasarc.gsfc.nasa.gov/software/fits/c

Note that the paths to this library and its include files needs to be
set in the Makefile.  The Makefile allows building the libbary.a
library and the executables.


Time is kept in three possible ways:
  The basic convention is:
    MJDTime t ;  The MJDTime struct is defined below;
                   a C++ class would be better.
                       typedef struct {
                           enum TimeSys ts ;
                           long MJDint ;
                           double MJDfr ;
                       } MJDTime ;
  The JPL ephemeris functions use:
    double time[2] ;  Where:
      time[0]   Integer part of JD
      time[1]   Fractional part of JD; -0.5 <= time[1] < 0.5
  MET uses a single double:
    double time ;  Where:
      time      Seconds or days since MJDREF, depending on the value of
                TIMEUNIT


baryFase
Usage:    baryFase [-d[ebug]]
The dialog asks the following questions:

Mission/observatory
  At this point only XTE, AXAF, and Geocenter are supported, but
  others could easily be added.
Orbit ephemeris file
  Only for missions/observatories that require such a file; if none is
  provided, baryFase defaults to Geocenter.
JPL ephemeris number
  200, 400, or default.
Input time scale
  Any legitimate time scale; default: TT.
Time unit
  "s" (second), "d" (day), or "m" (MJD day).
Reference MJD
  In two parts: integer and fractional, for precision reasons.
  Only if time unit is s or d.
  Note: time (below) can only be given as one double precision
  variable; if one uses MJD and precision is an issue, one may still
  want to specify time unit d and provide an MJD offset.
Source name
  For looking up timing ephemeris in pulsar daabase.
Binary orbit correction?
  "y" for binary pulsars.
Timing parameters
  If no source was specified, the timing ephemeris ay be given here.
  If only RA and Dec are given, no phase will be calculated.
Time
  Time for which the barycenter correction and absolute phase are to
  be calculated.  A simple <CR> or negative time will terminate the
  program.


axBary
Usage:    axBary -i orbitFile -f inputDataFile -o outFile [-ra RA]
                [-dec DEC] [-ref refFrame]

axBary is a multi-mission function that applies barycenter corrections
to HFWG-compliant FITS files.  It currently recognizes missions "XTE",
"AXAF", and "Geocenter" and has been provided with a simple "main"
wrapper that allows running axBary from the command line.
It copies inputDataFile into outFile, applying barycenter corrections
to all times, using orbit ephemeris file orbitFile, ra, dec, and
refFrame "FK5" (DE200) or "ICRS" (default; DE405).  When RA and Dec
are not provided, an effort is made to retrieve the information from
the FITS file, though this effort could be made more sophisticated,
taking into account other missions and telescopes.  J2000 positions
are required.

The core function is:
  int axbary (char *inFile, char *orbitFile, char *outfile,
              double ra, double dec, char *refFrame, int debug)
    Arguments:
      inFile    char*   Path of file to be corrected
      orbitFile char*   Path of orbit ephemeris file
      outFile   char*   Path of output file to be created
      ra        double  RA to be used for barycenter corrections
      dec       double  Dec to be used for barycenter corrections
                        If either ra or dec is out of bounds (0 -> 360,
                        -90 -> +90, respectively), RA_PNT and DEC_PNT
                        are used; if either is not available, RA_NOM
                        and DEC_NOM are used; if either of those is not
                        available, an error condition results
                        J2000 is required
      refFrame  char*   Reference frame to be used (FK5 or ICRS)
      debug     int     Debug message switch


Externally referenced functions:

dpleph.c

  int initephem (int ephnum, int *denum, double *c, double *radsol,
                 double *msol)
     Initialize ephemeris access.
     Arguments:
       Input
      ephnum   requested DE number of ephemeris (0: default)
       Output:
       denum   DE number of ephemeris sed
           c   speed of light (m/s)
      radsol   solar radius (light secs)
        msol   GM(solar), using light second as unit of length


  const double *dpleph (double *jd, int ntarg, int ncent)
      Return lunar/planetary ephemeris positions
      double[2]    jd   JD time for which ephemeris is requested
      int       ntarg   Target for which position is requested
      int       ncent   Center for which position is requested
    dpleph returns a pointer to an array containing the concatenated state
    vectors (X, Y, Z, VX, VY, VZ) of the target and the center with respect
    to the barycenter of the solar system at time jd[0]+jd[1] (in JD(TT)).

    The numbering convention for ntarg and ncent is:
      1 = Mercury            8 = Neptune
      2 = Venus              9 = Pluto
      3 = Earth             10 = Moon
      4 = Mars              11 = Sun
      5 = Jupiter           12 = Solar system barycenter
      6 = Saturn            13 = Earth-Moon barycenter
      7 = Uranus            14 = Nutations (longitude and obliquity; untested)
                            15 = Librations 
    This numbering scheme is 1-relative, to be consistent with the
    Fortran version.
  
    posn[12]   Interpolated quantities requested:
                 posn[0:2] ntarg position    posn[3:5]  ntarg velocity
                 posn[6:8] ncent position    posn[9:11] ncent velocity
                 units are lightseconds and lightseconds/s
  
                 posn[0:1] ntarg long/lat    posn[2:3]  ntarg rates
                 posn[0:2] ntarg Euler angs  posn[3:5]  ntarg rates
                 units are radians and radians/s (librations)


bary.c

  int baryinit (enum Observatory obs, char *fromts, char *tots,
                double mjdi, double mjdf, char *timeunit, int ephnum)
    Initialize barycenter correction environment; calls initephem.
    obs:      Observatory (Unknown, Geocenter, XTE, AXAF)
    fromts:   Time system of input
    tots:     Time system to convert to
    mjdi:     MJDREFI; use mission default if =-1.0
    mjdf:     MJDREFF
    timeunit: Unit of time (s or d)
    ephnum:   Ephemeris number requested (200, 405; 0: default)
    return:   denum; 0 upon error

  int timeinit (enum Observatory obs, char *fromts, char *tots,
                double mjdi, double mjdf, char *timeunit)
    Initialize time conversion environment.
    obs:      Observatory (Unknown, Geocenter, XTE, AXAF)
    fromts:   Time system of input
    tots:     Time system to convert to
    mjdi:     MJDREFI; use mission default if =-1.0
    mjdf:     MJDREFF
    timeunit: Unit of time (s or d)
    return:   0 upon error (for compatibility with baryinit)

  double barycorr (MJDTime *mjdTT, double *sourcedir, double *scposn)
    Return barycenter correction to TDB (in s)
    mjdTT:     Time of photon/pulse arrival in MJD(??)
    sourcedir: Direction cosines of source position
    scposn:    (Spacecraft) position vector in meters wrt geocenter

  double xbarycorr (double t, double *sourcedir, double *scposn)
    MET interface to barycorr
    t:         MET time of photon/pulse arrival
    sourcedir: Direction cosines of source position
    scposn:    Spacecraft position vector in meters

  void met2mjd (double t, MJDTime *mjdTT)
    met2mjd converts MET (seconds or days) to MJD day and fractional day.
    t:        MET time (s or d)
    mjdTT:    MJD (TT) day

  double mjd2met (MJDTime *mjdTT)
    mjd2met converts MJD day and fractional day to MET (seconds or days).
    mjdTT:    MJD (TT) day
    return:   MET time (s or d)

  double toTT (MJDTime *mjd)
    toTT converts mjd to TT
    mjd:      MJD (?) day, converted to MJD (TT)
    return:   Change in seconds

  double toUTC (MJDTime *mjd)
    toUTC converts mjd to UTC
    mjd:      MJD (?) day, converted to MJD (UTC)
    return:   Change in seconds

  const char *convertTime (MJDTime *mjd)
    convertTime converts mjd to a FITS-compliant DATE string
    mjd:      MJD (?) day
    returns YY/MM/DD or CCYY-MM-DDThh:mm:ss

  const char *fitsdate ()
    returns current UTC as CCYY-MM-DDThh:mm:ss

  int timeparms (char *insource, MJDTime *time, double roffset, int bin,
	         PsrTimPar *ptp, PsrBinPar *pbp, int debug)
    timeparms gets the timing parameters from the pulsar database.
    char *insource:  source name (J1234-5643, B1234-56, or 1234-56)
    MJDTime *time:   MJD(TB) time for which the parameters are requested
    double roffset:  extra offset to be applied to radio zero phase (in s)
    int bin:         binary pulsar?
    PsrTimPar *ptp:  pulsar timing parameters struct
    PsrBinPar *pbp:  pulsar binary orbit parameters struct
    int debug:       chatty

  double absphase (MJDTime *time, PsrTimPar *ptp, int binary, PsrBinPar *pbp,
        	   int bcorr, double *dir, char *oefile, int *error)
    absphase calculates absolute phase
    MJDTime *time:   MJD(TDB) time for which the parameters are requested
    PsrTimPar *ptp:  pulsar timing parameters struct
    int binary:      binary pulsar?
    PsrBinPar *pbp:  pulsar binary orbit parameters struct
    int bcorr:       apply barycenter correction?
    double dir[3]:   source position vector
    char *oefile:    orbit ephemeris filename
    int *error:      error (status) flag

  double xabsphase (double time,  PsrTimPar *ptp, int binary, PsrBinPar *pbp,
	            int bcorr, double *dir, char *oefile, int *error)
    MET interface to absphase

  void c200to405 (double *dir, int debug)
    c200to405 converts a direction cosine vector from the DE200 frame
    (nominally FK5, but not quite) to the DE405 frame (ICRS).
    double dir[3]    direction cosine vector (in and out)
    int debug        print the conversion?


scorbit.c

  void scorbitinit (enum Observatory)

  double *scorbit (char *filename, MJDTime *time, int *oerror)
    Return spacecraft position (in m, wrt geocenter)
    filename:  orbit ephemeris file
    time:      time for which the position is requested
    oerror:    non-zero if error occurred

  double *xscorbit (char *filename, double time, int *oerror)
    MET interface to scorbit
    filename:  orbit ephemeris file
    time:      time for which the position is requested
    oerror:    non-zero if error occurred


clock.c

  void clockinit (enum Observatory)

  double clockCorr (double time, double timezero, double *timeparms,
                    char *instrument)
    A function that determines the clock correction depending on mission.
      time         double   MET time
      timezero     double   Clock correction in header
      timeparms    double*  Clock parameters, the first of which is TIERABSO
      instrument   char*    INSTRUME (for instrument-dependent delays)
      timeparms[0] double*  Absolute clock error (<0 if unchanged)


phaseHist.c

  long** phaseHist (char *orbitFile, char *dataFile,
		    char *telescope, char *instrument,
		    char *radecsys, char *ephem,
		    char *source, double *ra, double *dec,
		    double roffset, double *dt0, int binary,
		    PsrTimPar *inptp, PsrBinPar *inpbp, int l1only,
		    int nTR, double *starts, double *stops,
		    double *tstart, double *tstop, double **exposure,
		    int accum, int *gainapp, char *datobs,
		    int *nphase, int *nchan, char *cpix, double *f0,
		    int nobcorr, int *nevts, int debug)
    phaseHist is a C functions that processes pulsar observations,
    calculating the phase and binning the events in phase/energy space.
    It requires bary.h, bary.c, and dpleph.c.
    The main interface is through phaseHist:
    long **phaseHist     Returns a two dimensional histogram
                           long image[*nchan][*nphase]
                           Upon entering phaseHist:
                             If the histogram exists :
                               If !accum :
                                 deallocate the histogram
                                 allocate a new histogram
                                 zero the histogram
                             Else :
                               allocate a new histogram
                               zero the histogram
                           If an error occurs, a NULL pointer is returned
      char *orbitFile    Path for orbit ephemeris file
      char *dataFile     Path for data file; supported data modes:
                           PCA GoodXenon
                           PCA Single-EA event mode, including CE_
                           PCA Pulsar fold mode
                           PCA Binned mode, including CB_, but
                                            excepting multi-PCU/layer
                           PCA Single-bit mode
                           HEXTE Event mode
    + char *telescope    Mission/telescope used for collecting data
    + char *instrument   Instrument used for collecting data
    + char *radecsys     FK5 or ICRS
    + char *ephem        JPL-DE200 or JPL-DE405
      char *source       Pulsar database source name; may be B or J name
                           It expects the radio pulsar database in file
                           $TIMING_DIR/psrtime.dat
    + double *ra         RA used for calculating the phase
    + double *dec        Dec used for calculating the phase
      double roffset     Additional radio offset (in s)
                           This is the dispersion measure correction
                           for Crab and Vela pulsars
      double *dt0        Additional clock correction (in s); return full corr.
      int binary         Binary pulsar?
                           It expects the binary orbit database in file
                           $TIMING_DIR/psrbin.dat.
      PsrTimPar *inptp   Input PTP, if radio database is not to be used
      PsrBinPar *inpbp   Input PBP, if radio database is not to be used
                           If radio DB is to be used, provide pointers
                           to struct with (inptp->ra < 0) or NULL pointers
                           These structs are defined in bary.h
      int l1only         Layer 1 only? (only for GoodXenon)
                           -1: only propane layer
                            0: all xenon layer
                            1: only layer 1
      int nTR            Number of time ranges to apply
      double *starts     double starts[nTR]: GTI start times
      double *stops      double stops[nTR]: GTI stop times
    + double *tstart     double tstart: start of data in MJD
    + double *tstop      double tstop: End of data in MJD
    + double **exposure  Total exposure time (always in s) for phase bin i
                           is in *exposure[i]
      int accum          Do not zero the histogram, but accumulate
    + int *gainapp       Value of GAINAPP keyword
    + char *datobs       Value of DATE-OBS keyword; only modified if != 0
   ++ int *nphase        Number of phase bins per period (default: 100)
   ++ int *nchan         Number of channel bins (default: 256)
   ++ char *cpix         Channel coordinate descriptor (default: none)
    + double *f0         Current frequency
      int nobcorr        Force no barycenter correction
    + int *nevts         Number of events accumulated
      int debug          Debug mode

      + output only
     ++ input/output



Time Conversion Recipes
-----------------------

The functions timeinit, met2mjd, toTT, toUTC, mjd2met, and convertTime
may be used to do some simple conversions and FITS-compliant
formatting:

  #include "bary.h"
  int main (int argc, char **argv) {
    char date[24] ;
    double t ;
    MJDTime mjd ;

For just formatting a MET (TT) time to FITS DATExxxx format:
    timeinit (AXAF, "TT", "TT", -1.0, 0.0, "s") ;
    met2mjd (t, &mjd) ;
    strcpy (date, convertTime(&mjd) ;

Converting a TT MET to a UTC DATExxxx string (not to be inserted into
FITS files with TIMESYS='TT'!):
    timeinit (AXAF, "TT", "UTC", -1.0, 0.0, "s") ;
    met2mjd (t, &mjd) ;
    toUTC (&mjd) ;
    strcpy (date, convertTime(&mjd) ;

Converting a MJD time to a DATExxxx string:
    timeinit (AXAF, "TT", "TT", 0.0, 0.0, "d") ;
    met2mjd (t, &mjd) ;
    strcpy (date, convertTime(&mjd) ;

Getting the current UTC date and time in FITS-compatible format:
    strcpy (date, fitsdate ()) ;
