sys/src/cmd/astro/sun.c

   1 #include "astro.h"
   2
   3 void
   4 sun(void)
   5 {
   6         double mven, merth, mmars, mjup, msat;
   7         double dmoon, mmoon, gmoon;
   8         double pturbb, pturbl, pturbr, lograd;
   9
  10         ecc = .01675104 - 4.180e-5 * capt - 1.26e-7*capt2;
  11         incl = 0;
  12         node = 0;
  13         argp = 281.220833 + .0000470684*eday + .000453*capt2
  14                  + .000003*capt3;
  15         mrad = 1;
  16         anom = 358.475845 + .9856002670*eday - .000150*capt2
  17                  - .000003*capt3;
  18         motion = .9856473354;
  19
  20         dmoon = 350.737681+12.1907491914*eday-.001436*capt2;
  21         gmoon = 11.250889 + 13.2293504490*eday - .003212*capt2;
  22         mmoon = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2;
  23         mven  = 212.448 + 1.602121635*eday;
  24         merth = 358.476 + 0.985600267*eday;
  25         mmars = 319.590 + .524024095*eday;
  26         mjup = 225.269 + .083082362*eday;
  27         msat  = 175.593 + .033450794*eday;
  28
  29         dmoon = fmod(dmoon, 360.)*radian;
  30         gmoon = fmod(gmoon, 360.)*radian;
  31         mmoon = fmod(mmoon, 360.)*radian;
  32         mven  *= radian;
  33         merth *= radian;
  34         mmars *= radian;
  35         mjup *= radian;
  36         msat  *= radian;
  37
  38         icosadd(sunfp, suncp);
  39         anom += cosadd(4, mmars, merth, mven, mjup)/3600.;
  40         anom += sinadd(5, mmars, merth, mven, mjup, .07884*capt)/3600.;
  41
  42         incl *= radian;
  43         node *= radian;
  44         argp *= radian;
  45         anom = fmod(anom, 360.)*radian;
  46
  47 /*
  48  *      computation of elliptic orbit
  49  */
  50
  51         lambda = anom + argp;
  52
  53         pturbl = (6910.057 - 17.240*capt - 0.052*capt2)*sin(anom)
  54                  + (72.338 - 0.361*capt) * sin(2.*anom)
  55                  + (1.054 - 0.001*capt) * sin(3.*anom)
  56                  + 0.018 * sin(4.*anom);
  57
  58         lambda += pturbl*radsec;
  59
  60         beta = 0.;
  61
  62         lograd = (30.57e-6 - 0.15e-6*capt)
  63                  - (7274.12e-6 - 18.14e-6*capt - 0.05e-6*capt2)*cos(anom)
  64                  - (91.38e-6 - 0.46e-6*capt) * cos(2.*anom)
  65                  - (1.45e-6 - 0.01e-6*capt) * cos(3.*anom)
  66                  - 0.02e-6 * cos(4.*anom);
  67
  68         pturbl = cosadd(5, mmars, merth, mven, mjup, msat);
  69         pturbl += sinadd(3, dmoon, mmoon, merth) + .9;
  70         pturbl *= radsec;
  71
  72         pturbb =  cosadd(3, merth, mven, mjup);
  73         pturbb += sinadd(3, gmoon, mmoon, dmoon);
  74         pturbb *= radsec;
  75
  76         pturbr =  cosadd(5, mmars, merth, mven, mjup, msat);
  77         pturbr += cosadd(3, dmoon, mmoon, merth);
  78
  79         lambda += pturbl;
  80         if(lambda > pipi)
  81                 lambda -= pipi;
  82
  83         beta += pturbb;
  84
  85         lograd = (lograd+pturbr) * 2.30258509;
  86         rad = 1 + lograd * (1 + lograd * (.5 + lograd/6));
  87
  88         motion *= radian*mrad*mrad/(rad*rad);
  89
  90         semi = 961.182;
  91         if(flags['o'])
  92                 semi = 959.63;
  93         mag = -26.5;
  94 }