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