Changeset 18712 for trunk/FACT++

Timestamp:

01/06/17 11:51:27 (8 years ago)

Author:

tbretz

Message:

Updated to PAL 0.9.7 (adds mainly light deflection to abberation which was previously missing)

Location:

trunk/FACT++/pal

Files:

• .gitignore (modified) (1 diff)
• component.xml (modified) (1 diff)
• configure.ac (modified) (3 diffs)
• pal.news (modified) (1 diff)
• palAmpqk.c (modified) (6 diffs)
• palMappa.c (modified) (1 diff)
• palMapqk.c (modified) (2 diffs)
• palMapqkz.c (modified) (4 diffs)
• palTest.c (modified) (5 diffs)

trunk/FACT++/pal/.gitignore

@@ -9 +9 @@
 *.o
 *.lo
+*.lof
 *.la
 *.dSYM

trunk/FACT++/pal/component.xml

@@ -4 +4 @@
 
 <component id="pal" support="S">
-  <version>0.9.3</version>
+  <version>0.9.7</version>
   <path>libext/pal</path>
   <description>Positional Astronomy Library</description>

trunk/FACT++/pal/configure.ac

@@ -3 +3 @@
 
 dnl    Initialisation: package name and version number
-AC_INIT([pal],[0.9.3],[starlink@jiscmail.ac.uk])
+AC_INIT([pal],[0.9.7],[starlink@jiscmail.ac.uk])
 AC_CONFIG_AUX_DIR([build-aux])
 
@@ -55 +55 @@
   AC_DEFINE( [HAVE_STAR_UTIL_H], [1], [Define to 1 if you have the <star/util.h> header file])
 else
+  AC_MSG_NOTICE([Building outside a Starlink environment])
 
   #   Allow ERFA/SOFA location to be specified using --with-erfa=$ERFA_DIR
@@ -118 +119 @@
 
   dnl  Disable document building regardless of --without-stardocs
-  _star_build_docs=:
+  _star_build_docs=false
 fi
 

trunk/FACT++/pal/pal.news

@@ -4 +4 @@
 SLALIB API. It is distributed under the GPL and uses the SOFA library wherever
 possible.
+
+V0.9.7
+
+- Enable light deflection in palAmpqk.
+
+V0.9.5
+
+- Add light deflection to palAmpqk. SLALIB always had it but for some reason
+  the relevant piece of code never got ported to PAL.
+
+V0.9.4
+
+- Add light deflection to palMapqkz (thanks to Scott Daniel)
+- Add test for palMapqk (thanks to Scott Daniel)
+- Correctly disable documentation build when running outside of a Starlink
+  environment.
 
 V0.9.3

trunk/FACT++/pal/palAmpqk.c

@@ -26 +26 @@
 *        (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)
@@ -44 +44 @@
 *     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}
 
@@ -52 +60 @@
 *        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}
 
@@ -57 +67 @@
 *     Copyright (C) 2000 Rutherford Appleton Laboratory
 *     Copyright (C) 2012 Science and Technology Facilities Council.
+*     Copyright (C) 2016 Tim Jenness
 *     All Rights Reserved.
 
@@ -91 +102 @@
    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];
    }
@@ -123 +137 @@
    }
 
+/* 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 );
 }

trunk/FACT++/pal/palMappa.c

@@ -26 +26 @@
 *        - (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
 

trunk/FACT++/pal/palMapqk.c

@@ -52 +52 @@
 *     routine can be used instead.
 
+*  Notes:
+*     - The reference frames and timescales used are post 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).
+*     - 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: Patrick T. Wallace
 *     TIMJ: Tim Jenness (JAC, Hawaii)
 *     {enter_new_authors_here}
-
-*  Notes:
-*     - The reference frames and timescales used are post IAU 2006.
 
 *  History:
@@ -101 +114 @@
 
 /* local constants */
-   const double VF = 0.21094502; /* Km/s to AU/year */
+   const double VF = 0.210945028; /* Km/s to AU/year */
 
 /* Local Variables: */

trunk/FACT++/pal/palMapqkz.c

@@ -25 +25 @@
 *        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)
@@ -50 +50 @@
 *     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:
@@ -100 +113 @@
 /* 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];
    }
 
@@ -111 +127 @@
    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] ) );
    }
 

trunk/FACT++/pal/palTest.c

@@ -276 +276 @@
   /* This is the palMapqk test */
   palAmp( 1.234, -0.567, 55927.0, 2010.0, &rm, &dm );
-  vvd( rm, 1.2335120411026936349, 1.0E-12, "palAmp", "rm", status );
-  vvd( dm, -0.56702908706930343907, 1.0E-12, "palAmp", "dm", status );
+  vvd( rm, 1.233512033578303857, 1.0E-12, "palAmp", "rm", status );
+  vvd( dm, -0.56702909748530827549, 1.0E-12, "palAmp", "dm", status );
 }
 
@@ -1270 +1270 @@
 }
 
+static void t_mapqk( int *status ) {
+    /* Test mapqk by taking the geocentric apparent positions of Arcturus
+       as downloaded from aa.usno.mil/data/docs/geocentric.php and trying
+       to calculate it from Arcturus' mean position, proper motion, parallax,
+       and radial velocity */
+
+    double amprms[21];
+    double ra_0, dec_0; /* mean position */
+    double ra_app, dec_app; /* geocentric apparent position */
+    double ra_test, dec_test;
+    double px, pm_ra, pm_dec, v_rad;
+    int j;
+    int iposn;
+
+    /*
+    The data below represents the position of Arcturus on
+    JD (UT) 2457000.375 as reported by
+    http://aa.usno.navy.mil/data/docs/geocentric.php
+    */
+
+    const char radec_0[] = "14 15 39.67207 19 10 56.673";
+    const char radec_app[] = "14 16 19.59 19 6 19.56";
+
+    iposn = 1;
+    palDafin(radec_0, &iposn, &ra_0, &j);
+    palDafin(radec_0, &iposn, &dec_0, &j);
+    ra_0 *= 15.0;
+
+    pm_ra = -1.0939*PAL__DAS2R;
+    pm_ra /= cos(dec_0);
+    pm_dec = -2.00006*PAL__DAS2R;
+    v_rad = -5.19;
+    px = 0.08883*PAL__DAS2R;
+
+    palMappa(2000.0, 56999.87537249177, amprms);  /* time is the TDB MJD calculated from
+                                                     a JD of 2457000.375 with astropy.time */
+
+    palMapqk(ra_0, dec_0, pm_ra, pm_dec, px, v_rad, amprms, &ra_test, &dec_test);
+
+    iposn = 1;
+    palDafin(radec_app, &iposn, &ra_app, &j);
+    palDafin(radec_app, &iposn, &dec_app, &j);
+    ra_app *= 15.0;
+
+    /* find the angular distance from the known mean position
+       to the calculated mean postion */
+    double dd;
+    dd = palDsep(ra_test, dec_test, ra_app, dec_app);
+    dd *= PAL__DR2AS;
+    vvd( dd, 0.0, 0.1, "palMapqk", "distance", status);
+}
+
 static void t_mapqkz( int *status ) {
-   double amprms[21],  ra, da;
+   double amprms[21],  ra, da, ra_c, da_c;
+
+   /* Run inputs through mapqk with zero proper motion, parallax
+      and radial velocity.  Then run the same inputs through mapqkz.
+      Verify that the results are the same */
    palMappa( 2010.0, 55927.0, amprms );
+   palMapqk(1.234, -0.567, 0.0, 0.0, 0.0, 0.0, amprms, &ra_c, &da_c);
    palMapqkz( 1.234, -0.567, amprms, &ra, &da );
-   vvd( ra, 1.2344879748414849807, 1.0E-12, "palMapqkz", "ra", status );
-   vvd( da, -0.56697099554368701746, 1.0E-12, "palMapqkz", "da", status );
-
-   /* Try the same with palMapqk and zero parallax and proper motion */
-   palMapqk( 1.234, -0.567, 0., 0., 0., 0., amprms, &ra, &da );
-   vvd( ra, 1.2344879748414849807, 1.0E-7, "palMapqkz", "ra", status );
-   vvd( da, -0.56697099554368701746, 1.0E-7, "palMapqkz", "da", status );
-
+   vvd( ra, ra_c, 1.0E-12, "palMapqkz", "ra", status );
+   vvd( da, da_c, 1.0E-12, "palMapqkz", "da", status );
 }
 
@@ -1288 +1339 @@
    palMappa( 2010.0, 55927.0, amprms );
    palAmpqk( 1.234, -0.567, amprms, &rm, &dm );
-   vvd( rm, 1.2335120411026936349, 1.0E-12, "palAmpqk", "rm", status );
-   vvd( dm, -0.56702908706930343907, 1.0E-12, "palAmpqk", "dm", status );
+   vvd( rm, 1.233512033578303857, 1.0E-12, "palAmpqk", "rm", status );
+   vvd( dm, -0.56702909748530827549, 1.0E-12, "palAmpqk", "dm", status );
 }
 
@@ -2028 +2079 @@
   t_map(&status);
   t_mappa(&status);
+  t_mapqk(&status);
   t_mapqkz(&status);
   t_moon(&status);
@@ -2055 +2107 @@
   return status;
 }
-