1 /* aberration, Jean Meeus, "Astronomical Algorithms", Willman-Bell, 1995;
2 * based on secular unperturbed Kepler orbit
3 *
4 * the corrections should be applied to ra/dec and lam/beta at the
5 * epoch of date.
6 */
7
8 #include <stdio.h>
9 #include <math.h>
10 #include <stdlib.h>
11
12 #include "astro.h"
13
14 #define ABERR_CONST (20.49552/3600./180.*PI) /* aberr const in rad */
15 #define AB_ECL_EOD 0
16 #define AB_EQ_EOD 1
17
18 static void ab_aux (double mj, double *x, double *y, double lsn, int mode);
19
20 /* apply aberration correction to ecliptical coordinates *lam and *bet
21 * (in radians) for a given time m and handily supplied longitude of sun,
22 * lsn (in radians)
23 */
24 void
ab_ecl(double mj,double lsn,double * lam,double * bet)25 ab_ecl (double mj, double lsn, double *lam, double *bet)
26 {
27 ab_aux(mj, lam, bet, lsn, AB_ECL_EOD);
28 }
29
30 /* apply aberration correction to equatoreal coordinates *ra and *dec
31 * (in radians) for a given time m and handily supplied longitude of sun,
32 * lsn (in radians)
33 */
34 void
ab_eq(double mj,double lsn,double * ra,double * dec)35 ab_eq (double mj, double lsn, double *ra, double *dec)
36 {
37 #if defined(USE_MEEUS_AB_EQ)
38
39 /* this claims to account for earth orbit excentricity and is also
40 * smooth clear to dec=90 but it does not work well backwards with
41 * ap_as()
42 */
43 ab_aux(mj, ra, dec, lsn, AB_EQ_EOD);
44
45 #else /* use Montenbruck */
46
47 /* this agrees with Meeus to within 0.2 arcsec until dec gets larger
48 * than about 89.9, then grows to 1as at 89.97. but it works very
49 * smoothly with ap_as
50 */
51 double x, y, z; /* equatorial rectangular coords */
52 double vx, vy, vz; /* aberration velocity in rectangular coords */
53 double L; /* helio long of earth */
54 double cL;
55 double r;
56
57
58 sphcart (*ra, *dec, 1.0, &x, &y, &z);
59
60 L = 2*PI*(0.27908 + 100.00214*(mj-J2000)/36525.0);
61 cL = cos(L);
62 vx = -0.994e-4*sin(L);
63 vy = 0.912e-4*cL;
64 vz = 0.395e-4*cL;
65 x += vx;
66 y += vy;
67 z += vz;
68
69 cartsph (x, y, z, ra, dec, &r);
70
71 #endif
72 }
73
74 /* because the e-terms are secular, keep the real transformation for both
75 * coordinate systems in here with the secular variables cached.
76 * mode == AB_ECL_EOD: x = lam, y = bet (ecliptical)
77 * mode == AB_EQ_EOD: x = ra, y = dec (equatoreal)
78 */
79 static void
ab_aux(double mj,double * x,double * y,double lsn,int mode)80 ab_aux (double mj, double *x, double *y, double lsn, int mode)
81 {
82 static double lastmj = -10000;
83 static double eexc; /* earth orbit excentricity */
84 static double leperi; /* ... and longitude of perihelion */
85 static char dirty = 1; /* flag for cached trig terms */
86
87 if (mj != lastmj) {
88 double T; /* centuries since J2000 */
89
90 T = (mj - J2000)/36525.;
91 eexc = 0.016708617 - (42.037e-6 + 0.1236e-6 * T) * T;
92 leperi = degrad(102.93735 + (0.71953 + 0.00046 * T) * T);
93 lastmj = mj;
94 dirty = 1;
95 }
96
97 switch (mode) {
98 case AB_ECL_EOD: /* ecliptical coords */
99 {
100 double *lam = x, *bet = y;
101 double dlsun, dlperi;
102
103 dlsun = lsn - *lam;
104 dlperi = leperi - *lam;
105
106 /* valid only for *bet != +-PI/2 */
107 *lam -= ABERR_CONST/cos(*bet) * (cos(dlsun) -
108 eexc*cos(dlperi));
109 *bet -= ABERR_CONST*sin(*bet) * (sin(dlsun) -
110 eexc*sin(dlperi));
111 }
112 break;
113
114 case AB_EQ_EOD: /* equatoreal coords */
115 {
116 double *ra = x, *dec = y;
117 double sr, cr, sd, cd, sls, cls;/* trig values coords */
118 static double cp, sp, ce, se; /* .. and perihel/eclipic */
119 double dra, ddec; /* changes in ra and dec */
120
121 if (dirty) {
122 double eps;
123
124 cp = cos(leperi);
125 sp = sin(leperi);
126 obliquity(mj, &eps);
127 se = sin(eps);
128 ce = cos(eps);
129 dirty = 0;
130 }
131
132 sr = sin(*ra);
133 cr = cos(*ra);
134 sd = sin(*dec);
135 cd = cos(*dec);
136 sls = sin(lsn);
137 cls = cos(lsn);
138
139 dra = ABERR_CONST/cd * ( -(cr * cls * ce + sr * sls) +
140 eexc * (cr * cp * ce + sr * sp));
141
142 ddec = se/ce * cd - sr * sd; /* tmp use */
143 ddec = ABERR_CONST * ( -(cls * ce * ddec + cr * sd * sls) +
144 eexc * (cp * ce * ddec + cr * sd * sp) );
145
146 *ra += dra;
147 *dec += ddec;
148 radecrange (ra, dec);
149 }
150 break;
151
152 default:
153 printf ("ab_aux: bad mode: %d\n", mode);
154 abort();
155 break;
156
157 } /* switch (mode) */
158 }
159
160