| [1457] | 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 |  */
 | 
|---|
| [2551] | 5 | #include <stdio.h>
 | 
|---|
| [1457] | 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 | 
 | 
|---|
| [2551] | 270 | /* given the modified JD, mj, find the nutation in obliquity, *deps, and
 | 
|---|
| [1457] | 271 |  * the nutation in longitude, *dpsi, each in radians.
 | 
|---|
 | 272 |  */
 | 
|---|
 | 273 | void
 | 
|---|
| [2551] | 274 | nutation (
 | 
|---|
 | 275 | double mj,
 | 
|---|
 | 276 | double *deps,   /* on input:  precision parameter in arc seconds */
 | 
|---|
 | 277 | double *dpsi)
 | 
|---|
| [1457] | 278 | {
 | 
|---|
| [2551] | 279 |         static double lastmj = -10000, lastdeps, lastdpsi;
 | 
|---|
| [1457] | 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 | 
 | 
|---|
| [2551] | 289 |         if (mj == lastmj) {
 | 
|---|
| [1457] | 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 | 
 | 
|---|
| [2551] | 306 |         T = (mj - J2000)/36525.;
 | 
|---|
| [1457] | 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 | 
 | 
|---|
| [2551] | 364 |         lastmj = mj;
 | 
|---|
| [1457] | 365 |         *deps = lastdeps;
 | 
|---|
 | 366 |         *dpsi = lastdpsi;
 | 
|---|
 | 367 | }
 | 
|---|
 | 368 | 
 | 
|---|
| [2551] | 369 | /* given the modified JD, mj, correct, IN PLACE, the right ascension *ra
 | 
|---|
| [1457] | 370 |  * and declination *dec (both in radians) for nutation.
 | 
|---|
 | 371 |  */
 | 
|---|
 | 372 | void
 | 
|---|
| [2551] | 373 | nut_eq (double mj, double *ra, double *dec)
 | 
|---|
| [1457] | 374 | {
 | 
|---|
| [2551] | 375 |         static double lastmj = -10000;
 | 
|---|
| [1457] | 376 |         static double a[3][3];          /* rotation matrix */
 | 
|---|
 | 377 |         double xold, yold, zold, x, y, z;
 | 
|---|
 | 378 | 
 | 
|---|
| [2551] | 379 |         if (mj != lastmj) {
 | 
|---|
| [1457] | 380 |             double epsilon, dpsi, deps;
 | 
|---|
 | 381 |             double se, ce, sp, cp, sede, cede;
 | 
|---|
 | 382 | 
 | 
|---|
| [2551] | 383 |             obliquity(mj, &epsilon);
 | 
|---|
 | 384 |             nutation(mj, &deps, &dpsi);
 | 
|---|
| [1457] | 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 | 
 | 
|---|
| [2551] | 429 |             lastmj = mj;
 | 
|---|
| [1457] | 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 */
 | 
|---|
| [3111] | 441 | static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2006-11-22 13:53:30 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};
 | 
|---|