Index: /trunk/FACT++/pal/.gitignore =================================================================== --- /trunk/FACT++/pal/.gitignore (revision 18711) +++ /trunk/FACT++/pal/.gitignore (revision 18712) @@ -9,4 +9,5 @@ *.o *.lo +*.lof *.la *.dSYM Index: /trunk/FACT++/pal/component.xml =================================================================== --- /trunk/FACT++/pal/component.xml (revision 18711) +++ /trunk/FACT++/pal/component.xml (revision 18712) @@ -4,5 +4,5 @@ - 0.9.3 + 0.9.7 libext/pal Positional Astronomy Library Index: /trunk/FACT++/pal/configure.ac =================================================================== --- /trunk/FACT++/pal/configure.ac (revision 18711) +++ /trunk/FACT++/pal/configure.ac (revision 18712) @@ -3,5 +3,5 @@ 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,4 +55,5 @@ AC_DEFINE( [HAVE_STAR_UTIL_H], [1], [Define to 1 if you have the header file]) else + AC_MSG_NOTICE([Building outside a Starlink environment]) # Allow ERFA/SOFA location to be specified using --with-erfa=$ERFA_DIR @@ -118,5 +119,5 @@ dnl Disable document building regardless of --without-stardocs - _star_build_docs=: + _star_build_docs=false fi Index: /trunk/FACT++/pal/pal.news =================================================================== --- /trunk/FACT++/pal/pal.news (revision 18711) +++ /trunk/FACT++/pal/pal.news (revision 18712) @@ -4,4 +4,20 @@ 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 Index: /trunk/FACT++/pal/palAmpqk.c =================================================================== --- /trunk/FACT++/pal/palAmpqk.c (revision 18711) +++ /trunk/FACT++/pal/palAmpqk.c (revision 18712) @@ -26,6 +26,6 @@ * (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,6 +44,14 @@ * 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,4 +60,6 @@ * 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,4 +67,5 @@ * Copyright (C) 2000 Rutherford Appleton Laboratory * Copyright (C) 2012 Science and Technology Facilities Council. +* Copyright (C) 2016 Tim Jenness * All Rights Reserved. @@ -91,9 +102,12 @@ 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,6 +137,23 @@ } +/* 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 ); } Index: /trunk/FACT++/pal/palMappa.c =================================================================== --- /trunk/FACT++/pal/palMappa.c (revision 18711) +++ /trunk/FACT++/pal/palMappa.c (revision 18712) @@ -26,7 +26,7 @@ * - (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 Index: /trunk/FACT++/pal/palMapqk.c =================================================================== --- /trunk/FACT++/pal/palMapqk.c (revision 18711) +++ /trunk/FACT++/pal/palMapqk.c (revision 18712) @@ -52,11 +52,24 @@ * 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,5 +114,5 @@ /* local constants */ - const double VF = 0.21094502; /* Km/s to AU/year */ + const double VF = 0.210945028; /* Km/s to AU/year */ /* Local Variables: */ Index: /trunk/FACT++/pal/palMapqkz.c =================================================================== --- /trunk/FACT++/pal/palMapqkz.c (revision 18711) +++ /trunk/FACT++/pal/palMapqkz.c (revision 18712) @@ -25,8 +25,8 @@ * 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,9 +50,22 @@ * 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,10 +113,13 @@ /* 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,10 +127,17 @@ 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] ) ); } Index: /trunk/FACT++/pal/palTest.c =================================================================== --- /trunk/FACT++/pal/palTest.c (revision 18711) +++ /trunk/FACT++/pal/palTest.c (revision 18712) @@ -276,6 +276,6 @@ /* 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,16 +1270,67 @@ } +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,6 +1339,6 @@ 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,4 +2079,5 @@ t_map(&status); t_mappa(&status); + t_mapqk(&status); t_mapqkz(&status); t_moon(&status); @@ -2055,3 +2107,2 @@ return status; } -