wmsun/SunRise.c

   1 #include <stdio.h>
   2 #include <math.h>
   3
   4 #define DegPerRad       57.29577951308232087680
   5 #define RadPerDeg        0.01745329251994329576
   6
   7 extern  double  Glon, SinGlat, CosGlat, TimeZone;
   8
   9 double    cosEPS = 0.91748;
  10 double    sinEPS = 0.39778;
  11 double    P2  = 6.283185307;
  12
  13 Interp(double ym, double y0, double yp, double *xe, double *ye, double *z1, double *z2, int *nz){
  14
  15     double      a, b, c, d, dx;
  16
  17     *nz = 0;
  18     a = 0.5*(ym+yp)-y0;
  19     b = 0.5*(yp-ym);
  20     c = y0;
  21     *xe = -b/(2.0*a);
  22     *ye = (a*(*xe) + b) * (*xe) + c;
  23     d = b*b - 4.0*a*c;
  24
  25     if (d >= 0){
  26         dx = 0.5*sqrt(d)/fabs(a);
  27         *z1 = *xe - dx;
  28         *z2 = *xe+dx;
  29         if (fabs(*z1) <= 1.0) *nz += 1;
  30         if (fabs(*z2) <= 1.0) *nz += 1;
  31         if (*z1 < -1.0) *z1 = *z2;
  32     }
  33
  34     return(0);
  35
  36
  37 }
  38
  39 SunRise(int year, int month, int day, double LocalHour, double *UTRise, double *UTSet){
  40
  41     double      UT, ym, y0, yp, SinH0;
  42     double      xe, ye, z1, z2, SinH(), hour24();
  43     int         Rise, Set, nz;
  44
  45     SinH0 = sin( -50.0/60.0 * RadPerDeg );
  46
  47
  48     UT = 1.0+TimeZone;
  49     *UTRise = -999.0;
  50     *UTSet = -999.0;
  51     Rise = Set = 0;
  52     ym = SinH(year, month, day, UT-1.0) - SinH0;
  53
  54     while ( (UT <= 24.0+TimeZone) ) {
  55
  56         y0 = SinH(year, month, day, UT) - SinH0;
  57         yp = SinH(year, month, day, UT+1.0) - SinH0;
  58
  59         Interp(ym, y0, yp, &xe, &ye, &z1, &z2, &nz);
  60
  61         switch(nz){
  62
  63                 case 0:
  64                         break;
  65                 case 1:
  66                         if (ym < 0.0){
  67                             *UTRise = UT + z1;
  68                             Rise = 1;
  69                         } else {
  70                             *UTSet = UT + z1;
  71                             Set = 1;
  72                         }
  73                         break;
  74                 case 2:
  75                         if (ye < 0.0){
  76                             *UTRise = UT + z2;
  77                             *UTSet = UT + z1;
  78                         } else {
  79                             *UTRise = UT + z1;
  80                             *UTSet = UT + z2;
  81                         }
  82                         Rise = 1;
  83                         Set = 1;
  84                         break;
  85         }
  86         ym = yp;
  87         UT += 2.0;
  88
  89     }
  90
  91     if (Rise){
  92         *UTRise -= TimeZone;
  93         *UTRise = hour24(*UTRise);
  94     } else {
  95         *UTRise = -999.0;
  96     }
  97
  98     if (Set){
  99         *UTSet -= TimeZone;
 100         *UTSet = hour24(*UTSet);
 101     } else {
 102         *UTSet = -999.0;
 103     }
 104
 105 }
 106
 107
 108 UTTohhmm(double UT, int *h, int *m){
 109
 110
 111     if (UT < 0.0) {
 112         *h = -1.0;
 113         *m = -1.0;
 114     } else {
 115         *h = (int)UT;
 116         *m = (int)((UT-(double)(*h))*60.0+0.5);
 117     }
 118
 119 }
 120
 121 double SinH(int year, int month, int day, double UT){
 122
 123     double      TU, frac(), jd();
 124     double      RA_Sun, DEC_Sun, gmst, lmst, Tau;
 125     double      M, DL, L, SL, X, Y, Z, RHO;
 126
 127
 128     TU = (jd(year, month, day, UT+62.0/3600.0) - 2451545.0)/36525.0;
 129
 130     M = P2*frac(0.993133 + 99.997361*TU);
 131     DL = 6893.0*sin(M) + 72.0*sin(2.0*M);
 132     L = P2*frac(0.7859453 + M/P2 + (6191.2*TU+DL)/1296e3);
 133     SL = sin(L);
 134     X = cos(L); Y = cosEPS*SL; Z = sinEPS*SL; RHO = sqrt(1.0-Z*Z);
 135     DEC_Sun = atan2(Z, RHO);
 136     RA_Sun = (48.0/P2)*atan(Y/(X+RHO));
 137     if (RA_Sun < 0) RA_Sun += 24.0;
 138
 139     RA_Sun = RA_Sun*15.0*RadPerDeg;
 140
 141     /*
 142      *  Compute Greenwich Mean Sidereal Time (gmst)
 143      */
 144     UT = 24.0*frac( UT/24.0 );
 145
 146     gmst = 6.697374558 + 1.0*UT + (8640184.812866+(0.093104-6.2e-6*TU)*TU)*TU/3600.0;
 147     lmst = 24.0*frac( (gmst-Glon/15.0) / 24.0 );
 148
 149     Tau = 15.0*lmst*RadPerDeg - RA_Sun;
 150     return( SinGlat*sin(DEC_Sun) + CosGlat*cos(DEC_Sun)*cos(Tau) );
 151
 152
 153 }
 154
 155
 156 /*
 157  *  Compute the Julian Day number for the given date.
 158  *  Julian Date is the number of days since noon of Jan 1 4713 B.C.
 159  */
 160 double jd(ny, nm, nd, UT)
 161 int ny, nm, nd;
 162 double UT;
 163 {
 164         double A, B, C, D, JD, day;
 165
 166         day = nd + UT/24.0;
 167
 168
 169         if ((nm == 1) || (nm == 2)){
 170                 ny = ny - 1;
 171                 nm = nm + 12;
 172         }
 173
 174         if (((double)ny+nm/12.0+day/365.25)>=(1582.0+10.0/12.0+15.0/365.25)){
 175                         A = ((int)(ny / 100.0));
 176                         B = 2.0 - A + (int)(A/4.0);
 177         }
 178         else{
 179                         B = 0.0;
 180         }
 181
 182         if (ny < 0.0){
 183                 C = (int)((365.25*(double)ny) - 0.75);
 184         }
 185         else{
 186                 C = (int)(365.25*(double)ny);
 187         }
 188
 189         D = (int)(30.6001*(double)(nm+1));
 190
 191
 192         JD = B + C + D + day + 1720994.5;
 193         return(JD);
 194
 195 }
 196
 197 double hour24(hour)
 198 double hour;
 199 {
 200         int n;
 201
 202         if (hour < 0.0){
 203                 n = (int)(hour/24.0) - 1;
 204                 return(hour-n*24.0);
 205         }
 206         else if (hour > 24.0){
 207                 n = (int)(hour/24.0);
 208                 return(hour-n*24.0);
 209         }
 210         else{
 211                 return(hour);
 212         }
 213 }
 214
 215 double frac(double x){
 216
 217     x -= (int)x;
 218     return( (x<0) ? x+1.0 : x );
 219
 220 }
 221