moon.c (8162B)
1 #include "astro.h" 2 3 double k1, k2, k3, k4; 4 double mnom, msun, noded, dmoon; 5 6 void 7 moon(void) 8 { 9 Moontab *mp; 10 double dlong, lsun, psun; 11 double eccm, eccs, chp, cpe; 12 double v0, t0, m0, j0; 13 double arg1, arg2, arg3, arg4, arg5, arg6, arg7; 14 double arg8, arg9, arg10; 15 double dgamma, k5, k6; 16 double lterms, sterms, cterms, nterms, pterms, spterms; 17 double gamma1, gamma2, gamma3, arglat; 18 double xmp, ymp, zmp; 19 double obl2; 20 21 /* 22 * the fundamental elements - all referred to the epoch of 23 * Jan 0.5, 1900 and to the mean equinox of date. 24 */ 25 26 dlong = 270.434164 + 13.1763965268*eday - .001133*capt2 27 + 2.e-6*capt3; 28 argp = 334.329556 + .1114040803*eday - .010325*capt2 29 - 12.e-6*capt3; 30 node = 259.183275 - .0529539222*eday + .002078*capt2 31 + 2.e-6*capt3; 32 lsun = 279.696678 + .9856473354*eday + .000303*capt2; 33 psun = 281.220833 + .0000470684*eday + .000453*capt2 34 + 3.e-6*capt3; 35 36 dlong = fmod(dlong, 360.); 37 argp = fmod(argp, 360.); 38 node = fmod(node, 360.); 39 lsun = fmod(lsun, 360.); 40 psun = fmod(psun, 360.); 41 42 eccm = 22639.550; 43 eccs = .01675104 - .00004180*capt; 44 incl = 18461.400; 45 cpe = 124.986; 46 chp = 3422.451; 47 48 /* 49 * some subsidiary elements - they are all longitudes 50 * and they are referred to the epoch 1/0.5 1900 and 51 * to the fixed mean equinox of 1850.0. 52 */ 53 54 v0 = 342.069128 + 1.6021304820*eday; 55 t0 = 98.998753 + 0.9856091138*eday; 56 m0 = 293.049675 + 0.5240329445*eday; 57 j0 = 237.352319 + 0.0830912295*eday; 58 59 /* 60 * the following are periodic corrections to the 61 * fundamental elements and constants. 62 * arg3 is the "Great Venus Inequality". 63 */ 64 65 arg1 = 41.1 + 20.2*(capt+.5); 66 arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5); 67 arg3 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5); 68 arg4 = node; 69 arg5 = node + 276.2 - 2.3*(capt+.5); 70 arg6 = 313.9 + 13.*t0 - 8.*v0; 71 arg7 = dlong - argp + 112.0 + 29.*t0 - 26.*v0; 72 arg8 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0; 73 arg9 = node + 290.1 - 0.9*(capt+.5); 74 arg10 = 115. + 38.5*(capt+.5); 75 arg1 *= radian; 76 arg2 *= radian; 77 arg3 *= radian; 78 arg4 *= radian; 79 arg5 *= radian; 80 arg6 *= radian; 81 arg7 *= radian; 82 arg8 *= radian; 83 arg9 *= radian; 84 arg10 *= radian; 85 86 dlong += 87 (0.84 *sin(arg1) 88 + 0.31 *sin(arg2) 89 + 14.27 *sin(arg3) 90 + 7.261*sin(arg4) 91 + 0.282*sin(arg5) 92 + 0.237*sin(arg6) 93 + 0.108*sin(arg7) 94 + 0.126*sin(arg8))/3600.; 95 96 argp += 97 (- 2.10 *sin(arg1) 98 - 0.118*sin(arg3) 99 - 2.076*sin(arg4) 100 - 0.840*sin(arg5) 101 - 0.593*sin(arg6))/3600.; 102 103 node += 104 (0.63*sin(arg1) 105 + 0.17*sin(arg3) 106 + 95.96*sin(arg4) 107 + 15.58*sin(arg5) 108 + 1.86*sin(arg9))/3600.; 109 110 t0 += 111 (- 6.40*sin(arg1) 112 - 1.89*sin(arg6))/3600.; 113 114 psun += 115 (6.40*sin(arg1) 116 + 1.89*sin(arg6))/3600.; 117 118 dgamma = - 4.318*cos(arg4) 119 - 0.698*cos(arg5) 120 - 0.083*cos(arg9); 121 122 j0 += 123 0.33*sin(arg10); 124 125 /* 126 * the following factors account for the fact that the 127 * eccentricity, solar eccentricity, inclination and 128 * parallax used by Brown to make up his coefficients 129 * are both wrong and out of date. Brown did the same 130 * thing in a different way. 131 */ 132 133 k1 = eccm/22639.500; 134 k2 = eccs/.01675104; 135 k3 = 1. + 2.708e-6 + .000108008*dgamma; 136 k4 = cpe/125.154; 137 k5 = chp/3422.700; 138 139 /* 140 * the principal arguments that are used to compute 141 * perturbations are the following differences of the 142 * fundamental elements. 143 */ 144 145 mnom = dlong - argp; 146 msun = lsun - psun; 147 noded = dlong - node; 148 dmoon = dlong - lsun; 149 150 /* 151 * solar terms in longitude 152 */ 153 154 lterms = 0.0; 155 mp = moontab; 156 for(;;) { 157 if(mp->f == 0.0) 158 break; 159 lterms += sinx(mp->f, 160 mp->c[0], mp->c[1], 161 mp->c[2], mp->c[3], 0.0); 162 mp++; 163 } 164 mp++; 165 166 /* 167 * planetary terms in longitude 168 */ 169 170 lterms += sinx(0.822, 0,0,0,0, t0-v0); 171 lterms += sinx(0.307, 0,0,0,0, 2.*t0-2.*v0+179.8); 172 lterms += sinx(0.348, 0,0,0,0, 3.*t0-2.*v0+272.9); 173 lterms += sinx(0.176, 0,0,0,0, 4.*t0-3.*v0+271.7); 174 lterms += sinx(0.092, 0,0,0,0, 5.*t0-3.*v0+199.); 175 lterms += sinx(0.129, 1,0,0,0, -t0+v0+180.); 176 lterms += sinx(0.152, 1,0,0,0, t0-v0); 177 lterms += sinx(0.127, 1,0,0,0, 3.*t0-3.*v0+180.); 178 lterms += sinx(0.099, 0,0,0,2, t0-v0); 179 lterms += sinx(0.136, 0,0,0,2, 2.*t0-2.*v0+179.5); 180 lterms += sinx(0.083, -1,0,0,2, -4.*t0+4.*v0+180.); 181 lterms += sinx(0.662, -1,0,0,2, -3.*t0+3.*v0+180.0); 182 lterms += sinx(0.137, -1,0,0,2, -2.*t0+2.*v0); 183 lterms += sinx(0.133, -1,0,0,2, t0-v0); 184 lterms += sinx(0.157, -1,0,0,2, 2.*t0-2.*v0+179.6); 185 lterms += sinx(0.079, -1,0,0,2, -8.*t0+6.*v0+162.6); 186 lterms += sinx(0.073, 2,0,0,-2, 3.*t0-3.*v0+180.); 187 lterms += sinx(0.643, 0,0,0,0, -t0+j0+178.8); 188 lterms += sinx(0.187, 0,0,0,0, -2.*t0+2.*j0+359.6); 189 lterms += sinx(0.087, 0,0,0,0, j0+289.9); 190 lterms += sinx(0.165, 0,0,0,0, -t0+2.*j0+241.5); 191 lterms += sinx(0.144, 1,0,0,0, t0-j0+1.0); 192 lterms += sinx(0.158, 1,0,0,0, -t0+j0+179.0); 193 lterms += sinx(0.190, 1,0,0,0, -2.*t0+2.*j0+180.0); 194 lterms += sinx(0.096, 1,0,0,0, -2.*t0+3.*j0+352.5); 195 lterms += sinx(0.070, 0,0,0,2, 2.*t0-2.*j0+180.); 196 lterms += sinx(0.167, 0,0,0,2, -t0+j0+178.5); 197 lterms += sinx(0.085, 0,0,0,2, -2.*t0+2.*j0+359.2); 198 lterms += sinx(1.137, -1,0,0,2, 2.*t0-2.*j0+180.3); 199 lterms += sinx(0.211, -1,0,0,2, -t0+j0+178.4); 200 lterms += sinx(0.089, -1,0,0,2, -2.*t0+2.*j0+359.2); 201 lterms += sinx(0.436, -1,0,0,2, 2.*t0-3.*j0+7.5); 202 lterms += sinx(0.240, 2,0,0,-2, -2.*t0+2.*j0+179.9); 203 lterms += sinx(0.284, 2,0,0,-2, -2.*t0+3.*j0+172.5); 204 lterms += sinx(0.195, 0,0,0,0, -2.*t0+2.*m0+180.2); 205 lterms += sinx(0.327, 0,0,0,0, -t0+2.*m0+224.4); 206 lterms += sinx(0.093, 0,0,0,0, -2.*t0+4.*m0+244.8); 207 lterms += sinx(0.073, 1,0,0,0, -t0+2.*m0+223.3); 208 lterms += sinx(0.074, 1,0,0,0, t0-2.*m0+306.3); 209 lterms += sinx(0.189, 0,0,0,0, node+180.); 210 211 /* 212 * solar terms in latitude 213 */ 214 215 sterms = 0; 216 for(;;) { 217 if(mp->f == 0) 218 break; 219 sterms += sinx(mp->f, 220 mp->c[0], mp->c[1], 221 mp->c[2], mp->c[3], 0); 222 mp++; 223 } 224 mp++; 225 226 cterms = 0; 227 for(;;) { 228 if(mp->f == 0) 229 break; 230 cterms += cosx(mp->f, 231 mp->c[0], mp->c[1], 232 mp->c[2], mp->c[3], 0); 233 mp++; 234 } 235 mp++; 236 237 nterms = 0; 238 for(;;) { 239 if(mp->f == 0) 240 break; 241 nterms += sinx(mp->f, 242 mp->c[0], mp->c[1], 243 mp->c[2], mp->c[3], 0); 244 mp++; 245 } 246 mp++; 247 248 /* 249 * planetary terms in latitude 250 */ 251 252 pterms = 253 sinx(0.215, 0,0,0,0, dlong); 254 255 /* 256 * solar terms in parallax 257 */ 258 259 spterms = 3422.700; 260 for(;;) { 261 if(mp->f == 0) 262 break; 263 spterms += cosx(mp->f, 264 mp->c[0], mp->c[1], 265 mp->c[2], mp->c[3], 0); 266 mp++; 267 } 268 269 /* 270 * planetary terms in parallax 271 */ 272 273 //spterms = spterms; 274 275 /* 276 * computation of longitude 277 */ 278 279 lambda = (dlong + lterms/3600.)*radian; 280 281 /* 282 * computation of latitude 283 */ 284 285 arglat = (noded + sterms/3600.)*radian; 286 gamma1 = 18519.700 * k3; 287 gamma2 = -6.241 * k3*k3*k3; 288 gamma3 = 0.004 * k3*k3*k3*k3*k3; 289 290 k6 = (gamma1 + cterms) / gamma1; 291 292 beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat) 293 + gamma3*sin(5.*arglat) + nterms) 294 + pterms; 295 if(flags['o']) 296 beta -= 0.6; 297 beta *= radsec; 298 299 /* 300 * computation of parallax 301 */ 302 303 spterms = k5 * spterms *radsec; 304 hp = spterms + (spterms*spterms*spterms)/6.; 305 306 rad = hp/radsec; 307 rp = 1.; 308 semi = .0799 + .272453*(hp/radsec); 309 if(dmoon < 0.) 310 dmoon += 360.; 311 mag = dmoon/360.; 312 313 /* 314 * change to equatorial coordinates 315 */ 316 317 lambda += phi; 318 obl2 = obliq + eps; 319 xmp = rp*cos(lambda)*cos(beta); 320 ymp = rp*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta)); 321 zmp = rp*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta)); 322 323 alpha = atan2(ymp, xmp); 324 delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp)); 325 meday = eday; 326 mhp = hp; 327 328 geo(); 329 } 330 331 double 332 sinx(double coef, int i, int j, int k, int m, double angle) 333 { 334 double x; 335 336 x = i*mnom + j*msun + k*noded + m*dmoon + angle; 337 x = coef*sin(x*radian); 338 if(i < 0) 339 i = -i; 340 for(; i>0; i--) 341 x *= k1; 342 if(j < 0) 343 j = -j; 344 for(; j>0; j--) 345 x *= k2; 346 if(k < 0) 347 k = -k; 348 for(; k>0; k--) 349 x *= k3; 350 if(m & 1) 351 x *= k4; 352 353 return x; 354 } 355 356 double 357 cosx(double coef, int i, int j, int k, int m, double angle) 358 { 359 double x; 360 361 x = i*mnom + j*msun + k*noded + m*dmoon + angle; 362 x = coef*cos(x*radian); 363 if(i < 0) 364 i = -i; 365 for(; i>0; i--) 366 x *= k1; 367 if(j < 0) 368 j = -j; 369 for(; j>0; j--) 370 x *= k2; 371 if(k < 0) 372 k = -k; 373 for(; k>0; k--) 374 x *= k3; 375 if(m & 1) 376 x *= k4; 377 378 return x; 379 }