| 1 | /* nutation (in IAU (1980) expression) and abberation; stern
 | 
|---|
| 2 |  * on an HP PA processor, this reproduces the Almanac nutation values
 | 
|---|
| 3 |  * (given to 0.001") EXACTLY over 750 days (1995 and 1996)
 | 
|---|
| 4 |  */
 | 
|---|
| 5 | #include <stdio.h>
 | 
|---|
| 6 | #include <math.h>
 | 
|---|
| 7 | 
 | 
|---|
| 8 | #include "astro.h"
 | 
|---|
| 9 | 
 | 
|---|
| 10 | #define NUT_SCALE       1e4
 | 
|---|
| 11 | #define NUT_SERIES      106
 | 
|---|
| 12 | #define NUT_MAXMUL      4
 | 
|---|
| 13 | #define SECPERCIRC      (3600.*360.)
 | 
|---|
| 14 | 
 | 
|---|
| 15 | /* Delaunay arguments, in arc seconds; they differ slightly from ELP82B */
 | 
|---|
| 16 | static double delaunay[5][4] = {
 | 
|---|
| 17 |     {485866.733,  1717915922.633, 31.310,  0.064}, /* M', moon mean anom */
 | 
|---|
| 18 |     {1287099.804, 129596581.224,  -0.577, -0.012}, /* M, sun mean anom */
 | 
|---|
| 19 |     {335778.877,  1739527263.137, -13.257, 0.011}, /* F, moon arg lat */
 | 
|---|
| 20 |     {1072261.307, 1602961601.328, -6.891,  0.019}, /* D, elong moon sun */
 | 
|---|
| 21 |     {450160.280,  -6962890.539,   7.455,   0.008}, /* Om, moon l asc node */
 | 
|---|
| 22 | };
 | 
|---|
| 23 | 
 | 
|---|
| 24 | /* multipliers for Delaunay arguments */
 | 
|---|
| 25 | static short multarg[NUT_SERIES][5] = {
 | 
|---|
| 26 |         /* bounds:  -2..3, -2..2, -2/0/2/4, -4..4, 0..2 */
 | 
|---|
| 27 |     {0, 0, 0, 0, 1},
 | 
|---|
| 28 |     {0, 0, 0, 0, 2},
 | 
|---|
| 29 |     {-2, 0, 2, 0, 1},
 | 
|---|
| 30 |     {2, 0, -2, 0, 0},
 | 
|---|
| 31 |     {-2, 0, 2, 0, 2},
 | 
|---|
| 32 |     {1, -1, 0, -1, 0},
 | 
|---|
| 33 |     {0, -2, 2, -2, 1},
 | 
|---|
| 34 |     {2, 0, -2, 0, 1},
 | 
|---|
| 35 |     {0, 0, 2, -2, 2},
 | 
|---|
| 36 |     {0, 1, 0, 0, 0},
 | 
|---|
| 37 |     {0, 1, 2, -2, 2},
 | 
|---|
| 38 |     {0, -1, 2, -2, 2},
 | 
|---|
| 39 |     {0, 0, 2, -2, 1},
 | 
|---|
| 40 |     {2, 0, 0, -2, 0},
 | 
|---|
| 41 |     {0, 0, 2, -2, 0},
 | 
|---|
| 42 |     {0, 2, 0, 0, 0},
 | 
|---|
| 43 |     {0, 1, 0, 0, 1},
 | 
|---|
| 44 |     {0, 2, 2, -2, 2},
 | 
|---|
| 45 |     {0, -1, 0, 0, 1},
 | 
|---|
| 46 |     {-2, 0, 0, 2, 1},
 | 
|---|
| 47 |     {0, -1, 2, -2, 1},
 | 
|---|
| 48 |     {2, 0, 0, -2, 1},
 | 
|---|
| 49 |     {0, 1, 2, -2, 1},
 | 
|---|
| 50 |     {1, 0, 0, -1, 0},
 | 
|---|
| 51 |     {2, 1, 0, -2, 0},
 | 
|---|
| 52 |     {0, 0, -2, 2, 1},
 | 
|---|
| 53 |     {0, 1, -2, 2, 0},
 | 
|---|
| 54 |     {0, 1, 0, 0, 2},
 | 
|---|
| 55 |     {-1, 0, 0, 1, 1},
 | 
|---|
| 56 |     {0, 1, 2, -2, 0},
 | 
|---|
| 57 |     {0, 0, 2, 0, 2},
 | 
|---|
| 58 |     {1, 0, 0, 0, 0},
 | 
|---|
| 59 |     {0, 0, 2, 0, 1},
 | 
|---|
| 60 |     {1, 0, 2, 0, 2},
 | 
|---|
| 61 |     {1, 0, 0, -2, 0},
 | 
|---|
| 62 |     {-1, 0, 2, 0, 2},
 | 
|---|
| 63 |     {0, 0, 0, 2, 0},
 | 
|---|
| 64 |     {1, 0, 0, 0, 1},
 | 
|---|
| 65 |     {-1, 0, 0, 0, 1},
 | 
|---|
| 66 |     {-1, 0, 2, 2, 2},
 | 
|---|
| 67 |     {1, 0, 2, 0, 1},
 | 
|---|
| 68 |     {0, 0, 2, 2, 2},
 | 
|---|
| 69 |     {2, 0, 0, 0, 0},
 | 
|---|
| 70 |     {1, 0, 2, -2, 2},
 | 
|---|
| 71 |     {2, 0, 2, 0, 2},
 | 
|---|
| 72 |     {0, 0, 2, 0, 0},
 | 
|---|
| 73 |     {-1, 0, 2, 0, 1},
 | 
|---|
| 74 |     {-1, 0, 0, 2, 1},
 | 
|---|
| 75 |     {1, 0, 0, -2, 1},
 | 
|---|
| 76 |     {-1, 0, 2, 2, 1},
 | 
|---|
| 77 |     {1, 1, 0, -2, 0},
 | 
|---|
| 78 |     {0, 1, 2, 0, 2},
 | 
|---|
| 79 |     {0, -1, 2, 0, 2},
 | 
|---|
| 80 |     {1, 0, 2, 2, 2},
 | 
|---|
| 81 |     {1, 0, 0, 2, 0},
 | 
|---|
| 82 |     {2, 0, 2, -2, 2},
 | 
|---|
| 83 |     {0, 0, 0, 2, 1},
 | 
|---|
| 84 |     {0, 0, 2, 2, 1},
 | 
|---|
| 85 |     {1, 0, 2, -2, 1},
 | 
|---|
| 86 |     {0, 0, 0, -2, 1},
 | 
|---|
| 87 |     {1, -1, 0, 0, 0},
 | 
|---|
| 88 |     {2, 0, 2, 0, 1},
 | 
|---|
| 89 |     {0, 1, 0, -2, 0},
 | 
|---|
| 90 |     {1, 0, -2, 0, 0},
 | 
|---|
| 91 |     {0, 0, 0, 1, 0},
 | 
|---|
| 92 |     {1, 1, 0, 0, 0},
 | 
|---|
| 93 |     {1, 0, 2, 0, 0},
 | 
|---|
| 94 |     {1, -1, 2, 0, 2},
 | 
|---|
| 95 |     {-1, -1, 2, 2, 2},
 | 
|---|
| 96 |     {-2, 0, 0, 0, 1},
 | 
|---|
| 97 |     {3, 0, 2, 0, 2},
 | 
|---|
| 98 |     {0, -1, 2, 2, 2},
 | 
|---|
| 99 |     {1, 1, 2, 0, 2},
 | 
|---|
| 100 |     {-1, 0, 2, -2, 1},
 | 
|---|
| 101 |     {2, 0, 0, 0, 1},
 | 
|---|
| 102 |     {1, 0, 0, 0, 2},
 | 
|---|
| 103 |     {3, 0, 0, 0, 0},
 | 
|---|
| 104 |     {0, 0, 2, 1, 2},
 | 
|---|
| 105 |     {-1, 0, 0, 0, 2},
 | 
|---|
| 106 |     {1, 0, 0, -4, 0},
 | 
|---|
| 107 |     {-2, 0, 2, 2, 2},
 | 
|---|
| 108 |     {-1, 0, 2, 4, 2},
 | 
|---|
| 109 |     {2, 0, 0, -4, 0},
 | 
|---|
| 110 |     {1, 1, 2, -2, 2},
 | 
|---|
| 111 |     {1, 0, 2, 2, 1},
 | 
|---|
| 112 |     {-2, 0, 2, 4, 2},
 | 
|---|
| 113 |     {-1, 0, 4, 0, 2},
 | 
|---|
| 114 |     {1, -1, 0, -2, 0},
 | 
|---|
| 115 |     {2, 0, 2, -2, 1},
 | 
|---|
| 116 |     {2, 0, 2, 2, 2},
 | 
|---|
| 117 |     {1, 0, 0, 2, 1},
 | 
|---|
| 118 |     {0, 0, 4, -2, 2},
 | 
|---|
| 119 |     {3, 0, 2, -2, 2},
 | 
|---|
| 120 |     {1, 0, 2, -2, 0},
 | 
|---|
| 121 |     {0, 1, 2, 0, 1},
 | 
|---|
| 122 |     {-1, -1, 0, 2, 1},
 | 
|---|
| 123 |     {0, 0, -2, 0, 1},
 | 
|---|
| 124 |     {0, 0, 2, -1, 2},
 | 
|---|
| 125 |     {0, 1, 0, 2, 0},
 | 
|---|
| 126 |     {1, 0, -2, -2, 0},
 | 
|---|
| 127 |     {0, -1, 2, 0, 1},
 | 
|---|
| 128 |     {1, 1, 0, -2, 1},
 | 
|---|
| 129 |     {1, 0, -2, 2, 0},
 | 
|---|
| 130 |     {2, 0, 0, 2, 0},
 | 
|---|
| 131 |     {0, 0, 2, 4, 2},
 | 
|---|
| 132 |     {0, 1, 0, 1, 0}
 | 
|---|
| 133 | };
 | 
|---|
| 134 | 
 | 
|---|
| 135 | /* amplitudes which  have secular terms; in 1/NUT_SCALE arc seconds
 | 
|---|
| 136 |  * {index, constant dPSI, T/10 in dPSI, constant in dEPS, T/10 in dEPS}
 | 
|---|
| 137 |  */
 | 
|---|
| 138 | static long ampsecul[][5] = {
 | 
|---|
| 139 |     {0  ,-171996 ,-1742 ,92025 ,89},
 | 
|---|
| 140 |     {1  ,2062    ,2     ,-895  ,5},
 | 
|---|
| 141 |     {8  ,-13187  ,-16   ,5736  ,-31},
 | 
|---|
| 142 |     {9  ,1426    ,-34   ,54    ,-1},
 | 
|---|
| 143 |     {10 ,-517    ,12    ,224   ,-6},
 | 
|---|
| 144 |     {11 ,217     ,-5    ,-95   ,3},
 | 
|---|
| 145 |     {12 ,129     ,1     ,-70   ,0},
 | 
|---|
| 146 |     {15 ,17      ,-1    ,0     ,0},
 | 
|---|
| 147 |     {17 ,-16     ,1     ,7     ,0},
 | 
|---|
| 148 |     {30 ,-2274   ,-2    ,977   ,-5},
 | 
|---|
| 149 |     {31 ,712     ,1     ,-7    ,0},
 | 
|---|
| 150 |     {32 ,-386    ,-4    ,200   ,0},
 | 
|---|
| 151 |     {33 ,-301    ,0     ,129   ,-1},
 | 
|---|
| 152 |     {37 ,63      ,1     ,-33   ,0},
 | 
|---|
| 153 |     {38 ,-58     ,-1    ,32    ,0},
 | 
|---|
| 154 |     /* termination */  { -1, }
 | 
|---|
| 155 | };
 | 
|---|
| 156 | 
 | 
|---|
| 157 | /* amplitudes which only have constant terms; same unit as above
 | 
|---|
| 158 |  * {dPSI, dEPS}
 | 
|---|
| 159 |  * indexes which are already in ampsecul[][] are zeroed
 | 
|---|
| 160 |  */
 | 
|---|
| 161 | static short ampconst[NUT_SERIES][2] = {
 | 
|---|
| 162 |     {0,0},
 | 
|---|
| 163 |     {0,0},
 | 
|---|
| 164 |     {46,-24},
 | 
|---|
| 165 |     {11,0},
 | 
|---|
| 166 |     {-3,1},
 | 
|---|
| 167 |     {-3,0},
 | 
|---|
| 168 |     {-2,1},
 | 
|---|
| 169 |     {1,0},
 | 
|---|
| 170 |     {0,0},
 | 
|---|
| 171 |     {0,0},
 | 
|---|
| 172 |     {0,0},
 | 
|---|
| 173 |     {0,0},
 | 
|---|
| 174 |     {0,0},
 | 
|---|
| 175 |     {48,1},
 | 
|---|
| 176 |     {-22,0},
 | 
|---|
| 177 |     {0,0},
 | 
|---|
| 178 |     {-15,9},
 | 
|---|
| 179 |     {0,0},
 | 
|---|
| 180 |     {-12,6},
 | 
|---|
| 181 |     {-6,3},
 | 
|---|
| 182 |     {-5,3},
 | 
|---|
| 183 |     {4,-2},
 | 
|---|
| 184 |     {4,-2},
 | 
|---|
| 185 |     {-4,0},
 | 
|---|
| 186 |     {1,0},
 | 
|---|
| 187 |     {1,0},
 | 
|---|
| 188 |     {-1,0},
 | 
|---|
| 189 |     {1,0},
 | 
|---|
| 190 |     {1,0},
 | 
|---|
| 191 |     {-1,0},
 | 
|---|
| 192 |     {0,0},
 | 
|---|
| 193 |     {0,0},
 | 
|---|
| 194 |     {0,0},
 | 
|---|
| 195 |     {0,0},
 | 
|---|
| 196 |     {-158,-1},
 | 
|---|
| 197 |     {123,-53},
 | 
|---|
| 198 |     {63,-2},
 | 
|---|
| 199 |     {0,0},
 | 
|---|
| 200 |     {0,0},
 | 
|---|
| 201 |     {-59,26},
 | 
|---|
| 202 |     {-51,27},
 | 
|---|
| 203 |     {-38,16},
 | 
|---|
| 204 |     {29,-1},
 | 
|---|
| 205 |     {29,-12},
 | 
|---|
| 206 |     {-31,13},
 | 
|---|
| 207 |     {26,-1},
 | 
|---|
| 208 |     {21,-10},
 | 
|---|
| 209 |     {16,-8},
 | 
|---|
| 210 |     {-13,7},
 | 
|---|
| 211 |     {-10,5},
 | 
|---|
| 212 |     {-7,0},
 | 
|---|
| 213 |     {7,-3},
 | 
|---|
| 214 |     {-7,3},
 | 
|---|
| 215 |     {-8,3},
 | 
|---|
| 216 |     {6,0},
 | 
|---|
| 217 |     {6,-3},
 | 
|---|
| 218 |     {-6,3},
 | 
|---|
| 219 |     {-7,3},
 | 
|---|
| 220 |     {6,-3},
 | 
|---|
| 221 |     {-5,3},
 | 
|---|
| 222 |     {5,0},
 | 
|---|
| 223 |     {-5,3},
 | 
|---|
| 224 |     {-4,0},
 | 
|---|
| 225 |     {4,0},
 | 
|---|
| 226 |     {-4,0},
 | 
|---|
| 227 |     {-3,0},
 | 
|---|
| 228 |     {3,0},
 | 
|---|
| 229 |     {-3,1},
 | 
|---|
| 230 |     {-3,1},
 | 
|---|
| 231 |     {-2,1},
 | 
|---|
| 232 |     {-3,1},
 | 
|---|
| 233 |     {-3,1},
 | 
|---|
| 234 |     {2,-1},
 | 
|---|
| 235 |     {-2,1},
 | 
|---|
| 236 |     {2,-1},
 | 
|---|
| 237 |     {-2,1},
 | 
|---|
| 238 |     {2,0},
 | 
|---|
| 239 |     {2,-1},
 | 
|---|
| 240 |     {1,-1},
 | 
|---|
| 241 |     {-1,0},
 | 
|---|
| 242 |     {1,-1},
 | 
|---|
| 243 |     {-2,1},
 | 
|---|
| 244 |     {-1,0},
 | 
|---|
| 245 |     {1,-1},
 | 
|---|
| 246 |     {-1,1},
 | 
|---|
| 247 |     {-1,1},
 | 
|---|
| 248 |     {1,0},
 | 
|---|
| 249 |     {1,0},
 | 
|---|
| 250 |     {1,-1},
 | 
|---|
| 251 |     {-1,0},
 | 
|---|
| 252 |     {-1,0},
 | 
|---|
| 253 |     {1,0},
 | 
|---|
| 254 |     {1,0},
 | 
|---|
| 255 |     {-1,0},
 | 
|---|
| 256 |     {1,0},
 | 
|---|
| 257 |     {1,0},
 | 
|---|
| 258 |     {-1,0},
 | 
|---|
| 259 |     {-1,0},
 | 
|---|
| 260 |     {-1,0},
 | 
|---|
| 261 |     {-1,0},
 | 
|---|
| 262 |     {-1,0},
 | 
|---|
| 263 |     {-1,0},
 | 
|---|
| 264 |     {-1,0},
 | 
|---|
| 265 |     {1,0},
 | 
|---|
| 266 |     {-1,0},
 | 
|---|
| 267 |     {1,0}
 | 
|---|
| 268 | };
 | 
|---|
| 269 | 
 | 
|---|
| 270 | /* given the modified JD, mj, find the nutation in obliquity, *deps, and
 | 
|---|
| 271 |  * the nutation in longitude, *dpsi, each in radians.
 | 
|---|
| 272 |  */
 | 
|---|
| 273 | void
 | 
|---|
| 274 | nutation (
 | 
|---|
| 275 | double mj,
 | 
|---|
| 276 | double *deps,   /* on input:  precision parameter in arc seconds */
 | 
|---|
| 277 | double *dpsi)
 | 
|---|
| 278 | {
 | 
|---|
| 279 |         static double lastmj = -10000, lastdeps, lastdpsi;
 | 
|---|
| 280 |         double T, T2, T3, T10;                  /* jul cent since J2000 */
 | 
|---|
| 281 |         double prec;                            /* series precis in arc sec */
 | 
|---|
| 282 |         int i, isecul;                          /* index in term table */
 | 
|---|
| 283 |         static double delcache[5][2*NUT_MAXMUL+1];
 | 
|---|
| 284 |                         /* cache for multiples of delaunay args
 | 
|---|
| 285 |                          * [M',M,F,D,Om][-min*x, .. , 0, .., max*x]
 | 
|---|
| 286 |                          * make static to have unfilled fields cleared on init
 | 
|---|
| 287 |                          */
 | 
|---|
| 288 | 
 | 
|---|
| 289 |         if (mj == lastmj) {
 | 
|---|
| 290 |             *deps = lastdeps;
 | 
|---|
| 291 |             *dpsi = lastdpsi;
 | 
|---|
| 292 |             return;
 | 
|---|
| 293 |         }
 | 
|---|
| 294 | 
 | 
|---|
| 295 |         prec = 0.0;
 | 
|---|
| 296 | 
 | 
|---|
| 297 | #if 0   /* this is if deps should contain a precision value */
 | 
|---|
| 298 |         prec =* deps;
 | 
|---|
| 299 |         if (prec < 0.0 || prec > 1.0)   /* accept only sane value */
 | 
|---|
| 300 |                 prec = 1.0;
 | 
|---|
| 301 | #endif
 | 
|---|
| 302 | 
 | 
|---|
| 303 |         /* augment for abundance of small terms */
 | 
|---|
| 304 |         prec *= NUT_SCALE/10;
 | 
|---|
| 305 | 
 | 
|---|
| 306 |         T = (mj - J2000)/36525.;
 | 
|---|
| 307 |         T2 = T * T;
 | 
|---|
| 308 |         T3 = T2 * T;
 | 
|---|
| 309 |         T10 = T/10.;
 | 
|---|
| 310 | 
 | 
|---|
| 311 |         /* calculate delaunay args and place in cache */
 | 
|---|
| 312 |         for (i = 0; i < 5; ++i) {
 | 
|---|
| 313 |             double x;
 | 
|---|
| 314 |             short j;
 | 
|---|
| 315 | 
 | 
|---|
| 316 |             x = delaunay[i][0] +
 | 
|---|
| 317 |                 delaunay[i][1] * T +
 | 
|---|
| 318 |                 delaunay[i][2] * T2 +
 | 
|---|
| 319 |                 delaunay[i][3] * T3;
 | 
|---|
| 320 | 
 | 
|---|
| 321 |             /* convert to radians */
 | 
|---|
| 322 |             x /= SECPERCIRC;
 | 
|---|
| 323 |             x -= floor(x);
 | 
|---|
| 324 |             x *= 2.*PI;
 | 
|---|
| 325 | 
 | 
|---|
| 326 |             /* fill cache table */
 | 
|---|
| 327 |             for (j = 0; j <= 2*NUT_MAXMUL; ++j)
 | 
|---|
| 328 |                 delcache[i][j] = (j - NUT_MAXMUL) * x;
 | 
|---|
| 329 |         }
 | 
|---|
| 330 | 
 | 
|---|
| 331 |         /* find dpsi and deps */
 | 
|---|
| 332 |         lastdpsi = lastdeps = 0.;
 | 
|---|
| 333 |         for (i = isecul = 0; i < NUT_SERIES ; ++i) {
 | 
|---|
| 334 |             double arg = 0., ampsin, ampcos;
 | 
|---|
| 335 |             short j;
 | 
|---|
| 336 | 
 | 
|---|
| 337 |             if (ampconst[i][0] || ampconst[i][1]) {
 | 
|---|
| 338 |                 /* take non-secular terms from simple array */
 | 
|---|
| 339 |                 ampsin = ampconst[i][0];
 | 
|---|
| 340 |                 ampcos = ampconst[i][1];
 | 
|---|
| 341 |             } else {
 | 
|---|
| 342 |                 /* secular terms from different array */
 | 
|---|
| 343 |                 ampsin = ampsecul[isecul][1] + ampsecul[isecul][2] * T10;
 | 
|---|
| 344 |                 ampcos = ampsecul[isecul][3] + ampsecul[isecul][4] * T10;
 | 
|---|
| 345 |                 ++isecul;
 | 
|---|
| 346 |             }
 | 
|---|
| 347 | 
 | 
|---|
| 348 |             for (j = 0; j < 5; ++j)
 | 
|---|
| 349 |                 arg += delcache[j][NUT_MAXMUL + multarg[i][j]];
 | 
|---|
| 350 | 
 | 
|---|
| 351 |             if (fabs(ampsin) >= prec)
 | 
|---|
| 352 |                 lastdpsi += ampsin * sin(arg);
 | 
|---|
| 353 | 
 | 
|---|
| 354 |             if (fabs(ampcos) >= prec)
 | 
|---|
| 355 |                 lastdeps += ampcos * cos(arg);
 | 
|---|
| 356 | 
 | 
|---|
| 357 |         }
 | 
|---|
| 358 | 
 | 
|---|
| 359 |         /* convert to radians.
 | 
|---|
| 360 |          */
 | 
|---|
| 361 |         lastdpsi = degrad(lastdpsi/3600./NUT_SCALE);
 | 
|---|
| 362 |         lastdeps = degrad(lastdeps/3600./NUT_SCALE);
 | 
|---|
| 363 | 
 | 
|---|
| 364 |         lastmj = mj;
 | 
|---|
| 365 |         *deps = lastdeps;
 | 
|---|
| 366 |         *dpsi = lastdpsi;
 | 
|---|
| 367 | }
 | 
|---|
| 368 | 
 | 
|---|
| 369 | /* given the modified JD, mj, correct, IN PLACE, the right ascension *ra
 | 
|---|
| 370 |  * and declination *dec (both in radians) for nutation.
 | 
|---|
| 371 |  */
 | 
|---|
| 372 | void
 | 
|---|
| 373 | nut_eq (double mj, double *ra, double *dec)
 | 
|---|
| 374 | {
 | 
|---|
| 375 |         static double lastmj = -10000;
 | 
|---|
| 376 |         static double a[3][3];          /* rotation matrix */
 | 
|---|
| 377 |         double xold, yold, zold, x, y, z;
 | 
|---|
| 378 | 
 | 
|---|
| 379 |         if (mj != lastmj) {
 | 
|---|
| 380 |             double epsilon, dpsi, deps;
 | 
|---|
| 381 |             double se, ce, sp, cp, sede, cede;
 | 
|---|
| 382 | 
 | 
|---|
| 383 |             obliquity(mj, &epsilon);
 | 
|---|
| 384 |             nutation(mj, &deps, &dpsi);
 | 
|---|
| 385 | 
 | 
|---|
| 386 |             /* the rotation matrix a applies the nutation correction to
 | 
|---|
| 387 |              * a vector of equatoreal coordinates Xeq to Xeq' by 3 subsequent
 | 
|---|
| 388 |              * rotations:  R1 - from equatoreal to ecliptic system by
 | 
|---|
| 389 |              * rotation of angle epsilon about x, R2 - rotate ecliptic
 | 
|---|
| 390 |              * system by -dpsi about its z, R3 - from ecliptic to equatoreal
 | 
|---|
| 391 |              * by rotation of angle -(epsilon + deps)
 | 
|---|
| 392 |              *
 | 
|---|
| 393 |              *  Xeq' = A * Xeq = R3 * R2 * R1 * Xeq
 | 
|---|
| 394 |              * 
 | 
|---|
| 395 |              *          [ 1       0          0    ]
 | 
|---|
| 396 |              * R1 =     [ 0   cos(eps)   sin(eps) ]
 | 
|---|
| 397 |              *          [ 0  - sin(eps)  cos(eps) ]
 | 
|---|
| 398 |              * 
 | 
|---|
| 399 |              *          [ cos(dpsi)  - sin(dpsi)  0 ]
 | 
|---|
| 400 |              * R2 =     [ sin(dpsi)   cos(dpsi)   0 ]
 | 
|---|
| 401 |              *          [      0           0      1 ]
 | 
|---|
| 402 |              * 
 | 
|---|
| 403 |              *          [ 1         0                 0         ]
 | 
|---|
| 404 |              * R3 =     [ 0  cos(eps + deps)  - sin(eps + deps) ]
 | 
|---|
| 405 |              *          [ 0  sin(eps + deps)   cos(eps + deps)  ]
 | 
|---|
| 406 |              * 
 | 
|---|
| 407 |              * for efficiency, here is a explicitely:
 | 
|---|
| 408 |              */
 | 
|---|
| 409 |             
 | 
|---|
| 410 |             se = sin(epsilon);
 | 
|---|
| 411 |             ce = cos(epsilon);
 | 
|---|
| 412 |             sp = sin(dpsi);
 | 
|---|
| 413 |             cp = cos(dpsi);
 | 
|---|
| 414 |             sede = sin(epsilon + deps);
 | 
|---|
| 415 |             cede = cos(epsilon + deps);
 | 
|---|
| 416 | 
 | 
|---|
| 417 |             a[0][0] = cp;
 | 
|---|
| 418 |             a[0][1] = -sp*ce;
 | 
|---|
| 419 |             a[0][2] = -sp*se;
 | 
|---|
| 420 | 
 | 
|---|
| 421 |             a[1][0] = cede*sp;
 | 
|---|
| 422 |             a[1][1] = cede*cp*ce+sede*se;
 | 
|---|
| 423 |             a[1][2] = cede*cp*se-sede*ce;
 | 
|---|
| 424 | 
 | 
|---|
| 425 |             a[2][0] = sede*sp;
 | 
|---|
| 426 |             a[2][1] = sede*cp*ce-cede*se;
 | 
|---|
| 427 |             a[2][2] = sede*cp*se+cede*ce;
 | 
|---|
| 428 | 
 | 
|---|
| 429 |             lastmj = mj;
 | 
|---|
| 430 |         }
 | 
|---|
| 431 | 
 | 
|---|
| 432 |         sphcart(*ra, *dec, 1.0, &xold, &yold, &zold);
 | 
|---|
| 433 |         x = a[0][0] * xold + a[0][1] * yold + a[0][2] * zold;
 | 
|---|
| 434 |         y = a[1][0] * xold + a[1][1] * yold + a[1][2] * zold;
 | 
|---|
| 435 |         z = a[2][0] * xold + a[2][1] * yold + a[2][2] * zold;
 | 
|---|
| 436 |         cartsph(x, y, z, ra, dec, &zold);       /* radius should be 1.0 */
 | 
|---|
| 437 |         if (*ra < 0.) *ra += 2.*PI;             /* make positive for display */
 | 
|---|
| 438 | }
 | 
|---|
| 439 | 
 | 
|---|
| 440 | /* For RCS Only -- Do Not Edit */
 | 
|---|
| 441 | static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2005-01-17 10:13:06 $ $Revision: 1.4 $ $Name: not supported by cvs2svn $"};
 | 
|---|