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
Interp(double ym,double y0,double yp,double * xe,double * ye,double * z1,double * z2,int * nz)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
SunRise(int year,int month,int day,double LocalHour,double * UTRise,double * UTSet)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
SinH(int year,int month,int day,double UT)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 */
jd(ny,nm,nd,UT)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
hour24(hour)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
frac(double x)206 double frac(double x){
207
208 x -= (int)x;
209 return( (x<0) ? x+1.0 : x );
210
211 }
212
213