wmmoonclock/src/MoonRise.c

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