forked from risacher/sunwait
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sunriset.c
executable file
·288 lines (241 loc) · 12.2 KB
/
sunriset.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
/*
** sunriset.c - computes Sun rise/set times, including twilights
** Written as DAYLEN.C, 1989-08-16
** Modified to SUNRISET.C, 1992-12-01
** (c) Paul Schlyter, 1989, 1992
** Released to the public domain by Paul Schlyter, December 1992
*/
/*
**
** Who When Ver What
** IFC 04-12-2014 0.5 General fixes to "sunwait for windows" and linux port
** IFC 08-12-2014 0.6 Add timezone for output of timings
** IFC 29-04-2015 0.7 Fix for timezone (esp near date line) and more debug
** IFC 2015-05-27 0.8 Resolve 'dodgy day' and cleanup
** DRR 2022-09-16 0.91 Move to C courtesy mstilkerich
**
*/
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "sunwait.h"
#include "sunriset.h"
/************************************************************************/
/* Note: Eastern longitude positive, Western longitude negative */
/* Northern latitude positive, Southern latitude negative */
/* */
/* >>> Longitude value IS critical in this function! <<< */
/* */
/* days = Days since 2000 plus fraction to local noon */
/* altit = the altitude which the Sun should cross */
/* Set to -35/60 degrees for rise/set, -6 degrees */
/* for civil, -12 degrees for nautical and -18 */
/* degrees for astronomical twilight. */
/* upper_limb: non-zero -> upper limb, zero -> center */
/* Set to non-zero (e.g. 1) when computing rise/set */
/* times, and to zero when computing start/end of */
/* twilight. */
/* */
/* Both times are relative to the specified altitude, */
/* and thus this function can be used to comupte */
/* various twilight times, as well as rise/set times */
/* Return Codes: */
/* 0 = sun rises/sets this day. Success. */
/* Times held at *trise and *tset. */
/* +1 = Midnight Sun. Fail. */
/* Sun above the specified "horizon" all 24 hours. */
/* *trise set to time when the sun is at south, */
/* minus 12 hours while *tset is set to the south */
/* time plus 12 hours. "Day" length = 24 hours */
/* -1 = Polar Night. Fail. */
/* Sun is below the specified "horizon" all 24hours. */
/* "Day" length = 0 hours, *trise and *tset are */
/* both set to the time when the sun is at south. */
/* */
/************************************************************************/
void sunriset (const runStruct *pRun, targetStruct *pTarget)
{
double sr; /* solar distance, astronomical units */
double sra; /* sun's right ascension */
double sdec; /* sun's declination */
double sradius; /* sun's apparent radius */
double siderealTime; /* local sidereal time */
double altitude; /* sun's altitude: angle to the sun relative to the mathematical (flat-earth) horizon */
double diurnalArc = 0.0; /* the diurnal arc, hours */
double southHour = 0.0; /* Hour UTC the sun is directly south (or north for southern Hemisphere) of lat/long position */
/* compute sideral time at 00:00 UTC of target day for this longitude. */
siderealTime = revolution (GMST0(pTarget->daysSince2000) + 180.0 + pRun->longitude); // 180 = 0 hour UTC is measured 180 degrees from dateline
/* compute sun's ra + decl at this moment */
sun_RA_dec (pTarget->daysSince2000, &sra, &sdec, &sr );
/* compute time when sun is directly south - in hours UTC. "12.00" == noon. "15" == 180degrees/12hours [degrees per hour] */
southHour = 12.0 - rev180 (siderealTime - sra)/15.0;
/* compute the sun's apparent radius, degrees */
sradius = 0.2666 / sr; // Apparent angular radius of sun is 0.2666/distance in AU (deg)
/* Do correction for upper limb ('top' of sun) only, for "daylight" sunrise or set. Otherwise calculate for centre of sun */
if (pTarget->twilightAngle == TWILIGHT_ANGLE_DAYLIGHT)
altitude = pTarget->twilightAngle - sradius;
else
altitude = pTarget->twilightAngle;
/* compute the diurnal arc that the sun traverses to reach the specified altitide altit: */
double cost = (sind(altitude) - sind(pRun->latitude) * sind(sdec)) / (cosd(pRun->latitude) * cosd(sdec));
if (cost > -1.0 && cost < 1.0)
diurnalArc = 2*acosd(cost)/15.0; /* Diurnal arc, hours */
else if (cost >= 1.0)
diurnalArc = 0.0; // Polar Night
else
diurnalArc = 24.0; // Midnight Sun
if (pRun->debug == ONOFF_ON)
{ printf ("Debug: %s: Sun directly south: %f UTC, Diurnal Arc = %f hours\n", __FILE__, southHour, diurnalArc);
printf ("Debug: %s: Days since 2000: %li\n", __FILE__, pTarget->daysSince2000);
if (diurnalArc >= 24.0) printf ("Debug: %s: No rise or set: Midnight Sun\n", __FILE__);
if (diurnalArc <= 0.0) printf ("Debug: %s: No rise or set: Polar Night\n", __FILE__);
}
// Error Check - just make sure odd things don't happen (causing trouble further on)
if (diurnalArc > 24.0) diurnalArc = 24.0;
if (diurnalArc < 0.0) diurnalArc = 0.0;
/* Apply values */
pTarget->southHourUTC = southHour;
pTarget->diurnalArc = diurnalArc;
}
void sunpos (const double d, double *lon, double *r)
/******************************************************/
/* Computes the Sun's ecliptic longitude and distance */
/* at an instant given in d, number of days since */
/* 2000 Jan 0.0. The Sun's ecliptic latitude is not */
/* computed, since it's always very near 0. */
/******************************************************/
{
double M, /* Mean anomaly of the Sun */
w, /* Mean longitude of perihelion */
/* Note: Sun's mean longitude = M + w */
e, /* Eccentricity of Earth's orbit */
E, /* Eccentric anomaly */
x, y, /* x, y coordinates in orbit */
v; /* True anomaly */
/* Compute mean elements */
M = revolution (356.0470 + 0.9856002585 * d);
w = 282.9404 + 4.70935E-5 * d;
e = 0.016709 - 1.151E-9 * d;
/* Compute true longitude and radius vector */
E = M + e * RADIAN_TO_DEGREE * sind(M) * (1.0 + e * cosd(M));
x = cosd (E) - e;
y = sqrt (1.0 - e*e) * sind(E);
*r = sqrt (x*x + y*y); /* Solar distance */
v = atan2d (y, x); /* True anomaly */
*lon = revolution (v + w); /* True solar longitude, made 0..360 degrees */
}
void sun_RA_dec (const double d, double *RA, double *dec, double *r)
{
double lon, obl_ecl;
double xs, ys;
double xe, ye, ze;
/* Compute Sun's ecliptical coordinates */
sunpos (d, &lon, r);
/* Compute ecliptic rectangular coordinates */
xs = *r * cosd(lon);
ys = *r * sind(lon);
/* Compute obliquity of ecliptic (inclination of Earth's axis) */
obl_ecl = 23.4393 - 3.563E-7 * d;
/* Convert to equatorial rectangular coordinates - x is unchanged */
xe = xs;
ye = ys * cosd(obl_ecl);
ze = ys * sind(obl_ecl);
/* Convert to spherical coordinates */
*RA = atan2d(ye, xe);
*dec = atan2d(ze, sqrt(xe*xe + ye*ye));
}
//
// Utility functions
//
// Reduce angle to within 0..359.999 degrees
double revolution (const double x)
{ double remainder = fmod (x, (double) 360.0);
return remainder < (double) 0.0 ? remainder + (double) 360.0 : remainder;
}
// Reduce angle to -179.999 to +180 degrees
double rev180 (const double x)
{ double y = revolution (x);
return y <= (double) 180.0 ? y : y - (double) 360.0;
}
// Fix angle to 0-359.999
double fixLongitude (const double x)
{ return revolution (x);
}
// Fix angle to 0-89.999 and -0.001 to -89.999
double fixLatitude (const double x)
{
// Make angle 0 to 359.9999
double y = revolution (x);
if (y <= (double) 90.0) ;
else if (y <= (double) 180.0) y = (double) 180.0 - y;
else if (y <= (double) 270.0) y = (double) 180.0 - y;
else if (y <= (double) 360.0) y = y - (double) 360.0;
// Linux compile of sunwait doesn't like 90, Windows is OK.
// Let's just wiggle things a little bit to make things OK.
if (y == (double) 90.0) y = (double) 89.9999999;
else if (y == (double) -90.0) y = (double) -89.9999999;
return y;
}
// Time must be between 0:00 amd 23:59
double fix24 (const double x)
{
double remainder = fmod (x, (double) 24.0);
return remainder < (double) 0.0 ? remainder + (double) 24.0 : remainder;
}
/*******************************************************************/
/* This function computes GMST0, the Greenwhich Mean Sidereal Time */
/* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at */
/* 0h UT). GMST is then the sidereal time at Greenwich at any */
/* time of the day. I've generalized GMST0 as well, and define it */
/* as: GMST0 = GMST - UT -- this allows GMST0 to be computed at */
/* other times than 0h UT as well. While this sounds somewhat */
/* contradictory, it is very practical: instead of computing */
/* GMST like: */
/* */
/* GMST = (GMST0) + UT * (366.2422/365.2422) */
/* */
/* where (GMST0) is the GMST last time UT was 0 hours, one simply */
/* computes: */
/* */
/* GMST = GMST0 + UT */
/* */
/* where GMST0 is the GMST "at 0h UT" but at the current moment! */
/* Defined in this way, GMST0 will increase with about 4 min a */
/* day. It also happens that GMST0 (in degrees, 1 hr = 15 degr) */
/* is equal to the Sun's mean longitude plus/minus 180 degrees! */
/* (if we neglect aberration, which amounts to 20 seconds of arc */
/* or 1.33 seconds of time) */
/* */
/*******************************************************************/
inline double GMST0 (const double d)
{
/* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr */
/* L = M + w, as defined in sunpos(). Since I'm too lazy to */
/* add these numbers, I'll let the C compiler do it for me. */
/* Any decent C compiler will add the constants at compile */
/* time, imposing no runtime or code overhead. */
return revolution ((180.0 + 356.0470 + 282.9404) + (0.9856002585 + 4.70935E-5) * d);
}
unsigned long daysSince2000 (const time_t *pTimet)
{
struct tm tmpTm;
myUtcTime (pTimet, &tmpTm);
unsigned int yearsSince2000 = tmpTm.tm_year - 100; // Get year, but tm_year starts from 1900
// Calucate number of leap days, but -
// yearsSince2000 - 1
// Don't include this year as tm_yday includes this year's leap day in the next bit
unsigned int leapDaysSince2000
= (unsigned int) floor ((yearsSince2000-1)/4) // Every evenly divisible 4 years is a leap-year
- (unsigned int) floor ((yearsSince2000-1)/100) // Except centuries
+ (unsigned int) floor ((yearsSince2000-1)/400) // Unless evenlt divisible by 400
+ 1; // 2000 itself was a leap year with the 400 rule (a fix for 0/400 == 0)
return (yearsSince2000 * 365) + leapDaysSince2000 + tmpTm.tm_yday;
}
/*
** Utility functions
*/
long myRound (const double d) { return d > 0.0 ? (int) (d + 0.5) : (int) (d - 0.5) ; }
long myTrunc (const double d) { return (d>0) ? (int) floor(d) : (int) ceil(d) ; }
double myAbs (const double d) { return (d>0) ? d : -d ; }
int hours (const double d) { return myTrunc (d); }
int minutes (const double d) { return myTrunc (fmod (myAbs (d) * 60, 60)); }