wmsun/SunRise.c

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