Update PAL files to PAL V0.9.7

ast.news

@@ -19,6 +19,9 @@ environment-independent.
 Main Changes in this Version
 ----------------------------
 
+- The PAL library files included in the AST distribution have been updated
+to PAL version 0.9.7.
+
 - Multiple identical NormMaps in series will now be simplified to a
 single NormMap.
 

pal.h

@@ -31,6 +31,8 @@
 *        If the PAL and ERFA libraries distributed within AST are being
 *        used, rename the public symbols defined by both to avoid clashes
 *        with any external PAL or ERFA library.
+*     11-JAN-2017 (DSB):
+*        Update to PAL V0.9.7.
 *     {enter_further_changes_here}
 
 *  Copyright:
@@ -118,7 +120,7 @@ void palAoppa ( double date, double dut, double elongm, double phim,
 
 void palAoppat ( double date, double aoprms[14] );
 
-void palAopqk ( double rap, double dap, double aoprms[14],
+void palAopqk ( double rap, double dap, const double aoprms[14],
                 double *aob, double *zob, double *hob,
                 double *dob, double *rob );
 
@@ -396,7 +398,7 @@ void palOap ( const char *type, double ob1, double ob2, double date,
               double rh, double wl, double tlr,
               double *rap, double *dap );
 
-void palOapqk ( const char *type, double ob1, double ob2, double aoprms[14],
+void palOapqk ( const char *type, double ob1, double ob2, const double aoprms[14],
                 double *rap, double *dap );
 
 int palObs( size_t n, const char * c,
@@ -557,6 +559,8 @@ void palV2tp ( float v[3], float v0[3], float *xi, float *eta, int *j );
 
 float palVdv ( float va[3], float vb[3] );
 
+int palVers ( char * verstring, size_t verlen );
+
 void palVn ( float v[3], float uv[3], float *vm );
 
 void palVxv ( float va[3], float vb[3], float vc[3] );

pal/pal.h

@@ -533,6 +533,8 @@ void palV2tp ( float v[3], float v0[3], float *xi, float *eta, int *j );
 
 float palVdv ( float va[3], float vb[3] );
 
+int palVers ( char * verstring, size_t verlen );
+
 void palVn ( float v[3], float uv[3], float *vm );
 
 void palVxv ( float va[3], float vb[3], float vc[3] );

pal/palAmpqk.c

@@ -25,8 +25,8 @@
 *        Star-independent mean-to-apparent parameters (see palMappa):
 *        (0)      time interval for proper motion (Julian years)
 *        (1-3)    barycentric position of the Earth (AU)
-*        (4-6)    not used
-*        (7)      not used
+*        (4-6)    heliocentric direction of the Earth (unit vector)
+*        (7)      (grav rad Sun)*2/(Sun-Earth distance)
 *        (8-10)   abv: barycentric Earth velocity in units of c
 *        (11)     sqrt(1-v*v) where v=modulus(abv)
 *        (12-20)  precession/nutation (3,3) matrix
@@ -43,19 +43,30 @@
 *     star-independent parameters can be obtained by calling the palMappa
 *     function.
 
+*  Note:
+*     Iterative techniques are used for the aberration and
+*     light deflection corrections so that the routines
+*     palAmp (or palAmpqk) and palMap (or palMapqk) are
+*     accurate inverses;  even at the edge of the Sun's disc
+*     the discrepancy is only about 1 nanoarcsecond.
+
 *  Authors:
 *     PTW: Pat Wallace (STFC)
+*     TIMJ: Tim Jenness
 *     {enter_new_authors_here}
 
 *  History:
 *     2012-02-13 (PTW):
 *        Initial version.
 *        Adapted with permission from the Fortran SLALIB library.
+*     2016-12-19 (TIMJ):
+*        Add in light deflection (was missed in the initial port).
 *     {enter_further_changes_here}
 
 *  Copyright:
 *     Copyright (C) 2000 Rutherford Appleton Laboratory
 *     Copyright (C) 2012 Science and Technology Facilities Council.
+*     Copyright (C) 2016 Tim Jenness
 *     All Rights Reserved.
 
 *  Licence:
@@ -90,11 +101,14 @@ void palAmpqk ( double ra, double da, double amprms[21], double *rm,
    double abv[3];              /* Earth velocity wrt SSB (c, FK5) */
    double p1[3], p2[3], p3[3]; /* work vectors */
    double ab1p1, p1dv, p1dvp1, w;
+   double gr2e, pde, pdep1, ehn[3], p[3];
    int i, j;
 
 /* Unpack some of the parameters */
+   gr2e = amprms[7];
    ab1  = amprms[11];
    for( i = 0; i < 3; i++ ) {
+      ehn[i] = amprms[i + 4];
       abv[i] = amprms[i + 8];
    }
 
@@ -122,7 +136,24 @@ void palAmpqk ( double ra, double da, double amprms[21], double *rm,
       }
    }
 
+/* Light deflection */
+   for( i = 0; i < 3; i++ ) {
+       p[i] = p1[i];
+   }
+   for( j = 0; j < 5; j++ ) {
+      pde = eraPdp( p, ehn );
+      pdep1 = 1.0 + pde;
+      w = pdep1 - gr2e*pde;
+      for( i = 0; i < 3; i++ ) {
+         p[i] = (pdep1*p1[i] - gr2e*ehn[i])/w;
+      }
+      eraPn( p, &w, p2 );
+      for( i = 0; i < 3; i++ ) {
+         p[i] = p2[i];
+      }
+   }
+
 /* Mean RA,Dec */
-   eraC2s( p1, rm, dm );
+   eraC2s( p, rm, dm );
    *rm = eraAnp( *rm );
 }

pal/palDat.c

@@ -1,10 +1,10 @@
 /*
 *+
 *  Name:
-*     palDtt
+*     palDat
 
 *  Purpose:
-*     Return offset between UTC and TT
+*     Return offset between UTC and TAI
 
 *  Language:
 *     Starlink ANSI C

pal/palFk524.c

@@ -142,7 +142,7 @@ void palFk524( double r2000, double d2000, double dr2000, double dd2000,
    int i, j;
 
 /* Small number to avoid arithmetic problems. */
-   static const double tiny = 1.0-30;
+   static const double tiny = 1.0E-30;
 
 /* Canonical constants (see references). Constant vector and matrix. */
    double a[ 6 ] = { -1.62557E-6,  -0.31919E-6, -0.13843E-6,

pal/palGalsup.c

@@ -41,7 +41,7 @@
 *       machine-readable version of the above catalogue,
 *       Contract NAS 5-26490.
 *
-*    (These two references give different values for the galactic
+*     (These two references give different values for the galactic
 *     longitude of the supergalactic origin.  Both are wrong;  the
 *     correct value is L2=137.37.)
 

pal/palMappa.c

@@ -25,9 +25,9 @@
 *        - (0)      time interval for proper motion (Julian years)
 *        - (1-3)    barycentric position of the Earth (AU)
 *        - (4-6)    heliocentric direction of the Earth (unit vector)
-*        - (7)      (grav rad Sun)*2/(Sun-Earth distance)
+*        - (7)      (Schwarzschild radius of Sun)/(Sun-Earth distance)
 *        - (8-10)   abv: barycentric Earth velocity in units of c
-*        - (11)     sqrt(1-v**2) where v=modulus(abv)
+*        - (11)     sqrt(1-v^2) where v=modulus(abv)
 *        - (12-20)  precession/nutation (3,3) matrix
 
 *  Description:

pal/palMapqkz.c

@@ -4,7 +4,7 @@
 *     palMapqkz
 
 *  Purpose:
-*     Quick mean to apparent place.
+*     Quick mean to apparent place (no proper motion or parallax).
 
 *  Language:
 *     Starlink ANSI C
@@ -24,10 +24,10 @@
 *     amprms = double[21] (Given)
 *        Star-independent mean-to-apparent parameters (see palMappa):
 *        (0-3)    not used
-*        (4-6)    not used
+*        (4-6)    heliocentric direction of the Earth (unit vector)
 *        (7)      not used
 *        (8-10)   abv: barycentric Earth velocity in units of c
-*        (11)     sqrt(1-v**2) where v=modulus(abv)
+*        (11)     sqrt(1-v^2) where v=modulus(abv)
 *        (12-20)  precession/nutation (3,3) matrix
 *     ra = double * (Returned)
 *        Apparent RA (radians).
@@ -49,11 +49,24 @@
 *
 *     The corresponding function for the case of non-zero parallax
 *     and proper motion is palMapqk.
-*
-*     The reference systems and timescales used are IAU 2006.
-*
-*     Strictly speaking, the function is not valid for solar-system
-*     sources, though the error will usually be extremely small.
+
+*  Notes:
+*     - The reference systems and timescales used are IAU 2006.
+*     - The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
+*       are referred to the mean equinox and equator of the epoch
+*       specified when generating the precession/nutation matrix
+*       amprms[12-20].  In the call to palMappa (q.v.) normally used
+*       to populate amprms, this epoch is the first argument (eq).
+*     - The vector amprms(4-6) is referred to the mean equinox and
+*       equator of epoch eq.
+*     - Strictly speaking, the routine is not valid for solar-system
+*       sources, though the error will usually be extremely small.
+*       However, to prevent gross errors in the case where the
+*       position of the Sun is specified, the gravitational
+*       deflection term is restrained within about 920 arcsec of the
+*       centre of the Sun's disc.  The term has a maximum value of
+*       about 1.85 arcsec at this radius, and decreases to zero as
+*       the centre of the disc is approached.
 
 *  Authors:
 *     PTW: Pat Wallace (STFC)
@@ -99,23 +112,33 @@ void palMapqkz ( double rm, double dm, double amprms[21], double *ra,
 
 /* Local Variables: */
    int i;
-   double ab1, abv[3], p[3], w, p1dv, p1dvp1, p2[3], p3[3];
+   double ab1, abv[3], p[3], w, p1dv, p2[3], p3[3];
+   double gr2e, pde, pdep1, ehn[3], p1[3];
 
 /* Unpack scalar and vector parameters. */
    ab1 = amprms[11];
+   gr2e = amprms[7];
    for( i = 0; i < 3; i++ ) {
       abv[i] = amprms[i+8];
+      ehn[i] = amprms[i+4];
    }
 
 /* Spherical to x,y,z. */
    eraS2c( rm, dm, p );
 
+/* Light deflection (restrained within the Sun's disc) */
+   pde = eraPdp( p, ehn );
+   pdep1 = pde + 1.0;
+   w = gr2e / ( pdep1 > 1.0e-5 ? pdep1 : 1.0e-5 );
+   for( i = 0; i < 3; i++) {
+      p1[i] = p[i] + w * ( ehn[i] - pde * p[i] );
+   }
+
 /* Aberration. */
-   p1dv = eraPdp( p, abv );
-   p1dvp1 = p1dv + 1.0;
+   p1dv = eraPdp( p1, abv );
    w = 1.0 + p1dv / ( ab1 + 1.0 );
    for( i = 0; i < 3; i++ ) {
-      p2[i] = ( ( ab1 * p[i] ) + ( w * abv[i] ) ) / p1dvp1;
+      p2[i] = ( ( ab1 * p1[i] ) + ( w * abv[i] ) );
    }
 
 /* Precession and nutation. */