wmmoonclock/src/MoonRise.c

   1 #include <math.h>
   2 #include "MoonRise.h"
   3 #include "Moon.h"
   4
   5 #define DegPerRad       57.29577951308232087680
   6 #define RadPerDeg        0.01745329251994329576
   7
   8 extern  double  Glon, SinGlat, CosGlat, TimeZone;
   9
  10
  11 void MoonRise(int year, int month, int day, double LocalHour, double *UTRise, double *UTSet){
  12
  13     double      UT, ym, y0, yp, SinH0;
  14     double      xe, ye, z1, z2, SinH(), hour24();
  15     int         Rise, Set, nz;
  16
  17     SinH0 = sin( 8.0/60.0 * RadPerDeg );
  18
  19
  20     UT = 1.0+TimeZone;
  21     *UTRise = -999.0;
  22     *UTSet = -999.0;
  23     Rise = Set = 0;
  24     ym = SinH(year, month, day, UT-1.0) - SinH0;
  25
  26     while ( (UT <= 24.0+TimeZone) ) {
  27
  28         y0 = SinH(year, month, day, UT) - SinH0;
  29         yp = SinH(year, month, day, UT+1.0) - SinH0;
  30
  31         Interp(ym, y0, yp, &xe, &ye, &z1, &z2, &nz);
  32
  33         switch(nz){
  34
  35                 case 0:
  36                         break;
  37                 case 1:
  38                         if (ym < 0.0){
  39                             *UTRise = UT + z1;
  40                             Rise = 1;
  41                         } else {
  42                             *UTSet = UT + z1;
  43                             Set = 1;
  44                         }
  45                         break;
  46                 case 2:
  47                         if (ye < 0.0){
  48                             *UTRise = UT + z2;
  49                             *UTSet = UT + z1;
  50                         } else {
  51                             *UTRise = UT + z1;
  52                             *UTSet = UT + z2;
  53                         }
  54                         Rise = 1;
  55                         Set = 1;
  56                         break;
  57         }
  58         ym = yp;
  59         UT += 2.0;
  60
  61     }
  62
  63     if (Rise){
  64         *UTRise -= TimeZone;
  65         *UTRise = hour24(*UTRise);
  66     } else {
  67         *UTRise = -999.0;
  68     }
  69
  70     if (Set){
  71         *UTSet -= TimeZone;
  72         *UTSet = hour24(*UTSet);
  73     } else {
  74         *UTSet = -999.0;
  75     }
  76
  77 }
  78
  79
  80 void UTTohhmm(double UT, int *h, int *m){
  81
  82
  83     if (UT < 0.0) {
  84         *h = -1.0;
  85         *m = -1.0;
  86     } else {
  87         *h = (int)UT;
  88         *m = (int)((UT-(double)(*h))*60.0+0.5);
  89         if (*m == 60) {
  90             /*
  91              *  If it was 23:60 this should become 24:00
  92              *  I prefer this designation to 00:00. So dont reset h to 0 when it goes above 24.
  93              */
  94             *h += 1;
  95             *m = 0;
  96         }
  97     }
  98
  99 }
 100
 101
 102
 103
 104
 105
 106 void Interp(double ym, double y0, double yp, double *xe, double *ye, double *z1, double *z2, int *nz){
 107
 108     double      a, b, c, d, dx;
 109
 110     *nz = 0;
 111     a = 0.5*(ym+yp)-y0;
 112     b = 0.5*(yp-ym);
 113     c = y0;
 114     *xe = -b/(2.0*a);
 115     *ye = (a*(*xe) + b) * (*xe) + c;
 116     d = b*b - 4.0*a*c;
 117
 118     if (d >= 0){
 119         dx = 0.5*sqrt(d)/fabs(a);
 120         *z1 = *xe - dx;
 121         *z2 = *xe+dx;
 122         if (fabs(*z1) <= 1.0) *nz += 1;
 123         if (fabs(*z2) <= 1.0) *nz += 1;
 124         if (*z1 < -1.0) *z1 = *z2;
 125     }
 126
 127
 128
 129 }
 130
 131
 132
 133
 134 double SinH(int year, int month, int day, double UT){
 135
 136     double      TU, frac(), jd();
 137     double      RA_Moon, DEC_Moon, gmst, lmst, Tau, Moon();
 138     double      angle2pi();
 139
 140     TU = (jd(year, month, day, UT) - 2451545.0)/36525.0;
 141
 142     /* this is more accurate, but wasteful for this -- use low res approx.
 143     TU2 = TU*TU;
 144     TU3 = TU2*TU;
 145     Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE);
 146     LambdaMoon *= RadPerDeg;
 147     BetaMoon *= RadPerDeg;
 148     epsilon = (23.43929167 - 0.013004166*TU - 1.6666667e-7*TU2 - 5.0277777778e-7*TU3)*RadPerDeg;
 149     RA_Moon  = angle2pi(atan2(sin(LambdaMoon)*cos(epsilon)-tan(BetaMoon)*sin(epsilon), cos(LambdaMoon)));
 150     DEC_Moon = asin( sin(BetaMoon)*cos(epsilon) + cos(BetaMoon)*sin(epsilon)*sin(LambdaMoon));
 151     */
 152
 153     MiniMoon(TU, &RA_Moon, &DEC_Moon);
 154     RA_Moon *= 15.0*RadPerDeg;
 155     DEC_Moon *= RadPerDeg;
 156
 157     /*
 158      *  Compute Greenwich Mean Sidereal Time (gmst)
 159      */
 160     UT = 24.0*frac( UT/24.0 );
 161     /* this is for the ephemeris meridian???
 162     gmst = 6.697374558 + 1.0027379093*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
 163     */
 164     gmst = UT + 6.697374558 + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
 165     lmst = 24.0*frac( (gmst-Glon/15.0) / 24.0 );
 166
 167     Tau = 15.0*lmst*RadPerDeg - RA_Moon;
 168     return( SinGlat*sin(DEC_Moon) + CosGlat*cos(DEC_Moon)*cos(Tau) );
 169
 170
 171 }
 172