sys/src/cmd/astro/moon.c

   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 }