/* Precession of the equinox and ecliptic * from epoch Julian date J to or from J2000.0 * * Program by Steve Moshier. */ #define WILLIAMS 1 /* James G. Williams, "Contributions to the Earth's obliquity rate, precession, and nutation," Astron. J. 108, 711-724 (1994) */ #define SIMON 0 /* J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze', G. Francou, and J. Laskar, "Numerical Expressions for precession formulae and mean elements for the Moon and the planets," Astronomy and Astrophysics 282, 663-683 (1994) */ #define IAU 0 /* IAU Coefficients are from: * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando, * "Expressions for the Precession Quantities Based upon the IAU * (1976) System of Astronomical Constants," Astronomy and * Astrophysics 58, 1-16 (1977). */ #define LASKAR 0 /* Newer formulas that cover a much longer time span are from: * J. Laskar, "Secular terms of classical planetary theories * using the results of general theory," Astronomy and Astrophysics * 157, 59070 (1986). * * See also: * P. Bretagnon and G. Francou, "Planetary theories in rectangular * and spherical variables. VSOP87 solutions," Astronomy and * Astrophysics 202, 309-315 (1988). * * Laskar's expansions are said by Bretagnon and Francou * to have "a precision of about 1" over 10000 years before * and after J2000.0 in so far as the precession constants p^0_A * and epsilon^0_A are perfectly known." * * Bretagnon and Francou's expansions for the node and inclination * of the ecliptic were derived from Laskar's data but were truncated * after the term in T**6. I have recomputed these expansions from * Laskar's data, retaining powers up to T**10 in the result. * * The following table indicates the differences between the result * of the IAU formula and Laskar's formula using four different test * vectors, checking at J2000 plus and minus the indicated number * of years. * * Years Arc * from J2000 Seconds * ---------- ------- * 0 0 * 100 .006 * 200 .006 * 500 .015 * 1000 .28 * 2000 6.4 * 3000 38. * 10000 9400. */ #define DOUBLE double #if __STDC__ double cos (double); double sin (double); int epsiln (double); #else double cos(), sin(); extern int epsiln(); #endif #define COS cos #define SIN sin extern DOUBLE J2000; /* = 2451545.0, 2000 January 1.5 */ extern DOUBLE STR; /* = 4.8481368110953599359e-6 radians per arc second */ extern DOUBLE coseps, sineps; /* see epsiln.c */ /* In WILLIAMS and SIMON, Laskar's terms of order higher than t^4 have been retained, because Simon et al mention that the solution is the same except for the lower order terms. */ #if WILLIAMS static DOUBLE pAcof[] = { #if 1 /* Corrections to Williams (1994) introduced in DE403. */ -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3, -0.235316, 0.076, 110.5414, 50287.91959 #else -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3, -0.235316, 0.076, 110.5407, 50287.70000 #endif }; #endif #if SIMON /* Precession coefficients from Simon et al: */ static DOUBLE pAcof[] = { -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3, -0.235316, 0.07732, 111.2022, 50288.200 }; #endif #if LASKAR /* Precession coefficients taken from Laskar's paper: */ static DOUBLE pAcof[] = { -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3, -0.235316, 0.07732, 111.1971, 50290.966 }; #endif #if WILLIAMS /* Pi from Williams' 1994 paper, in radians. No change in DE403. */ static DOUBLE nodecof[] = { 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10, -3.54e-9, -1.8103e-7, 1.26e-7, 7.436169e-5, -0.04207794833, 3.052115282424}; /* pi from Williams' 1994 paper, in radians. No change in DE403. */ static DOUBLE inclcof[] = { 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, -5.4000441e-11, 1.32115526e-9, -6.012e-7, -1.62442e-5, 0.00227850649, 0.0 }; #endif #if SIMON static DOUBLE nodecof[] = { 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10, -3.54e-9, -1.8103e-7, 2.579e-8, 7.4379679e-5, -0.0420782900, 3.0521126906}; static DOUBLE inclcof[] = { 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, -5.4000441e-11, 1.32115526e-9, -5.99908e-7, -1.624383e-5, 0.002278492868, 0.0 }; #endif #if LASKAR /* Node and inclination of the earth's orbit computed from * Laskar's data as done in Bretagnon and Francou's paper. * Units are radians. */ static DOUBLE nodecof[] = { 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10, -3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5, -0.042078604317, 3.052112654975 }; static DOUBLE inclcof[] = { 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, -5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5, 0.002278495537, 0.0 }; #endif /* Subroutine arguments: * * R = rectangular equatorial coordinate vector to be precessed. * The result is written back into the input vector. * J = Julian date * direction = * Precess from J to J2000: direction = 1 * Precess from J2000 to J: direction = -1 * Note that if you want to precess from J1 to J2, you would * first go from J1 to J2000, then call the program again * to go from J2000 to J2. */ int precess( R, J, direction ) DOUBLE R[], J; int direction; { DOUBLE A, B, T, pA, W, z; DOUBLE x[3]; DOUBLE *p; int i; #if IAU DOUBLE sinth, costh, sinZ, cosZ, sinz, cosz, Z, TH; #endif if( J == J2000 ) return(0); /* Each precession angle is specified by a polynomial in * T = Julian centuries from J2000.0. See AA page B18. */ T = (J - J2000)/36525.0; #if IAU /* Use IAU formula only for a few centuries, if at all. */ if( FABS(T) > Two ) goto laskar; Z = (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR; z = (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR; TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR; sinth = SIN(TH); costh = COS(TH); sinZ = SIN(Z); cosZ = COS(Z); sinz = SIN(z); cosz = COS(z); A = cosZ*costh; B = sinZ*costh; if( direction < 0 ) { /* From J2000.0 to J */ x[0] = (A*cosz - sinZ*sinz)*R[0] - (B*cosz + cosZ*sinz)*R[1] - sinth*cosz*R[2]; x[1] = (A*sinz + sinZ*cosz)*R[0] - (B*sinz - cosZ*cosz)*R[1] - sinth*sinz*R[2]; x[2] = cosZ*sinth*R[0] - sinZ*sinth*R[1] + costh*R[2]; } else { /* From J to J2000.0 */ x[0] = (A*cosz - sinZ*sinz)*R[0] + (A*sinz + sinZ*cosz)*R[1] + cosZ*sinth*R[2]; x[1] = -(B*cosz + cosZ*sinz)*R[0] - (B*sinz - cosZ*cosz)*R[1] - sinZ*sinth*R[2]; x[2] = -sinth*cosz*R[0] - sinth*sinz*R[1] + costh*R[2]; } goto done; laskar: #endif /* IAU */ /* Implementation by elementary rotations using Laskar's expansions. * First rotate about the x axis from the initial equator * to the ecliptic. (The input is equatorial.) */ if( direction == 1 ) epsiln( J ); /* To J2000 */ else epsiln( J2000 ); /* From J2000 */ x[0] = R[0]; z = coseps*R[1] + sineps*R[2]; x[2] = -sineps*R[1] + coseps*R[2]; x[1] = z; /* Precession in longitude */ T /= 10.0; /* thousands of years */ p = pAcof; pA = *p++; for( i=0; i<9; i++ ) pA = pA * T + *p++; pA *= STR * T; /* Node of the moving ecliptic on the J2000 ecliptic. */ p = nodecof; W = *p++; for( i=0; i<10; i++ ) W = W * T + *p++; /* Rotate about z axis to the node. */ if( direction == 1 ) z = W + pA; else z = W; B = COS(z); A = SIN(z); z = B * x[0] + A * x[1]; x[1] = -A * x[0] + B * x[1]; x[0] = z; /* Rotate about new x axis by the inclination of the moving * ecliptic on the J2000 ecliptic. */ p = inclcof; z = *p++; for( i=0; i<10; i++ ) z = z * T + *p++; if( direction == 1 ) z = -z; B = COS(z); A = SIN(z); z = B * x[1] + A * x[2]; x[2] = -A * x[1] + B * x[2]; x[1] = z; /* Rotate about new z axis back from the node. */ if( direction == 1 ) z = -W; else z = -W - pA; B = COS(z); A = SIN(z); z = B * x[0] + A * x[1]; x[1] = -A * x[0] + B * x[1]; x[0] = z; /* Rotate about x axis to final equator. */ if( direction == 1 ) epsiln( J2000 ); else epsiln( J ); z = coseps * x[1] - sineps * x[2]; x[2] = sineps * x[1] + coseps * x[2]; x[1] = z; #if IAU done: #endif for( i=0; i<3; i++ ) R[i] = x[i]; return(0); }