/[ascend]/trunk/models/johnpye/datareader/sun.c
ViewVC logotype

Contents of /trunk/models/johnpye/datareader/sun.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 816 - (show annotations) (download) (as text)
Tue Aug 15 14:40:33 2006 UTC (13 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 32042 byte(s)
Adding the 'sun.c' and 'sun.h' files back into ASCEND so that we don't need dependency on libradtran.
Added detection of Scrollkeeper, the freedesktop.org documentation manager.
Added placeholder SConscript file for building the user's manual.
1 /*--------------------------------------------------------------------
2 * $Id: sun.c,v 1.6 2004/07/05 06:57:56 pa2h Exp $
3 *
4 * This file is part of libRadtran.
5 * Copyright (c) 1997-2001 by Arve Kylling and Bernhard Mayer.
6 *
7 * ######### Contact info: http://www.libradtran.org #########
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place - Suite 330,
22 * Boston, MA 02111-1307, USA.
23 *--------------------------------------------------------------------*/
24
25
26 /* sun.a - Solar zenith and azimuth calculations @ti@ */
27 /* */
28 /* Most of the formulas are adopted from: */
29 /* Iqbal, Muhammad: "An Introduction to Solar Radiation", */
30 /* Academic Press, Inc., 1983 */
31 /* (page numbers in the comments refer to this book). */
32 /* */
33 /* Time is specified in seconds from midnight, */
34 /* angles are specified in degrees. */
35
36 #include <math.h>
37 #include <stdio.h>
38 #include <string.h>
39 #include "sun.h"
40
41
42
43 /* @42c@ */
44 /* @c42@ */
45 /****************************************************************************/
46 /* In order to use the functions provided by the sun library, @42_10c@ */
47 /* #include <sun.h> in your source code and link with libRadtran_c.a. */
48 /* */
49 /* @strong{Example:} */
50 /* Example for a source file: */
51 /* @example */
52 /* */
53 /* ... */
54 /* #include "../src_c/sun.h" */
55 /* ... */
56 /* */
57 /* @end example */
58 /* */
59 /* Linking of the executable, using the GNU compiler gcc: */
60 /* @example */
61 /* */
62 /* gcc -o test test.c -lRadtran_c -L../lib */
63 /* */
64 /* @end example @c42_10@ */
65 /****************************************************************************/
66
67
68 /****************************************************************************/
69 /* The sun library provides functions for solar zenith and azimuth @42_20c@ */
70 /* angle and sun-earth-distance calculations. All formulas have been taken */
71 /* from Iqbal, "An introduction to solar radiation". */
72 /* @c42_20@ */
73 /****************************************************************************/
74
75
76 /* define _PI_ */
77 #define _PI_ 3.1415926
78
79 /* internal functions */
80 static double dayangle (int day);
81
82
83 /***********************************************************************************/
84 /* Function: dayangle */
85 /* Description: */
86 /* Internal function to calculate the dayangle according to Iqbal, pg. 3 */
87 /* */
88 /* Parameters: */
89 /* int day: day of year. */
90 /* */
91 /* Return value: */
92 /* double dayangle, 0..2pi */
93 /* */
94 /* Example: */
95 /* Files: */
96 /* Known bugs: */
97 /* Author: */
98 /* */
99 /***********************************************************************************/
100
101 static double dayangle (int day)
102 {
103 return 2.0 * _PI_ * (double) (day-1) / 365.0;
104 }
105
106
107
108 /***********************************************************************************/
109 /* Function: eccentricity @42_30i@ */
110 /* Description: */
111 /* Calculate the eccentricity correction factor E0 = (r0/r)**2 according to */
112 /* Iqbal, page 3. This factor, when multiplied with the irradiance, accounts */
113 /* for the annual variation of the sun-earth-distance. */
114 /* */
115 /* Parameters: */
116 /* int day: day of year (leap day is usually @strong{not} counted. */
117 /* */
118 /* Return value: */
119 /* The eccentricity (double) for the specified day. */
120 /* */
121 /* Example: */
122 /* Files: */
123 /* Known bugs: */
124 /* Author: */
125 /* @i42_30@ */
126 /***********************************************************************************/
127
128 double eccentricity (int day)
129 {
130 double E0=0, angle=0;
131
132 angle = dayangle (day);
133
134 E0 = 1.000110 + 0.034221 * cos(angle) + 0.001280 * sin(angle)
135 + 0.000719 * cos(2*angle) + 0.000077 * sin(2*angle);
136
137 return E0;
138 }
139
140
141
142 /***********************************************************************************/
143 /* Function: declination @42_30i@ */
144 /* Description: */
145 /* Calculate the declination for a specified day (Iqbal, page 7). */
146 /* */
147 /* Parameters: */
148 /* int day: day of year (leap day is usually @strong{not} counted. */
149 /* */
150 /* Return value: */
151 /* The declination in degrees (double) for the specified day. */
152 /* */
153 /* Example: */
154 /* Files: */
155 /* Known bugs: */
156 /* Author: */
157 /* @i42_30@ */
158 /***********************************************************************************/
159
160 double declination (int day)
161 {
162 double delta=0, angle=0;
163
164 angle = dayangle (day);
165
166 delta = 0.006918 - 0.399912 * cos (angle) + 0.070257 * sin (angle)
167 - 0.006758 * cos (2*angle) + 0.000907 * sin (2*angle)
168 - 0.002697 * cos (3*angle) + 0.00148 * sin (3*angle);
169
170 delta *= 180.0/_PI_;
171
172 return delta;
173 }
174
175
176
177 /***********************************************************************************/
178 /* Function: equation_of_time @42_30i@ */
179 /* Description: */
180 /* Calculate the equation of time for a specified day (Iqbal, page 11). */
181 /* */
182 /* Parameters: */
183 /* int day: day of year (leap day is usually @strong{not} counted. */
184 /* */
185 /* Return value: */
186 /* The equation of time in seconds (double) for the specified day. */
187 /* */
188 /* Example: */
189 /* Files: */
190 /* Known bugs: */
191 /* Author: */
192 /* @i42_30@ */
193 /***********************************************************************************/
194
195 int equation_of_time (int day)
196 {
197 double angle=0, et=0;
198
199 angle = dayangle (day);
200
201 et = (0.000075 + 0.001868 * cos(angle) - 0.032077 * sin(angle)
202 -0.014615 * cos(2*angle) - 0.04089 * sin(2*angle)) * 13750.8;
203
204 return (int) (et+0.5);
205 }
206
207
208
209 /***********************************************************************************/
210 /* Function: LAT @42_30i@ */
211 /* Description: */
212 /* Calculate the local apparent time for a given standard time and location. */
213 /* */
214 /* Parameters: */
215 /* int time_std: Standard time [seconds since midnight]. */
216 /* int day: day of year (leap day is usually @strong{not} counted. */
217 /* double longitude: Longitude [degrees] (West positive). */
218 /* double long_std: Standard longitude [degrees]. */
219 /* */
220 /* Return value: */
221 /* The local apparent time in seconds since midnight (double). */
222 /* */
223 /* Example: */
224 /* Files: */
225 /* Known bugs: */
226 /* Author: */
227 /* @i42_30@ */
228 /***********************************************************************************/
229
230 int LAT (int time_std, int day, double longitude, double long_std)
231 {
232 int lat=0;
233
234 lat = time_std + (int) (240.0 * (long_std-longitude))
235 + equation_of_time (day);
236
237 return lat;
238 }
239
240
241 /************************************************/
242 /* convert local apparent time to standard time */
243 /************************************************/
244
245 int standard_time (int lat, int day, double longitude, double long_std)
246 {
247 return lat - (int) (240.0 * (long_std-longitude)) - equation_of_time (day);
248 }
249
250
251 /****************************/
252 /* hour angle omega; pg. 15 */
253 /****************************/
254
255 double hour_angle (int time)
256 {
257 double omega=0;
258
259 omega = _PI_ * (1.0 - ((double) time) / 43200.0);
260
261 omega *= 180.0/_PI_;
262 return omega;
263 }
264
265
266
267 /***********************************************************************************/
268 /* Function: solar_zenith @42_30i@ */
269 /* Description: */
270 /* Calculate the solar zenith angle for a given time and location. */
271 /* */
272 /* Parameters: */
273 /* int time: Standard time [seconds since midnight]. */
274 /* int day: day of year (leap day is usually @strong{not} counted. */
275 /* double latitude: Latitude [degrees] (North positive). */
276 /* double longitude: Longitude [degrees] (West positive). */
277 /* double long_std: Standard longitude [degrees]. */
278 /* */
279 /* Return value: */
280 /* The solar zenith angle [degrees]. */
281 /* */
282 /* Example: */
283 /* Files: */
284 /* Known bugs: */
285 /* Author: */
286 /* @i42_30@ */
287 /***********************************************************************************/
288
289 double solar_zenith (int time, int day,
290 double latitude, double longitude, double long_std)
291 {
292 double theta=0, omega=0, delta=0, phi = latitude*_PI_/180.0;
293 int lat = LAT (time, day, longitude, long_std);
294
295 delta = _PI_/180.0*declination (day);
296 omega = _PI_/180.0*hour_angle (lat);
297
298 theta = acos(sin(delta) * sin(phi) + cos(delta) * cos(phi) * cos(omega));
299
300 theta *= (180.0/_PI_); /* convert to degrees */
301 return theta;
302 }
303
304
305
306
307 /***********************************************************************************/
308 /* Function: solar_azimuth @42_30i@ */
309 /* Description: */
310 /* Calculate the solar azimuth angle for a given time and location. */
311 /* */
312 /* Parameters: */
313 /* int time: Standard time [seconds since midnight]. */
314 /* int day: day of year (leap day is usually @strong{not} counted. */
315 /* double latitude: Latitude [degrees] (North positive). */
316 /* double longitude: Longitude [degrees] (West positive). */
317 /* double long_std: Standard longitude [degrees]. */
318 /* */
319 /* Return value: */
320 /* The solar azimuth angle [degrees]. */
321 /* */
322 /* Example: */
323 /* Files: */
324 /* Known bugs: */
325 /* Author: */
326 /* @i42_30@ */
327 /***********************************************************************************/
328
329 double solar_azimuth (int time, int day,
330 double latitude, double longitude, double long_std)
331 {
332 double theta=0, omega=0, delta=0, phi = latitude*_PI_/180.0, psi=0, cospsi=0;
333 int lat = LAT (time, day, longitude, long_std);
334
335 delta = _PI_/180.0*declination (day);
336 omega = _PI_/180.0*hour_angle (lat);
337 theta = _PI_/180.0*solar_zenith (time, day, latitude, longitude, long_std);
338
339
340 cospsi = (cos(theta)*sin(phi) - sin(delta)) / sin(theta) / cos (phi);
341 /* allow tiny roundoff errors */
342 if (cospsi>1.0 && cospsi<1.0+1e-6)
343 cospsi=1.0;
344
345 if (cospsi<-1.0 && cospsi>-1.0-1e-6)
346 cospsi=-1.0;
347
348 psi = -acos(cospsi);
349
350 /* adjust sign */
351 if (lat>43200 || lat<0)
352 psi=-psi;
353
354 psi *= (180.0/_PI_); /* convert to degrees */
355
356 return psi;
357 }
358
359
360
361 /***********************************************************************************/
362 /* Function: day_of_year @42_30i@ */
363 /* Description: */
364 /* Calculate the day of year for given date (leap days are not considered. */
365 /* */
366 /* Parameters: */
367 /* int day: Day of month (1..31). */
368 /* int month: Month (1..12). */
369 /* */
370 /* Return value: */
371 /* The day of year (int); -1 if error. */
372 /* */
373 /* Example: */
374 /* Files: */
375 /* Known bugs: */
376 /* Author: */
377 /* @i42_30@ */
378 /***********************************************************************************/
379
380 int day_of_year (int day, int month)
381 {
382 char mon = 0;
383 int doy = -1;
384
385 if (month<1 || month > 12)
386 return (-1);
387
388 if (day<1 || day > 31)
389 return (-1);
390
391 mon = (char) month;
392
393 switch (mon) {
394 case 1:
395 doy = day + 0;
396 break;
397 case 2:
398 doy = day + 31;
399 break;
400 case 3:
401 doy = day + 59;
402 break;
403 case 4:
404 doy = day + 90;
405 break;
406 case 5:
407 doy = day + 120;
408 break;
409 case 6:
410 doy = day + 151;
411 break;
412 case 7:
413 doy = day + 181;
414 break;
415 case 8:
416 doy = day + 212;
417 break;
418 case 9:
419 doy = day + 243;
420 break;
421 case 10:
422 doy = day + 273;
423 break;
424 case 11:
425 doy = day + 304;
426 break;
427 case 12:
428 doy = day + 334;
429 break;
430 default:
431 doy = -1;
432 }
433
434 return (doy);
435 }
436
437
438
439
440 /**************************************************************************/
441 /* convert time to string */
442 /* memory for timestr must be allocated by programmer (at least 10 bytes) */
443 /**************************************************************************/
444
445 char *time2str (char *timestr, int hour, int min, int sec)
446 {
447 char hourstr[3] = "";
448 char minstr[3] = "";
449 char secstr[3] = "";
450
451 strcpy (timestr, "");
452
453 if (hour<0 || hour>24)
454 return NULL;
455
456 if (min<0 || min>60)
457 return NULL;
458
459 if (sec<0 || sec>60)
460 return NULL;
461
462
463
464 if (hour<10)
465 sprintf (hourstr, "0%d", hour);
466 else
467 sprintf (hourstr, "%d", hour);
468
469 if (min<10)
470 sprintf (minstr, "0%d", min);
471 else
472 sprintf (minstr, "%d", min);
473
474 if (sec<10)
475 sprintf (secstr, "0%d", sec);
476 else
477 sprintf (secstr, "%d", sec);
478
479 sprintf (timestr, "%s:%s:%s", hourstr, minstr, secstr);
480
481 return timestr;
482 }
483
484
485
486
487 /***********************************************************************************/
488 /* Function: zenith2time @42_30i@ */
489 /* Description: */
490 /* Calculate the times for a given solar zenith angle, day of year and location. */
491 /* */
492 /* Parameters: */
493 /* int day: day of year */
494 /* double zenith_angle: Solar zenith angle [degrees]. */
495 /* double latitude: Latitude [degrees] (North positive). */
496 /* double longitude: Longitude [degrees] (West positive). */
497 /* double long_std: Standard longitude [degrees]. */
498 /* int *time1: 1st time of occurence. */
499 /* int *time2: 2nd time of occurence. */
500 /* */
501 /* Return value: */
502 /* 0 if o.k., <0 if error. */
503 /* */
504 /* Example: */
505 /* Files: */
506 /* Known bugs: */
507 /* Author: */
508 /* @i42_30@ */
509 /***********************************************************************************/
510
511 int zenith2time (int day,
512 double zenith_angle,
513 double latitude,
514 double longitude,
515 double long_std,
516 int *time1,
517 int *time2)
518 {
519 double delta = _PI_/180.0*declination (day);
520 double phi = _PI_/180.0*latitude;
521 double theta = _PI_/180.0*zenith_angle;
522 double cos_omega = (cos(theta) - sin(delta)*sin(phi)) /
523 cos(delta) / cos(phi);
524 double omega1=0, omega2=0;
525 int lat1=0, lat2=0;
526
527 if (fabs(cos_omega) > 1.0)
528 return ERROR_NO_ZENITH;
529
530 omega1 = acos (cos_omega);
531 omega2 = 0.0-omega1;
532
533 lat1 = 43200*(1-omega1/_PI_);
534 lat2 = 43200*(1-omega2/_PI_);
535
536 *time1 = standard_time (lat1, day, longitude, long_std);
537 *time2 = standard_time (lat2, day, longitude, long_std);
538
539 return 0;
540 }
541
542
543
544
545
546 /***********************************************************************************/
547 /* Function: Gregorian2Julian @42_30i@ */
548 /* Description: */
549 /* Convert from Gregorian day (day, month, year) to Julian day (by the */
550 /* astronomical definition). This function, in combination with */
551 /* Julian2Gregorian() is very useful to answer questions like "which date is */
552 /* 666 days after December 31, 1999?" Algorithm from */
553 /* H.F. Fliegel and T.C. Van Flandern, "A Machine Algorithm for Processing */
554 /* Calendar Dates", Communications of the Association for Computing Machinery */
555 /* (CACM), Vol. 11, No. 10, 657, 1968. */
556 /* */
557 /* Parameters: */
558 /* int d: Day of month (1..31). */
559 /* int m: Month (1..12). */
560 /* int y: Year (attention: full year required, 1999 instead of 99) */
561 /* int *jd: The Julian day, to be calculated. */
562 /* */
563 /* Return value: */
564 /* 0 if o.k., <0 if error. */
565 /* */
566 /* Example: */
567 /* Files: */
568 /* Known bugs: */
569 /* Author: */
570 /* @i42_30@ */
571 /***********************************************************************************/
572
573 int Gregorian2Julian (int d, int m, int y, int *jd)
574 {
575 char ystr[3]="";
576
577 /* reset output */
578 *jd =0;
579
580 if (d<1||d>31)
581 return -1;
582
583 if (m<1||m>12)
584 return -1;
585
586 if (y>0&&y<100) {
587 if (y<10)
588 sprintf (ystr, "0%d", y);
589 else
590 sprintf (ystr, "%d", y);
591
592 fprintf (stderr, "Warning: Gregorian2Julian() has been called for year %s\n", ystr);
593 fprintf (stderr, "Maybe you meant 19%s or 20%s?\n", ystr, ystr);
594 }
595
596
597 *jd = (1461 * (y + 4800 + (m - 14) / 12)) / 4 +
598 (367 * (m - 2 - 12 * ((m - 14) / 12))) / 12 -
599 (3 * ((y + 4900 + (m - 14) / 12) / 100)) / 4 + d - 32075;
600
601 return 0; /* if o.k. */
602 }
603
604
605
606 /***********************************************************************************/
607 /* Function: Julian2Gregorian @42_30i@ */
608 /* Description: */
609 /* Convert from Julian day (by the astronomical definition) to Gregorian day */
610 /* (day, month, year) to . This function, in combination with */
611 /* Gregorian2Julian() is very useful to answer questions like "which date is */
612 /* 666 days after December 31, 1999?" Algorithm from */
613 /* H.F. Fliegel and T.C. Van Flandern, "A Machine Algorithm for Processing */
614 /* Calendar Dates", Communications of the Association for Computing Machinery */
615 /* (CACM), Vol. 11, No. 10, 657, 1968. */
616 /* */
617 /* Parameters: */
618 /* int *d: Day of month (1..31), to be calculated. */
619 /* int *m: Month (1..12), to be calculated. */
620 /* int *y: Year, to be calculated. */
621 /* int jd: The Julian day. */
622 /* */
623 /* Return value: */
624 /* 0 if o.k., <0 if error. */
625 /* */
626 /* Example: */
627 /* Files: */
628 /* Known bugs: */
629 /* Author: */
630 /* @i42_30@ */
631 /***********************************************************************************/
632
633 int Julian2Gregorian(int *d, int *m, int *y, int jd)
634 {
635 int l=0, n=0, i=0, j=0;
636
637 l = jd + 68569;
638 n = ( 4 * l ) / 146097;
639 l = l - ( 146097 * n + 3 ) / 4;
640 i = ( 4000 * ( l + 1 ) ) / 1461001;
641 l = l - ( 1461 * i ) / 4 + 31;
642 j = ( 80 * l ) / 2447;
643 *d = l - ( 2447 * j ) / 80;
644 l = j / 11;
645 *m = j + 2 - ( 12 * l );
646 *y = 100 * ( n - 49 ) + i + l;
647
648 return 0;
649 }
650
651
652 /***********************************************************************************/
653 /* Function: location @42_30i@ */
654 /* Description: */
655 /* Return latitude, longitude, and standard longitude for a given location. */
656 /* */
657 /* Parameters: */
658 /* double *latitude: Latitude (North positive). */
659 /* double *longitude: Longitude (West positive). */
660 /* double *long_std: Standard longitude (West positive). */
661 /* char *location: String identifying the location. */
662 /* */
663 /* Return value: */
664 /* 0 if o.k., <0 if error. */
665 /* */
666 /* Example: */
667 /* Files: */
668 /* Known bugs: */
669 /* Author: */
670 /* @i42_30@ */
671 /***********************************************************************************/
672
673 int location (char *locstr, double *latitude, double *longitude, double *long_std)
674 {
675
676 if (!strcasecmp(locstr, "ifu")) {
677 *latitude = 47.48;
678 *longitude = -11.07;
679 *long_std = -15.00;
680
681 return 0;
682 }
683
684 if (!strcasecmp(locstr, "dlrop")) {
685 *latitude = 48.088;
686 *longitude = -11.280;
687 *long_std = -15.000;
688
689 return 0;
690 }
691
692
693 return -1;
694 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22