3 SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
4 the length of the day at any date and latitude
6 Written as DAYLEN.C, 1989-08-16
8 Modified to SUNRISET.C, 1992-12-01
10 (c) Paul Schlyter, 1989, 1992
12 Released to the public domain by Paul Schlyter, December 1992
22 /* The "workhorse" function for sun rise/set times */
24 int __sunriset__( int year, int month, int day, double lon, double lat,
25 double altit, int upper_limb, double *trise, double *tset )
26 /***************************************************************************/
27 /* Note: year,month,date = calendar date, 1801-2099 only. */
28 /* Eastern longitude positive, Western longitude negative */
29 /* Northern latitude positive, Southern latitude negative */
30 /* The longitude value IS critical in this function! */
31 /* altit = the altitude which the Sun should cross */
32 /* Set to -35/60 degrees for rise/set, -6 degrees */
33 /* for civil, -12 degrees for nautical and -18 */
34 /* degrees for astronomical twilight. */
35 /* upper_limb: non-zero -> upper limb, zero -> center */
36 /* Set to non-zero (e.g. 1) when computing rise/set */
37 /* times, and to zero when computing start/end of */
39 /* *rise = where to store the rise time */
40 /* *set = where to store the set time */
41 /* Both times are relative to the specified altitude, */
42 /* and thus this function can be used to comupte */
43 /* various twilight times, as well as rise/set times */
44 /* Return value: 0 = sun rises/sets this day, times stored at */
45 /* *trise and *tset. */
46 /* +1 = sun above the specified "horizon" 24 hours. */
47 /* *trise set to time when the sun is at south, */
48 /* minus 12 hours while *tset is set to the south */
49 /* time plus 12 hours. "Day" length = 24 hours */
50 /* -1 = sun is below the specified "horizon" 24 hours */
51 /* "Day" length = 0 hours, *trise and *tset are */
52 /* both set to the time when the sun is at south. */
54 /**********************************************************************/
56 double d, /* Days since 2000 Jan 0.0 (negative before) */
57 sr, /* Solar distance, astronomical units */
58 sRA, /* Sun's Right Ascension */
59 sdec, /* Sun's declination */
60 sradius, /* Sun's apparent radius */
62 tsouth, /* Time when Sun is at south */
63 sidtime; /* Local sidereal time */
65 int rc = 0; /* Return cde from function - usually 0 */
67 /* Compute d of 12h local mean solar time */
68 d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
70 /* Compute local sideral time of this moment */
71 sidtime = revolution( GMST0(d) + 180.0 + lon );
73 /* Compute Sun's RA + Decl at this moment */
74 sun_RA_dec( d, &sRA, &sdec, &sr );
76 /* Compute time when Sun is at south - in hours UT */
77 tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
79 /* Compute the Sun's apparent radius, degrees */
80 sradius = 0.2666 / sr;
82 /* Do correction to upper limb, if necessary */
86 /* Compute the diurnal arc that the Sun traverses to reach */
87 /* the specified altitide altit: */
90 cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
91 ( cosd(lat) * cosd(sdec) );
93 rc = -1, t = 0.0; /* Sun always below altit */
94 else if ( cost <= -1.0 )
95 rc = +1, t = 12.0; /* Sun always above altit */
97 t = acosd(cost)/15.0; /* The diurnal arc, hours */
100 /* Store rise and set times - in hours UT */
109 /* The "workhorse" function */
112 double __daylen__( int year, int month, int day, double lon, double lat,
113 double altit, int upper_limb )
114 /**********************************************************************/
115 /* Note: year,month,date = calendar date, 1801-2099 only. */
116 /* Eastern longitude positive, Western longitude negative */
117 /* Northern latitude positive, Southern latitude negative */
118 /* The longitude value is not critical. Set it to the correct */
119 /* longitude if you're picky, otherwise set to to, say, 0.0 */
120 /* The latitude however IS critical - be sure to get it correct */
121 /* altit = the altitude which the Sun should cross */
122 /* Set to -35/60 degrees for rise/set, -6 degrees */
123 /* for civil, -12 degrees for nautical and -18 */
124 /* degrees for astronomical twilight. */
125 /* upper_limb: non-zero -> upper limb, zero -> center */
126 /* Set to non-zero (e.g. 1) when computing day length */
127 /* and to zero when computing day+twilight length. */
128 /**********************************************************************/
130 double d, /* Days since 2000 Jan 0.0 (negative before) */
131 obl_ecl, /* Obliquity (inclination) of Earth's axis */
132 sr, /* Solar distance, astronomical units */
133 slon, /* True solar longitude */
134 sin_sdecl, /* Sine of Sun's declination */
135 cos_sdecl, /* Cosine of Sun's declination */
136 sradius, /* Sun's apparent radius */
139 /* Compute d of 12h local mean solar time */
140 d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
142 /* Compute obliquity of ecliptic (inclination of Earth's axis) */
143 obl_ecl = 23.4393 - 3.563E-7 * d;
145 /* Compute Sun's position */
146 sunpos( d, &slon, &sr );
148 /* Compute sine and cosine of Sun's declination */
149 sin_sdecl = sind(obl_ecl) * sind(slon);
150 cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
152 /* Compute the Sun's apparent radius, degrees */
153 sradius = 0.2666 / sr;
155 /* Do correction to upper limb, if necessary */
159 /* Compute the diurnal arc that the Sun traverses to reach */
160 /* the specified altitide altit: */
163 cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
164 ( cosd(lat) * cos_sdecl );
166 t = 0.0; /* Sun always below altit */
167 else if ( cost <= -1.0 )
168 t = 24.0; /* Sun always above altit */
169 else t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
175 /* This function computes the Sun's position at any instant */
177 void sunpos( double d, double *lon, double *r )
178 /******************************************************/
179 /* Computes the Sun's ecliptic longitude and distance */
180 /* at an instant given in d, number of days since */
181 /* 2000 Jan 0.0. The Sun's ecliptic latitude is not */
182 /* computed, since it's always very near 0. */
183 /******************************************************/
185 double M, /* Mean anomaly of the Sun */
186 w, /* Mean longitude of perihelion */
187 /* Note: Sun's mean longitude = M + w */
188 e, /* Eccentricity of Earth's orbit */
189 E, /* Eccentric anomaly */
190 x, y, /* x, y coordinates in orbit */
191 v; /* True anomaly */
193 /* Compute mean elements */
194 M = revolution( 356.0470 + 0.9856002585 * d );
195 w = 282.9404 + 4.70935E-5 * d;
196 e = 0.016709 - 1.151E-9 * d;
198 /* Compute true longitude and radius vector */
199 E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
201 y = sqrt( 1.0 - e*e ) * sind(E);
202 *r = sqrt( x*x + y*y ); /* Solar distance */
203 v = atan2d( y, x ); /* True anomaly */
204 *lon = v + w; /* True solar longitude */
206 *lon -= 360.0; /* Make it 0..360 degrees */
209 void sun_RA_dec( double d, double *RA, double *dec, double *r )
215 /* Compute Sun's ecliptical coordinates */
216 sunpos( d, &lon, r );
218 /* Compute ecliptic rectangular coordinates */
221 zs = 0; /* because the Sun is always in the ecliptic plane! */
223 /* Compute obliquity of ecliptic (inclination of Earth's axis) */
224 obl_ecl = 23.4393 - 3.563E-7 * d;
226 /* Convert to equatorial rectangular coordinates - x is unchanged */
228 ye = ys * cosd(obl_ecl);
229 ze = ys * sind(obl_ecl);
231 /* Convert to spherical coordinates */
232 *RA = atan2d( ye, xe );
233 *dec = atan2d( ze, sqrt(xe*xe + ye*ye) );
238 /******************************************************************/
239 /* This function reduces any angle to within the first revolution */
240 /* by subtracting or adding even multiples of 360.0 until the */
241 /* result is >= 0.0 and < 360.0 */
242 /******************************************************************/
244 #define INV360 ( 1.0 / 360.0 )
246 double revolution( double x )
247 /*****************************************/
248 /* Reduce angle to within 0..360 degrees */
249 /*****************************************/
251 return( x - 360.0 * floor( x * INV360 ) );
254 double rev180( double x )
255 /*********************************************/
256 /* Reduce angle to within -180..+180 degrees */
257 /*********************************************/
259 return( x - 360.0 * floor( x * INV360 + 0.5 ) );
263 /*******************************************************************/
264 /* This function computes GMST0, the Greenwhich Mean Sidereal Time */
265 /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at */
266 /* 0h UT). GMST is then the sidereal time at Greenwich at any */
267 /* time of the day. I've generelized GMST0 as well, and define it */
268 /* as: GMST0 = GMST - UT -- this allows GMST0 to be computed at */
269 /* other times than 0h UT as well. While this sounds somewhat */
270 /* contradictory, it is very practical: instead of computing */
273 /* GMST = (GMST0) + UT * (366.2422/365.2422) */
275 /* where (GMST0) is the GMST last time UT was 0 hours, one simply */
278 /* GMST = GMST0 + UT */
280 /* where GMST0 is the GMST "at 0h UT" but at the current moment! */
281 /* Defined in this way, GMST0 will increase with about 4 min a */
282 /* day. It also happens that GMST0 (in degrees, 1 hr = 15 degr) */
283 /* is equal to the Sun's mean longitude plus/minus 180 degrees! */
284 /* (if we neglect aberration, which amounts to 20 seconds of arc */
285 /* or 1.33 seconds of time) */
287 /*******************************************************************/
289 double GMST0( double d )
292 /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr */
293 /* L = M + w, as defined in sunpos(). Since I'm too lazy to */
294 /* add these numbers, I'll let the C compiler do it for me. */
295 /* Any decent C compiler will add the constants at compile */
296 /* time, imposing no runtime or code overhead. */
297 sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
298 ( 0.9856002585 + 4.70935E-5 ) * d );