/[zanavi_public1]/navit/navit/sunriset.c
ZANavi

Contents of /navit/navit/sunriset.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2 - (show annotations) (download)
Fri Oct 28 21:19:04 2011 UTC (12 years, 5 months ago) by zoff99
File MIME type: text/plain
File size: 12799 byte(s)
import files
1 /*
2
3 SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
4 the length of the day at any date and latitude
5
6 Written as DAYLEN.C, 1989-08-16
7
8 Modified to SUNRISET.C, 1992-12-01
9
10 (c) Paul Schlyter, 1989, 1992
11
12 Released to the public domain by Paul Schlyter, December 1992
13
14 */
15
16
17 #include <stdio.h>
18 #include <math.h>
19
20 #include "sunriset.h"
21
22 /* The "workhorse" function for sun rise/set times */
23
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 */
38 /* twilight. */
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. */
53 /* */
54 /**********************************************************************/
55 {
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 */
61 t, /* Diurnal arc */
62 tsouth, /* Time when Sun is at south */
63 sidtime; /* Local sidereal time */
64
65 int rc = 0; /* Return cde from function - usually 0 */
66
67 /* Compute d of 12h local mean solar time */
68 d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
69
70 /* Compute local sideral time of this moment */
71 sidtime = revolution( GMST0(d) + 180.0 + lon );
72
73 /* Compute Sun's RA + Decl at this moment */
74 sun_RA_dec( d, &sRA, &sdec, &sr );
75
76 /* Compute time when Sun is at south - in hours UT */
77 tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
78
79 /* Compute the Sun's apparent radius, degrees */
80 sradius = 0.2666 / sr;
81
82 /* Do correction to upper limb, if necessary */
83 if ( upper_limb )
84 altit -= sradius;
85
86 /* Compute the diurnal arc that the Sun traverses to reach */
87 /* the specified altitide altit: */
88 {
89 double cost;
90 cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
91 ( cosd(lat) * cosd(sdec) );
92 if ( cost >= 1.0 )
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 */
96 else
97 t = acosd(cost)/15.0; /* The diurnal arc, hours */
98 }
99
100 /* Store rise and set times - in hours UT */
101 *trise = tsouth - t;
102 *tset = tsouth + t;
103
104 return rc;
105 } /* __sunriset__ */
106
107
108
109 /* The "workhorse" function */
110
111
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 /**********************************************************************/
129 {
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 */
137 t; /* Diurnal arc */
138
139 /* Compute d of 12h local mean solar time */
140 d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
141
142 /* Compute obliquity of ecliptic (inclination of Earth's axis) */
143 obl_ecl = 23.4393 - 3.563E-7 * d;
144
145 /* Compute Sun's position */
146 sunpos( d, &slon, &sr );
147
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 );
151
152 /* Compute the Sun's apparent radius, degrees */
153 sradius = 0.2666 / sr;
154
155 /* Do correction to upper limb, if necessary */
156 if ( upper_limb )
157 altit -= sradius;
158
159 /* Compute the diurnal arc that the Sun traverses to reach */
160 /* the specified altitide altit: */
161 {
162 double cost;
163 cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
164 ( cosd(lat) * cos_sdecl );
165 if ( cost >= 1.0 )
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 */
170 }
171 return t;
172 } /* __daylen__ */
173
174
175 /* This function computes the Sun's position at any instant */
176
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 /******************************************************/
184 {
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 */
192
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;
197
198 /* Compute true longitude and radius vector */
199 E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
200 x = cosd(E) - e;
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 */
205 if ( *lon >= 360.0 )
206 *lon -= 360.0; /* Make it 0..360 degrees */
207 }
208
209 void sun_RA_dec( double d, double *RA, double *dec, double *r )
210 {
211 double lon, obl_ecl;
212 double xs, ys, zs;
213 double xe, ye, ze;
214
215 /* Compute Sun's ecliptical coordinates */
216 sunpos( d, &lon, r );
217
218 /* Compute ecliptic rectangular coordinates */
219 xs = *r * cosd(lon);
220 ys = *r * sind(lon);
221 zs = 0; /* because the Sun is always in the ecliptic plane! */
222
223 /* Compute obliquity of ecliptic (inclination of Earth's axis) */
224 obl_ecl = 23.4393 - 3.563E-7 * d;
225
226 /* Convert to equatorial rectangular coordinates - x is unchanged */
227 xe = xs;
228 ye = ys * cosd(obl_ecl);
229 ze = ys * sind(obl_ecl);
230
231 /* Convert to spherical coordinates */
232 *RA = atan2d( ye, xe );
233 *dec = atan2d( ze, sqrt(xe*xe + ye*ye) );
234
235 } /* sun_RA_dec */
236
237
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 /******************************************************************/
243
244 #define INV360 ( 1.0 / 360.0 )
245
246 double revolution( double x )
247 /*****************************************/
248 /* Reduce angle to within 0..360 degrees */
249 /*****************************************/
250 {
251 return( x - 360.0 * floor( x * INV360 ) );
252 } /* revolution */
253
254 double rev180( double x )
255 /*********************************************/
256 /* Reduce angle to within -180..+180 degrees */
257 /*********************************************/
258 {
259 return( x - 360.0 * floor( x * INV360 + 0.5 ) );
260 } /* revolution */
261
262
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 */
271 /* GMST like: */
272 /* */
273 /* GMST = (GMST0) + UT * (366.2422/365.2422) */
274 /* */
275 /* where (GMST0) is the GMST last time UT was 0 hours, one simply */
276 /* computes: */
277 /* */
278 /* GMST = GMST0 + UT */
279 /* */
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) */
286 /* */
287 /*******************************************************************/
288
289 double GMST0( double d )
290 {
291 double sidtim0;
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 );
299 return sidtim0;
300 } /* GMST0 */

   
Visit the ZANavi Wiki