Changeset 18712 for trunk/FACT++/pal/palMapqkz.c

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)

File:

• trunk/FACT++/pal/palMapqkz.c (modified) (4 diffs)

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] ) );
    }