/[ascend]/trunk/models/johnpye/nrel/spa.c
ViewVC logotype

Annotation of /trunk/models/johnpye/nrel/spa.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2583 - (hide annotations) (download) (as text)
Tue Apr 10 06:02:09 2012 UTC (10 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 35006 byte(s)
Note about code source.
1 jpye 2583 /*
2     THIS CODE COMES FROM APPENDIX A.6 OF REPORT NREL/TP-560-34302
3     WHICH WAS OBTAINED FROM
4     http://www.osti.gov/bridge/servlets/purl/15003974-iP3z6k/native/15003974.PDF
5 jpye 2582
6 jpye 2583 THE CODE HAS BEEN MODIFIED TO
7     (1) convert local functions to 'static' to reduce symbol exports in object code
8     (2) add functionality to allow Julian Day as input, instead of datatime fields.
9     */
10    
11 jpye 2582 /////////////////////////////////////////////
12     // Solar Position Algorithm (SPA) //
13     // for //
14     // Solar Radiation Application //
15     // //
16     // May 12, 2003 //
17     // //
18     // Filename: SPA.C //
19     // //
20     // Afshin Michael Andreas //
21     // Afshin.Andreas@NREL.gov (303)384-6383 //
22     // //
23     // Measurement & Instrumentation Team //
24     // Solar Radiation Research Laboratory //
25     // National Renewable Energy Laboratory //
26     // 1617 Cole Blvd, Golden, CO 80401 //
27     /////////////////////////////////////////////
28    
29     /////////////////////////////////////////////
30     // See the SPA.H header file for usage //
31     // //
32     // This code is based on the NREL //
33     // technical report "Solar Position //
34     // Algorithm for Solar Radiation //
35     // Application" by I. Reda & A. Andreas //
36     /////////////////////////////////////////////
37    
38     ///////////////////////////////////////////////////////////////////////////////////////////////
39     //
40 jpye 2583 // NOTICE
41 jpye 2582 //
42 jpye 2583 //This solar position algorithm for solar radiation applications (the "data") was produced by
43     //the National Renewable Energy Laboratory ("NREL"), which is operated by the Midwest Research
44     //Institute ("MRI") under Contract No. DE-AC36-99-GO10337 with the U.S. Department of Energy
45     //(the "Government").
46 jpye 2582 //
47 jpye 2583 //Reference herein, directly or indirectly to any specific commercial product, process, or
48     //service by trade name, trademark, manufacturer, or otherwise, does not constitute or imply
49     //its endorsement, recommendation, or favoring by the Government, MRI or NREL.
50 jpye 2582 //
51 jpye 2583 //THESE DATA ARE PROVIDED "AS IS" AND NEITHER THE GOVERNMENT, MRI, NREL NOR ANY OF THEIR
52     //EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, INCLUDING THE WARRANTIES OF MERCHANTABILITY
53     //AND FITNESS FOR A PARTICULAR PURPOSE, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE
54     //ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY SUCH INFORMATION DISCLOSED IN THE ALGORITHM, OR
55     //OF ANY APPARATUS, PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE
56     //PRIVATELY OWNED RIGHTS.
57 jpye 2582 //
58     ///////////////////////////////////////////////////////////////////////////////////////////////
59    
60 jpye 2583
61 jpye 2582 ///////////////////////////////////////////////////////////////////////////////////////////////
62     // Revised 27-FEB-2004 Andreas
63     // Added bounds check on inputs and return value for spa_calculate().
64     // Revised 10-MAY-2004 Andreas
65     // Changed temperature bound check minimum from -273.15 to -273 degrees C.
66     // Revised 17-JUN-2004 Andreas
67     // Corrected a problem that caused a bogus sunrise/set/transit on the equinox.
68     // Revised 18-JUN-2004 Andreas
69     // Added a "function" input variable that allows the selecting of desired outputs.
70     // Revised 21-JUN-2004 Andreas
71     // Added 3 new intermediate output values to SPA structure (srha, ssha, & sta).
72     // Revised 23-JUN-2004 Andreas
73     // Enumerations for "function" were renamed and 2 were added.
74     // Prevented bound checks on inputs that are not used (based on function).
75     // Revised 01-SEP-2004 Andreas
76     // Changed a local variable from integer to double.
77     // Revised 12-JUL-2005 Andreas
78     // Put a limit on the EOT calculation, so that the result is between -20 and 20.
79     // Revised 26-OCT-2005 Andreas
80     // Set the atmos. refraction correction to zero, when sun is below horizon.
81     // Made atmos_refract input a requirement for all "functions".
82     // Changed atmos_refract bound check from +/- 10 to +/- 5 degrees.
83     // Revised 07-NOV-2006 Andreas
84     // Corrected 3 earth periodic terms in the L_TERMS array.
85     // Corrected 2 earth periodic terms in the R_TERMS array.
86     // Revised 10-NOV-2006 Andreas
87     // Corrected a constant used to calculate topocentric sun declination.
88     // Put a limit on observer hour angle, so result is between 0 and 360.
89     // Revised 13-NOV-2006 Andreas
90     // Corrected calculation of topocentric sun declination.
91     // Converted all floating point inputs in spa structure to doubles.
92     // Revised 27-FEB-2007 Andreas
93     // Minor correction made as to when atmos. refraction correction is set to zero.
94     // Revised 21-JAN-2008 Andreas
95     // Minor change to two variable declarations.
96     // Revised 12-JAN-2009 Andreas
97     // Changed timezone bound check from +/-12 to +/-18 hours.
98     // Revised 14-JAN-2009 Andreas
99     // Corrected a constant used to calculate ecliptic mean obliquity.
100     ///////////////////////////////////////////////////////////////////////////////////////////////
101    
102     #include <math.h>
103     #include "spa.h"
104    
105     #define PI 3.1415926535897932384626433832795028841971
106     #define SUN_RADIUS 0.26667
107    
108     #define L_COUNT 6
109     #define B_COUNT 2
110     #define R_COUNT 5
111     #define Y_COUNT 63
112    
113     #define L_MAX_SUBCOUNT 64
114     #define B_MAX_SUBCOUNT 5
115     #define R_MAX_SUBCOUNT 40
116    
117     enum {TERM_A, TERM_B, TERM_C, TERM_COUNT};
118     enum {TERM_X0, TERM_X1, TERM_X2, TERM_X3, TERM_X4, TERM_X_COUNT};
119     enum {TERM_PSI_A, TERM_PSI_B, TERM_EPS_C, TERM_EPS_D, TERM_PE_COUNT};
120     enum {JD_MINUS, JD_ZERO, JD_PLUS, JD_COUNT};
121     enum {SUN_TRANSIT, SUN_RISE, SUN_SET, SUN_COUNT};
122    
123     #define TERM_Y_COUNT TERM_X_COUNT
124    
125     const int l_subcount[L_COUNT] = {64,34,20,7,3,1};
126     const int b_subcount[B_COUNT] = {5,2};
127     const int r_subcount[R_COUNT] = {40,10,6,2,1};
128    
129     ///////////////////////////////////////////////////
130     /// Earth Periodic Terms
131     ///////////////////////////////////////////////////
132     const double L_TERMS[L_COUNT][L_MAX_SUBCOUNT][TERM_COUNT]=
133     {
134     {
135     {175347046.0,0,0},
136     {3341656.0,4.6692568,6283.07585},
137     {34894.0,4.6261,12566.1517},
138     {3497.0,2.7441,5753.3849},
139     {3418.0,2.8289,3.5231},
140     {3136.0,3.6277,77713.7715},
141     {2676.0,4.4181,7860.4194},
142     {2343.0,6.1352,3930.2097},
143     {1324.0,0.7425,11506.7698},
144     {1273.0,2.0371,529.691},
145     {1199.0,1.1096,1577.3435},
146     {990,5.233,5884.927},
147     {902,2.045,26.298},
148     {857,3.508,398.149},
149     {780,1.179,5223.694},
150     {753,2.533,5507.553},
151     {505,4.583,18849.228},
152     {492,4.205,775.523},
153     {357,2.92,0.067},
154     {317,5.849,11790.629},
155     {284,1.899,796.298},
156     {271,0.315,10977.079},
157     {243,0.345,5486.778},
158     {206,4.806,2544.314},
159     {205,1.869,5573.143},
160     {202,2.458,6069.777},
161     {156,0.833,213.299},
162     {132,3.411,2942.463},
163     {126,1.083,20.775},
164     {115,0.645,0.98},
165     {103,0.636,4694.003},
166     {102,0.976,15720.839},
167     {102,4.267,7.114},
168     {99,6.21,2146.17},
169     {98,0.68,155.42},
170     {86,5.98,161000.69},
171     {85,1.3,6275.96},
172     {85,3.67,71430.7},
173     {80,1.81,17260.15},
174     {79,3.04,12036.46},
175     {75,1.76,5088.63},
176     {74,3.5,3154.69},
177     {74,4.68,801.82},
178     {70,0.83,9437.76},
179     {62,3.98,8827.39},
180     {61,1.82,7084.9},
181     {57,2.78,6286.6},
182     {56,4.39,14143.5},
183     {56,3.47,6279.55},
184     {52,0.19,12139.55},
185     {52,1.33,1748.02},
186     {51,0.28,5856.48},
187     {49,0.49,1194.45},
188     {41,5.37,8429.24},
189     {41,2.4,19651.05},
190     {39,6.17,10447.39},
191     {37,6.04,10213.29},
192     {37,2.57,1059.38},
193     {36,1.71,2352.87},
194     {36,1.78,6812.77},
195     {33,0.59,17789.85},
196     {30,0.44,83996.85},
197     {30,2.74,1349.87},
198     {25,3.16,4690.48}
199     },
200     {
201     {628331966747.0,0,0},
202     {206059.0,2.678235,6283.07585},
203     {4303.0,2.6351,12566.1517},
204     {425.0,1.59,3.523},
205     {119.0,5.796,26.298},
206     {109.0,2.966,1577.344},
207     {93,2.59,18849.23},
208     {72,1.14,529.69},
209     {68,1.87,398.15},
210     {67,4.41,5507.55},
211     {59,2.89,5223.69},
212     {56,2.17,155.42},
213     {45,0.4,796.3},
214     {36,0.47,775.52},
215     {29,2.65,7.11},
216     {21,5.34,0.98},
217     {19,1.85,5486.78},
218     {19,4.97,213.3},
219     {17,2.99,6275.96},
220     {16,0.03,2544.31},
221     {16,1.43,2146.17},
222     {15,1.21,10977.08},
223     {12,2.83,1748.02},
224     {12,3.26,5088.63},
225     {12,5.27,1194.45},
226     {12,2.08,4694},
227     {11,0.77,553.57},
228     {10,1.3,6286.6},
229     {10,4.24,1349.87},
230     {9,2.7,242.73},
231     {9,5.64,951.72},
232     {8,5.3,2352.87},
233     {6,2.65,9437.76},
234     {6,4.67,4690.48}
235     },
236     {
237     {52919.0,0,0},
238     {8720.0,1.0721,6283.0758},
239     {309.0,0.867,12566.152},
240     {27,0.05,3.52},
241     {16,5.19,26.3},
242     {16,3.68,155.42},
243     {10,0.76,18849.23},
244     {9,2.06,77713.77},
245     {7,0.83,775.52},
246     {5,4.66,1577.34},
247     {4,1.03,7.11},
248     {4,3.44,5573.14},
249     {3,5.14,796.3},
250     {3,6.05,5507.55},
251     {3,1.19,242.73},
252     {3,6.12,529.69},
253     {3,0.31,398.15},
254     {3,2.28,553.57},
255     {2,4.38,5223.69},
256     {2,3.75,0.98}
257     },
258     {
259     {289.0,5.844,6283.076},
260     {35,0,0},
261     {17,5.49,12566.15},
262     {3,5.2,155.42},
263     {1,4.72,3.52},
264     {1,5.3,18849.23},
265     {1,5.97,242.73}
266     },
267     {
268     {114.0,3.142,0},
269     {8,4.13,6283.08},
270     {1,3.84,12566.15}
271     },
272     {
273     {1,3.14,0}
274     }
275     };
276    
277     const double B_TERMS[B_COUNT][B_MAX_SUBCOUNT][TERM_COUNT]=
278     {
279     {
280     {280.0,3.199,84334.662},
281     {102.0,5.422,5507.553},
282     {80,3.88,5223.69},
283     {44,3.7,2352.87},
284     {32,4,1577.34}
285     },
286     {
287     {9,3.9,5507.55},
288     {6,1.73,5223.69}
289     }
290     };
291    
292     const double R_TERMS[R_COUNT][R_MAX_SUBCOUNT][TERM_COUNT]=
293     {
294     {
295     {100013989.0,0,0},
296     {1670700.0,3.0984635,6283.07585},
297     {13956.0,3.05525,12566.1517},
298     {3084.0,5.1985,77713.7715},
299     {1628.0,1.1739,5753.3849},
300     {1576.0,2.8469,7860.4194},
301     {925.0,5.453,11506.77},
302     {542.0,4.564,3930.21},
303     {472.0,3.661,5884.927},
304     {346.0,0.964,5507.553},
305     {329.0,5.9,5223.694},
306     {307.0,0.299,5573.143},
307     {243.0,4.273,11790.629},
308     {212.0,5.847,1577.344},
309     {186.0,5.022,10977.079},
310     {175.0,3.012,18849.228},
311     {110.0,5.055,5486.778},
312     {98,0.89,6069.78},
313     {86,5.69,15720.84},
314     {86,1.27,161000.69},
315     {65,0.27,17260.15},
316     {63,0.92,529.69},
317     {57,2.01,83996.85},
318     {56,5.24,71430.7},
319     {49,3.25,2544.31},
320     {47,2.58,775.52},
321     {45,5.54,9437.76},
322     {43,6.01,6275.96},
323     {39,5.36,4694},
324     {38,2.39,8827.39},
325     {37,0.83,19651.05},
326     {37,4.9,12139.55},
327     {36,1.67,12036.46},
328     {35,1.84,2942.46},
329     {33,0.24,7084.9},
330     {32,0.18,5088.63},
331     {32,1.78,398.15},
332     {28,1.21,6286.6},
333     {28,1.9,6279.55},
334     {26,4.59,10447.39}
335     },
336     {
337     {103019.0,1.10749,6283.07585},
338     {1721.0,1.0644,12566.1517},
339     {702.0,3.142,0},
340     {32,1.02,18849.23},
341     {31,2.84,5507.55},
342     {25,1.32,5223.69},
343     {18,1.42,1577.34},
344     {10,5.91,10977.08},
345     {9,1.42,6275.96},
346     {9,0.27,5486.78}
347     },
348     {
349     {4359.0,5.7846,6283.0758},
350     {124.0,5.579,12566.152},
351     {12,3.14,0},
352     {9,3.63,77713.77},
353     {6,1.87,5573.14},
354     {3,5.47,18849.23}
355     },
356     {
357     {145.0,4.273,6283.076},
358     {7,3.92,12566.15}
359     },
360     {
361     {4,2.56,6283.08}
362     }
363     };
364    
365     ////////////////////////////////////////////////////////////////
366     /// Periodic Terms for the nutation in longitude and obliquity
367     ////////////////////////////////////////////////////////////////
368    
369     const int Y_TERMS[Y_COUNT][TERM_Y_COUNT]=
370     {
371     {0,0,0,0,1},
372     {-2,0,0,2,2},
373     {0,0,0,2,2},
374     {0,0,0,0,2},
375     {0,1,0,0,0},
376     {0,0,1,0,0},
377     {-2,1,0,2,2},
378     {0,0,0,2,1},
379     {0,0,1,2,2},
380     {-2,-1,0,2,2},
381     {-2,0,1,0,0},
382     {-2,0,0,2,1},
383     {0,0,-1,2,2},
384     {2,0,0,0,0},
385     {0,0,1,0,1},
386     {2,0,-1,2,2},
387     {0,0,-1,0,1},
388     {0,0,1,2,1},
389     {-2,0,2,0,0},
390     {0,0,-2,2,1},
391     {2,0,0,2,2},
392     {0,0,2,2,2},
393     {0,0,2,0,0},
394     {-2,0,1,2,2},
395     {0,0,0,2,0},
396     {-2,0,0,2,0},
397     {0,0,-1,2,1},
398     {0,2,0,0,0},
399     {2,0,-1,0,1},
400     {-2,2,0,2,2},
401     {0,1,0,0,1},
402     {-2,0,1,0,1},
403     {0,-1,0,0,1},
404     {0,0,2,-2,0},
405     {2,0,-1,2,1},
406     {2,0,1,2,2},
407     {0,1,0,2,2},
408     {-2,1,1,0,0},
409     {0,-1,0,2,2},
410     {2,0,0,2,1},
411     {2,0,1,0,0},
412     {-2,0,2,2,2},
413     {-2,0,1,2,1},
414     {2,0,-2,0,1},
415     {2,0,0,0,1},
416     {0,-1,1,0,0},
417     {-2,-1,0,2,1},
418     {-2,0,0,0,1},
419     {0,0,2,2,1},
420     {-2,0,2,0,1},
421     {-2,1,0,2,1},
422     {0,0,1,-2,0},
423     {-1,0,1,0,0},
424     {-2,1,0,0,0},
425     {1,0,0,0,0},
426     {0,0,1,2,0},
427     {0,0,-2,2,2},
428     {-1,-1,1,0,0},
429     {0,1,1,0,0},
430     {0,-1,1,2,2},
431     {2,-1,-1,2,2},
432     {0,0,3,2,2},
433     {2,-1,0,2,2},
434     };
435    
436     const double PE_TERMS[Y_COUNT][TERM_PE_COUNT]={
437     {-171996,-174.2,92025,8.9},
438     {-13187,-1.6,5736,-3.1},
439     {-2274,-0.2,977,-0.5},
440     {2062,0.2,-895,0.5},
441     {1426,-3.4,54,-0.1},
442     {712,0.1,-7,0},
443     {-517,1.2,224,-0.6},
444     {-386,-0.4,200,0},
445     {-301,0,129,-0.1},
446     {217,-0.5,-95,0.3},
447     {-158,0,0,0},
448     {129,0.1,-70,0},
449     {123,0,-53,0},
450     {63,0,0,0},
451     {63,0.1,-33,0},
452     {-59,0,26,0},
453     {-58,-0.1,32,0},
454     {-51,0,27,0},
455     {48,0,0,0},
456     {46,0,-24,0},
457     {-38,0,16,0},
458     {-31,0,13,0},
459     {29,0,0,0},
460     {29,0,-12,0},
461     {26,0,0,0},
462     {-22,0,0,0},
463     {21,0,-10,0},
464     {17,-0.1,0,0},
465     {16,0,-8,0},
466     {-16,0.1,7,0},
467     {-15,0,9,0},
468     {-13,0,7,0},
469     {-12,0,6,0},
470     {11,0,0,0},
471     {-10,0,5,0},
472     {-8,0,3,0},
473     {7,0,-3,0},
474     {-7,0,0,0},
475     {-7,0,3,0},
476     {-7,0,3,0},
477     {6,0,0,0},
478     {6,0,-3,0},
479     {6,0,-3,0},
480     {-6,0,3,0},
481     {-6,0,3,0},
482     {5,0,0,0},
483     {-5,0,3,0},
484     {-5,0,3,0},
485     {-5,0,3,0},
486     {4,0,0,0},
487     {4,0,0,0},
488     {4,0,0,0},
489     {-4,0,0,0},
490     {-4,0,0,0},
491     {-4,0,0,0},
492     {3,0,0,0},
493     {-3,0,0,0},
494     {-3,0,0,0},
495     {-3,0,0,0},
496     {-3,0,0,0},
497     {-3,0,0,0},
498     {-3,0,0,0},
499     {-3,0,0,0},
500     };
501    
502     ///////////////////////////////////////////////
503    
504     double rad2deg(double radians)
505     {
506     return (180.0/PI)*radians;
507     }
508    
509     double deg2rad(double degrees)
510     {
511     return (PI/180.0)*degrees;
512     }
513    
514     double limit_degrees(double degrees)
515     {
516     double limited;
517    
518     degrees /= 360.0;
519     limited = 360.0*(degrees-floor(degrees));
520     if (limited < 0) limited += 360.0;
521    
522     return limited;
523     }
524    
525     double limit_degrees180pm(double degrees)
526     {
527     double limited;
528    
529     degrees /= 360.0;
530     limited = 360.0*(degrees-floor(degrees));
531     if (limited < -180.0) limited += 360.0;
532     else if (limited > 180.0) limited -= 360.0;
533    
534     return limited;
535     }
536    
537     double limit_degrees180(double degrees)
538     {
539     double limited;
540    
541     degrees /= 180.0;
542     limited = 180.0*(degrees-floor(degrees));
543     if (limited < 0) limited += 180.0;
544    
545     return limited;
546     }
547    
548     double limit_zero2one(double value)
549     {
550     double limited;
551    
552     limited = value - floor(value);
553     if (limited < 0) limited += 1.0;
554    
555     return limited;
556     }
557    
558     double limit_minutes(double minutes)
559     {
560     double limited=minutes;
561    
562     if (limited < -20.0) limited += 1440.0;
563     else if (limited > 20.0) limited -= 1440.0;
564    
565     return limited;
566     }
567    
568     double dayfrac_to_local_hr(double dayfrac, double timezone)
569     {
570     return 24.0*limit_zero2one(dayfrac + timezone/24.0);
571     }
572    
573     double third_order_polynomial(double a, double b, double c, double d, double x)
574     {
575     return ((a*x + b)*x + c)*x + d;
576     }
577    
578     ///////////////////////////////////////////////////////////////////////////////////////////////
579     int validate_inputs(spa_data *spa)
580     {
581     if ((spa->year < -2000) || (spa->year > 6000)) return 1;
582     if ((spa->month < 1 ) || (spa->month > 12 )) return 2;
583     if ((spa->day < 1 ) || (spa->day > 31 )) return 3;
584     if ((spa->hour < 0 ) || (spa->hour > 24 )) return 4;
585     if ((spa->minute < 0 ) || (spa->minute > 59 )) return 5;
586     if ((spa->second < 0 ) || (spa->second > 59 )) return 6;
587     if ((spa->pressure < 0 ) || (spa->pressure > 5000)) return 12;
588     if ((spa->temperature <= -273) || (spa->temperature > 6000)) return 13;
589     if ((spa->hour == 24 ) && (spa->minute > 0 )) return 5;
590     if ((spa->hour == 24 ) && (spa->second > 0 )) return 6;
591    
592     if (fabs(spa->delta_t) > 8000 ) return 7;
593     if (fabs(spa->timezone) > 18 ) return 8;
594     if (fabs(spa->longitude) > 180 ) return 9;
595     if (fabs(spa->latitude) > 90 ) return 10;
596     if (fabs(spa->atmos_refract) > 5 ) return 16;
597     if ( spa->elevation < -6500000) return 11;
598    
599     if ((spa->function == SPA_ZA_INC) || (spa->function == SPA_ALL))
600     {
601     if (fabs(spa->slope) > 360) return 14;
602     if (fabs(spa->azm_rotation) > 360) return 15;
603     }
604    
605     return 0;
606     }
607     ///////////////////////////////////////////////////////////////////////////////////////////////
608     double julian_day (int year, int month, int day, int hour, int minute, int second, double tz)
609     {
610     double day_decimal, julian_day, a;
611    
612     day_decimal = day + (hour - tz + (minute + second/60.0)/60.0)/24.0;
613    
614     if (month < 3) {
615     month += 12;
616     year--;
617     }
618    
619     julian_day = floor(365.25*(year+4716.0)) + floor(30.6001*(month+1)) + day_decimal - 1524.5;
620    
621     if (julian_day > 2299160.0) {
622     a = floor(year/100);
623     julian_day += (2 - a + floor(a/4));
624     }
625    
626     return julian_day;
627     }
628    
629     double julian_century(double jd)
630     {
631     return (jd-2451545.0)/36525.0;
632     }
633    
634     double julian_ephemeris_day(double jd, double delta_t)
635     {
636     return jd+delta_t/86400.0;
637     }
638    
639     double julian_ephemeris_century(double jde)
640     {
641     return (jde - 2451545.0)/36525.0;
642     }
643    
644     double julian_ephemeris_millennium(double jce)
645     {
646     return (jce/10.0);
647     }
648    
649     double earth_periodic_term_summation(const double terms[][TERM_COUNT], int count, double jme)
650     {
651     int i;
652     double sum=0;
653    
654     for (i = 0; i < count; i++)
655     sum += terms[i][TERM_A]*cos(terms[i][TERM_B]+terms[i][TERM_C]*jme);
656    
657     return sum;
658     }
659    
660     double earth_values(double term_sum[], int count, double jme)
661     {
662     int i;
663     double sum=0;
664    
665     for (i = 0; i < count; i++)
666     sum += term_sum[i]*pow(jme, i);
667    
668     sum /= 1.0e8;
669    
670     return sum;
671     }
672    
673     double earth_heliocentric_longitude(double jme)
674     {
675     double sum[L_COUNT];
676     int i;
677    
678     for (i = 0; i < L_COUNT; i++)
679     sum[i] = earth_periodic_term_summation(L_TERMS[i], l_subcount[i], jme);
680    
681     return limit_degrees(rad2deg(earth_values(sum, L_COUNT, jme)));
682    
683     }
684    
685     double earth_heliocentric_latitude(double jme)
686     {
687     double sum[B_COUNT];
688     int i;
689    
690     for (i = 0; i < B_COUNT; i++)
691     sum[i] = earth_periodic_term_summation(B_TERMS[i], b_subcount[i], jme);
692    
693     return rad2deg(earth_values(sum, B_COUNT, jme));
694    
695     }
696    
697     double earth_radius_vector(double jme)
698     {
699     double sum[R_COUNT];
700     int i;
701    
702     for (i = 0; i < R_COUNT; i++)
703     sum[i] = earth_periodic_term_summation(R_TERMS[i], r_subcount[i], jme);
704    
705     return earth_values(sum, R_COUNT, jme);
706    
707     }
708    
709     double geocentric_longitude(double l)
710     {
711     double theta = l + 180.0;
712    
713     if (theta >= 360.0) theta -= 360.0;
714    
715     return theta;
716     }
717    
718     double geocentric_latitude(double b)
719     {
720     return -b;
721     }
722    
723     double mean_elongation_moon_sun(double jce)
724     {
725     return third_order_polynomial(1.0/189474.0, -0.0019142, 445267.11148, 297.85036, jce);
726     }
727    
728     double mean_anomaly_sun(double jce)
729     {
730     return third_order_polynomial(-1.0/300000.0, -0.0001603, 35999.05034, 357.52772, jce);
731     }
732    
733     double mean_anomaly_moon(double jce)
734     {
735     return third_order_polynomial(1.0/56250.0, 0.0086972, 477198.867398, 134.96298, jce);
736     }
737    
738     double argument_latitude_moon(double jce)
739     {
740     return third_order_polynomial(1.0/327270.0, -0.0036825, 483202.017538, 93.27191, jce);
741     }
742    
743     double ascending_longitude_moon(double jce)
744     {
745     return third_order_polynomial(1.0/450000.0, 0.0020708, -1934.136261, 125.04452, jce);
746     }
747    
748     double xy_term_summation(int i, double x[TERM_X_COUNT])
749     {
750     int j;
751     double sum=0;
752    
753     for (j = 0; j < TERM_Y_COUNT; j++)
754     sum += x[j]*Y_TERMS[i][j];
755    
756     return sum;
757     }
758    
759     void nutation_longitude_and_obliquity(double jce, double x[TERM_X_COUNT], double *del_psi,
760     double *del_epsilon)
761     {
762     int i;
763     double xy_term_sum, sum_psi=0, sum_epsilon=0;
764    
765     for (i = 0; i < Y_COUNT; i++) {
766     xy_term_sum = deg2rad(xy_term_summation(i, x));
767     sum_psi += (PE_TERMS[i][TERM_PSI_A] + jce*PE_TERMS[i][TERM_PSI_B])*sin(xy_term_sum);
768     sum_epsilon += (PE_TERMS[i][TERM_EPS_C] + jce*PE_TERMS[i][TERM_EPS_D])*cos(xy_term_sum);
769     }
770    
771     *del_psi = sum_psi / 36000000.0;
772     *del_epsilon = sum_epsilon / 36000000.0;
773     }
774    
775     double ecliptic_mean_obliquity(double jme)
776     {
777     double u = jme/10.0;
778    
779     return 84381.448 + u*(-4680.93 + u*(-1.55 + u*(1999.25 + u*(-51.38 + u*(-249.67 +
780     u*( -39.05 + u*( 7.12 + u*( 27.87 + u*( 5.79 + u*2.45)))))))));
781     }
782    
783     double ecliptic_true_obliquity(double delta_epsilon, double epsilon0)
784     {
785     return delta_epsilon + epsilon0/3600.0;
786     }
787    
788     double aberration_correction(double r)
789     {
790     return -20.4898 / (3600.0*r);
791     }
792    
793     double apparent_sun_longitude(double theta, double delta_psi, double delta_tau)
794     {
795     return theta + delta_psi + delta_tau;
796     }
797    
798     double greenwich_mean_sidereal_time (double jd, double jc)
799     {
800     return limit_degrees(280.46061837 + 360.98564736629 * (jd - 2451545.0) +
801     jc*jc*(0.000387933 - jc/38710000.0));
802     }
803    
804     double greenwich_sidereal_time (double nu0, double delta_psi, double epsilon)
805     {
806     return nu0 + delta_psi*cos(deg2rad(epsilon));
807     }
808    
809     double geocentric_sun_right_ascension(double lamda, double epsilon, double beta)
810     {
811     double lamda_rad = deg2rad(lamda);
812     double epsilon_rad = deg2rad(epsilon);
813    
814     return limit_degrees(rad2deg(atan2(sin(lamda_rad)*cos(epsilon_rad) -
815     tan(deg2rad(beta))*sin(epsilon_rad), cos(lamda_rad))));
816     }
817    
818     double geocentric_sun_declination(double beta, double epsilon, double lamda)
819     {
820     double beta_rad = deg2rad(beta);
821     double epsilon_rad = deg2rad(epsilon);
822    
823     return rad2deg(asin(sin(beta_rad)*cos(epsilon_rad) +
824     cos(beta_rad)*sin(epsilon_rad)*sin(deg2rad(lamda))));
825     }
826    
827     double observer_hour_angle(double nu, double longitude, double alpha_deg)
828     {
829     return limit_degrees(nu + longitude - alpha_deg);
830     }
831    
832     double sun_equatorial_horizontal_parallax(double r)
833     {
834     return 8.794 / (3600.0 * r);
835     }
836    
837     void sun_right_ascension_parallax_and_topocentric_dec(double latitude, double elevation,
838     double xi, double h, double delta, double *delta_alpha, double *delta_prime)
839     {
840     double delta_alpha_rad;
841     double lat_rad = deg2rad(latitude);
842     double xi_rad = deg2rad(xi);
843     double h_rad = deg2rad(h);
844     double delta_rad = deg2rad(delta);
845     double u = atan(0.99664719 * tan(lat_rad));
846     double y = 0.99664719 * sin(u) + elevation*sin(lat_rad)/6378140.0;
847     double x = cos(u) + elevation*cos(lat_rad)/6378140.0;
848    
849     delta_alpha_rad = atan2( - x*sin(xi_rad) *sin(h_rad),
850     cos(delta_rad) - x*sin(xi_rad) *cos(h_rad));
851    
852     *delta_prime = rad2deg(atan2((sin(delta_rad) - y*sin(xi_rad))*cos(delta_alpha_rad),
853     cos(delta_rad) - x*sin(xi_rad) *cos(h_rad)));
854    
855     *delta_alpha = rad2deg(delta_alpha_rad);
856     }
857    
858     double topocentric_sun_right_ascension(double alpha_deg, double delta_alpha)
859     {
860     return alpha_deg + delta_alpha;
861     }
862    
863     double topocentric_local_hour_angle(double h, double delta_alpha)
864     {
865     return h - delta_alpha;
866     }
867    
868     double topocentric_elevation_angle(double latitude, double delta_prime, double h_prime)
869     {
870     double lat_rad = deg2rad(latitude);
871     double delta_prime_rad = deg2rad(delta_prime);
872    
873     return rad2deg(asin(sin(lat_rad)*sin(delta_prime_rad) +
874     cos(lat_rad)*cos(delta_prime_rad) * cos(deg2rad(h_prime))));
875     }
876    
877     double atmospheric_refraction_correction(double pressure, double temperature,
878     double atmos_refract, double e0)
879     {
880     double del_e = 0;
881    
882     if (e0 >= -1*(SUN_RADIUS + atmos_refract))
883     del_e = (pressure / 1010.0) * (283.0 / (273.0 + temperature)) *
884     1.02 / (60.0 * tan(deg2rad(e0 + 10.3/(e0 + 5.11))));
885    
886     return del_e;
887     }
888    
889     double topocentric_elevation_angle_corrected(double e0, double delta_e)
890     {
891     return e0 + delta_e;
892     }
893    
894     double topocentric_zenith_angle(double e)
895     {
896     return 90.0 - e;
897     }
898    
899     double topocentric_azimuth_angle_neg180_180(double h_prime, double latitude, double delta_prime)
900     {
901     double h_prime_rad = deg2rad(h_prime);
902     double lat_rad = deg2rad(latitude);
903    
904     return rad2deg(atan2(sin(h_prime_rad),
905     cos(h_prime_rad)*sin(lat_rad) - tan(deg2rad(delta_prime))*cos(lat_rad)));
906     }
907    
908     double topocentric_azimuth_angle_zero_360(double azimuth180)
909     {
910     return azimuth180 + 180.0;
911     }
912    
913     double surface_incidence_angle(double zenith, double azimuth180, double azm_rotation,
914     double slope)
915     {
916     double zenith_rad = deg2rad(zenith);
917     double slope_rad = deg2rad(slope);
918    
919     return rad2deg(acos(cos(zenith_rad)*cos(slope_rad) +
920     sin(slope_rad )*sin(zenith_rad) * cos(deg2rad(azimuth180 - azm_rotation))));
921     }
922    
923     double sun_mean_longitude(double jme)
924     {
925     return limit_degrees(280.4664567 + jme*(360007.6982779 + jme*(0.03032028 +
926     jme*(1/49931.0 + jme*(-1/15300.0 + jme*(-1/2000000.0))))));
927     }
928    
929     double eot(double m, double alpha, double del_psi, double epsilon)
930     {
931     return limit_minutes(4.0*(m - 0.0057183 - alpha + del_psi*cos(deg2rad(epsilon))));
932     }
933    
934     double approx_sun_transit_time(double alpha_zero, double longitude, double nu)
935     {
936     return (alpha_zero - longitude - nu) / 360.0;
937     }
938    
939     double sun_hour_angle_at_rise_set(double latitude, double delta_zero, double h0_prime)
940     {
941     double h0 = -99999;
942     double latitude_rad = deg2rad(latitude);
943     double delta_zero_rad = deg2rad(delta_zero);
944     double argument = (sin(deg2rad(h0_prime)) - sin(latitude_rad)*sin(delta_zero_rad)) /
945     (cos(latitude_rad)*cos(delta_zero_rad));
946    
947     if (fabs(argument) <= 1) h0 = limit_degrees180(rad2deg(acos(argument)));
948    
949     return h0;
950     }
951    
952     void approx_sun_rise_and_set(double *m_rts, double h0)
953     {
954     double h0_dfrac = h0/360.0;
955    
956     m_rts[SUN_RISE] = limit_zero2one(m_rts[SUN_TRANSIT] - h0_dfrac);
957     m_rts[SUN_SET] = limit_zero2one(m_rts[SUN_TRANSIT] + h0_dfrac);
958     m_rts[SUN_TRANSIT] = limit_zero2one(m_rts[SUN_TRANSIT]);
959     }
960    
961     double rts_alpha_delta_prime(double *ad, double n)
962     {
963     double a = ad[JD_ZERO] - ad[JD_MINUS];
964     double b = ad[JD_PLUS] - ad[JD_ZERO];
965    
966     if (fabs(a) >= 2.0) a = limit_zero2one(a);
967     if (fabs(b) >= 2.0) b = limit_zero2one(b);
968    
969     return ad[JD_ZERO] + n * (a + b + (b-a)*n)/2.0;
970     }
971    
972     double rts_sun_altitude(double latitude, double delta_prime, double h_prime)
973     {
974     double latitude_rad = deg2rad(latitude);
975     double delta_prime_rad = deg2rad(delta_prime);
976    
977     return rad2deg(asin(sin(latitude_rad)*sin(delta_prime_rad) +
978     cos(latitude_rad)*cos(delta_prime_rad)*cos(deg2rad(h_prime))));
979     }
980    
981     double sun_rise_and_set(double *m_rts, double *h_rts, double *delta_prime, double latitude,
982     double *h_prime, double h0_prime, int sun)
983     {
984     return m_rts[sun] + (h_rts[sun] - h0_prime) /
985     (360.0*cos(deg2rad(delta_prime[sun]))*cos(deg2rad(latitude))*sin(deg2rad(h_prime[sun])));
986     }
987    
988     ////////////////////////////////////////////////////////////////////////////////////////////////
989     // Calculate required SPA parameters to get the right ascension (alpha) and declination (delta)
990     // Note: JD must be already calculated and in structure
991     ////////////////////////////////////////////////////////////////////////////////////////////////
992     void calculate_geocentric_sun_right_ascension_and_declination(spa_data *spa)
993     {
994     double x[TERM_X_COUNT];
995    
996     spa->jc = julian_century(spa->jd);
997    
998     spa->jde = julian_ephemeris_day(spa->jd, spa->delta_t);
999     spa->jce = julian_ephemeris_century(spa->jde);
1000     spa->jme = julian_ephemeris_millennium(spa->jce);
1001    
1002     spa->l = earth_heliocentric_longitude(spa->jme);
1003     spa->b = earth_heliocentric_latitude(spa->jme);
1004     spa->r = earth_radius_vector(spa->jme);
1005    
1006     spa->theta = geocentric_longitude(spa->l);
1007     spa->beta = geocentric_latitude(spa->b);
1008    
1009     x[TERM_X0] = spa->x0 = mean_elongation_moon_sun(spa->jce);
1010     x[TERM_X1] = spa->x1 = mean_anomaly_sun(spa->jce);
1011     x[TERM_X2] = spa->x2 = mean_anomaly_moon(spa->jce);
1012     x[TERM_X3] = spa->x3 = argument_latitude_moon(spa->jce);
1013     x[TERM_X4] = spa->x4 = ascending_longitude_moon(spa->jce);
1014    
1015     nutation_longitude_and_obliquity(spa->jce, x, &(spa->del_psi), &(spa->del_epsilon));
1016    
1017     spa->epsilon0 = ecliptic_mean_obliquity(spa->jme);
1018     spa->epsilon = ecliptic_true_obliquity(spa->del_epsilon, spa->epsilon0);
1019    
1020     spa->del_tau = aberration_correction(spa->r);
1021     spa->lamda = apparent_sun_longitude(spa->theta, spa->del_psi, spa->del_tau);
1022     spa->nu0 = greenwich_mean_sidereal_time (spa->jd, spa->jc);
1023     spa->nu = greenwich_sidereal_time (spa->nu0, spa->del_psi, spa->epsilon);
1024    
1025     spa->alpha = geocentric_sun_right_ascension(spa->lamda, spa->epsilon, spa->beta);
1026     spa->delta = geocentric_sun_declination(spa->beta, spa->epsilon, spa->lamda);
1027     }
1028    
1029     ////////////////////////////////////////////////////////////////////////
1030     // Calculate Equation of Time (EOT) and Sun Rise, Transit, & Set (RTS)
1031     ////////////////////////////////////////////////////////////////////////
1032    
1033     void calculate_eot_and_sun_rise_transit_set(spa_data *spa)
1034     {
1035     spa_data sun_rts;
1036     double nu, m, h0, n;
1037     double alpha[JD_COUNT], delta[JD_COUNT];
1038     double m_rts[SUN_COUNT], nu_rts[SUN_COUNT], h_rts[SUN_COUNT];
1039     double alpha_prime[SUN_COUNT], delta_prime[SUN_COUNT], h_prime[SUN_COUNT];
1040     double h0_prime = -1*(SUN_RADIUS + spa->atmos_refract);
1041     int i;
1042    
1043     sun_rts = *spa;
1044     m = sun_mean_longitude(spa->jme);
1045     spa->eot = eot(m, spa->alpha, spa->del_psi, spa->epsilon);
1046    
1047     sun_rts.hour = sun_rts.minute = sun_rts.second = 0;
1048     sun_rts.timezone = 0.0;
1049    
1050     sun_rts.jd = julian_day (sun_rts.year, sun_rts.month, sun_rts.day,
1051     sun_rts.hour, sun_rts.minute, sun_rts.second, sun_rts.timezone);
1052    
1053     calculate_geocentric_sun_right_ascension_and_declination(&sun_rts);
1054     nu = sun_rts.nu;
1055    
1056     sun_rts.delta_t = 0;
1057     sun_rts.jd--;
1058     for (i = 0; i < JD_COUNT; i++) {
1059     calculate_geocentric_sun_right_ascension_and_declination(&sun_rts);
1060     alpha[i] = sun_rts.alpha;
1061     delta[i] = sun_rts.delta;
1062     sun_rts.jd++;
1063     }
1064    
1065     m_rts[SUN_TRANSIT] = approx_sun_transit_time(alpha[JD_ZERO], spa->longitude, nu);
1066     h0 = sun_hour_angle_at_rise_set(spa->latitude, delta[JD_ZERO], h0_prime);
1067    
1068     if (h0 >= 0) {
1069    
1070     approx_sun_rise_and_set(m_rts, h0);
1071    
1072     for (i = 0; i < SUN_COUNT; i++) {
1073    
1074     nu_rts[i] = nu + 360.985647*m_rts[i];
1075    
1076     n = m_rts[i] + spa->delta_t/86400.0;
1077     alpha_prime[i] = rts_alpha_delta_prime(alpha, n);
1078     delta_prime[i] = rts_alpha_delta_prime(delta, n);
1079    
1080     h_prime[i] = limit_degrees180pm(nu_rts[i] + spa->longitude - alpha_prime[i]);
1081    
1082     h_rts[i] = rts_sun_altitude(spa->latitude, delta_prime[i], h_prime[i]);
1083     }
1084    
1085     spa->srha = h_prime[SUN_RISE];
1086     spa->ssha = h_prime[SUN_SET];
1087     spa->sta = h_rts[SUN_TRANSIT];
1088    
1089     spa->suntransit = dayfrac_to_local_hr(m_rts[SUN_TRANSIT] - h_prime[SUN_TRANSIT] / 360.0,
1090     spa->timezone);
1091    
1092     spa->sunrise = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime,
1093     spa->latitude, h_prime, h0_prime, SUN_RISE), spa->timezone);
1094    
1095     spa->sunset = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime,
1096     spa->latitude, h_prime, h0_prime, SUN_SET), spa->timezone);
1097    
1098     } else spa->srha= spa->ssha= spa->sta= spa->suntransit= spa->sunrise= spa->sunset= -99999;
1099    
1100     }
1101    
1102     ///////////////////////////////////////////////////////////////////////////////////////////
1103     // Calculate all SPA parameters and put into structure
1104     // Note: All inputs values (listed in header file) must already be in structure
1105     ///////////////////////////////////////////////////////////////////////////////////////////
1106     int spa_calculate(spa_data *spa)
1107     {
1108     int result;
1109    
1110     result = validate_inputs(spa);
1111    
1112     if (result == 0)
1113     {
1114     spa->jd = julian_day (spa->year, spa->month, spa->day,
1115     spa->hour, spa->minute, spa->second, spa->timezone);
1116    
1117     calculate_geocentric_sun_right_ascension_and_declination(spa);
1118    
1119     spa->h = observer_hour_angle(spa->nu, spa->longitude, spa->alpha);
1120     spa->xi = sun_equatorial_horizontal_parallax(spa->r);
1121    
1122     sun_right_ascension_parallax_and_topocentric_dec(spa->latitude, spa->elevation, spa->xi,
1123     spa->h, spa->delta, &(spa->del_alpha), &(spa->delta_prime));
1124    
1125     spa->alpha_prime = topocentric_sun_right_ascension(spa->alpha, spa->del_alpha);
1126     spa->h_prime = topocentric_local_hour_angle(spa->h, spa->del_alpha);
1127    
1128     spa->e0 = topocentric_elevation_angle(spa->latitude, spa->delta_prime, spa->h_prime);
1129     spa->del_e = atmospheric_refraction_correction(spa->pressure, spa->temperature,
1130     spa->atmos_refract, spa->e0);
1131     spa->e = topocentric_elevation_angle_corrected(spa->e0, spa->del_e);
1132    
1133     spa->zenith = topocentric_zenith_angle(spa->e);
1134     spa->azimuth180 = topocentric_azimuth_angle_neg180_180(spa->h_prime, spa->latitude,
1135     spa->delta_prime);
1136     spa->azimuth = topocentric_azimuth_angle_zero_360(spa->azimuth180);
1137    
1138     if ((spa->function == SPA_ZA_INC) || (spa->function == SPA_ALL))
1139     spa->incidence = surface_incidence_angle(spa->zenith, spa->azimuth180,
1140     spa->azm_rotation, spa->slope);
1141    
1142     if ((spa->function == SPA_ZA_RTS) || (spa->function == SPA_ALL))
1143     calculate_eot_and_sun_rise_transit_set(spa);
1144     }
1145    
1146     return result;
1147     }
1148     ///////////////////////////////////////////////////////////////////////////////////////////

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