/[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 2599 - (hide annotations) (download) (as text)
Wed Apr 18 18:05:24 2012 UTC (10 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 36370 byte(s)
Add 'energy_per_area' atom type.
Model solar/incident.a4c shows TMY3 and NREL SPA working together.
Fixed error in bounds check for JD in SPA.
Fix typo in relation_util.c.
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 jpye 2584 static const int l_subcount[L_COUNT] = {64,34,20,7,3,1};
126     static const int b_subcount[B_COUNT] = {5,2};
127     static const int r_subcount[R_COUNT] = {40,10,6,2,1};
128 jpye 2582
129     ///////////////////////////////////////////////////
130     /// Earth Periodic Terms
131     ///////////////////////////////////////////////////
132 jpye 2584 static const double L_TERMS[L_COUNT][L_MAX_SUBCOUNT][TERM_COUNT]=
133 jpye 2582 {
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 jpye 2584 static const double B_TERMS[B_COUNT][B_MAX_SUBCOUNT][TERM_COUNT]=
278 jpye 2582 {
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 jpye 2584 static const double R_TERMS[R_COUNT][R_MAX_SUBCOUNT][TERM_COUNT]=
293 jpye 2582 {
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 jpye 2584 static const int Y_TERMS[Y_COUNT][TERM_Y_COUNT]=
370 jpye 2582 {
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 jpye 2584 static const double PE_TERMS[Y_COUNT][TERM_PE_COUNT]={
437 jpye 2582 {-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 jpye 2584 static double rad2deg(double radians)
505 jpye 2582 {
506     return (180.0/PI)*radians;
507     }
508    
509 jpye 2584 static double deg2rad(double degrees)
510 jpye 2582 {
511     return (PI/180.0)*degrees;
512     }
513    
514 jpye 2584 static double limit_degrees(double degrees)
515 jpye 2582 {
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 jpye 2584 static double limit_degrees180pm(double degrees)
526 jpye 2582 {
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 jpye 2584 static double limit_degrees180(double degrees)
538 jpye 2582 {
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 jpye 2584 static double limit_zero2one(double value)
549 jpye 2582 {
550     double limited;
551    
552     limited = value - floor(value);
553     if (limited < 0) limited += 1.0;
554    
555     return limited;
556     }
557    
558 jpye 2584 static double limit_minutes(double minutes)
559 jpye 2582 {
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 jpye 2584 static double dayfrac_to_local_hr(double dayfrac, double timezone)
569 jpye 2582 {
570     return 24.0*limit_zero2one(dayfrac + timezone/24.0);
571     }
572    
573 jpye 2584 static double third_order_polynomial(double a, double b, double c, double d, double x)
574 jpye 2582 {
575     return ((a*x + b)*x + c)*x + d;
576     }
577    
578     ///////////////////////////////////////////////////////////////////////////////////////////////
579 jpye 2584 static int validate_inputs(spa_data *spa)
580 jpye 2582 {
581 jpye 2584 switch(spa->function){
582     case SPA_ZA_JD:
583     case SPA_ZA_INC_JD:
584     case SPA_ZA_RTS_JD:
585 jpye 2599 if ((spa->jd < 990557.5) || (spa->jd > 3912879.500000)) return 16;
586 jpye 2584 break;
587     default:
588     if ((spa->year < -2000) || (spa->year > 6000)) return 1;
589     if ((spa->month < 1 ) || (spa->month > 12 )) return 2;
590     if ((spa->day < 1 ) || (spa->day > 31 )) return 3;
591     if ((spa->hour < 0 ) || (spa->hour > 24 )) return 4;
592     if ((spa->minute < 0 ) || (spa->minute > 59 )) return 5;
593     if ((spa->second < 0 ) || (spa->second > 59 )) return 6;
594     }
595    
596 jpye 2582 if ((spa->pressure < 0 ) || (spa->pressure > 5000)) return 12;
597     if ((spa->temperature <= -273) || (spa->temperature > 6000)) return 13;
598     if ((spa->hour == 24 ) && (spa->minute > 0 )) return 5;
599     if ((spa->hour == 24 ) && (spa->second > 0 )) return 6;
600    
601     if (fabs(spa->delta_t) > 8000 ) return 7;
602     if (fabs(spa->timezone) > 18 ) return 8;
603     if (fabs(spa->longitude) > 180 ) return 9;
604     if (fabs(spa->latitude) > 90 ) return 10;
605     if (fabs(spa->atmos_refract) > 5 ) return 16;
606     if ( spa->elevation < -6500000) return 11;
607    
608     if ((spa->function == SPA_ZA_INC) || (spa->function == SPA_ALL))
609     {
610     if (fabs(spa->slope) > 360) return 14;
611     if (fabs(spa->azm_rotation) > 360) return 15;
612     }
613    
614     return 0;
615     }
616     ///////////////////////////////////////////////////////////////////////////////////////////////
617 jpye 2591 double julian_day (int year, int month, int day, int hour, int minute, int second, double tz)
618 jpye 2582 {
619     double day_decimal, julian_day, a;
620    
621 jpye 2599 //fprintf(stderr,"Julian day for %d/%d/%d %d:%d:%d (%+0.2f) is ",day,month,year,hour,minute,second,tz);
622    
623 jpye 2582 day_decimal = day + (hour - tz + (minute + second/60.0)/60.0)/24.0;
624    
625     if (month < 3) {
626     month += 12;
627     year--;
628     }
629    
630     julian_day = floor(365.25*(year+4716.0)) + floor(30.6001*(month+1)) + day_decimal - 1524.5;
631    
632     if (julian_day > 2299160.0) {
633     a = floor(year/100);
634     julian_day += (2 - a + floor(a/4));
635     }
636    
637 jpye 2599 //fprintf(stderr,"%f\n",julian_day);
638 jpye 2582 return julian_day;
639     }
640    
641 jpye 2584 static double julian_century(double jd)
642 jpye 2582 {
643     return (jd-2451545.0)/36525.0;
644     }
645    
646 jpye 2584 static double julian_ephemeris_day(double jd, double delta_t)
647 jpye 2582 {
648     return jd+delta_t/86400.0;
649     }
650    
651 jpye 2584 static double julian_ephemeris_century(double jde)
652 jpye 2582 {
653     return (jde - 2451545.0)/36525.0;
654     }
655    
656 jpye 2584 static double julian_ephemeris_millennium(double jce)
657 jpye 2582 {
658     return (jce/10.0);
659     }
660    
661 jpye 2584 static double earth_periodic_term_summation(const double terms[][TERM_COUNT], int count, double jme)
662 jpye 2582 {
663     int i;
664     double sum=0;
665    
666     for (i = 0; i < count; i++)
667     sum += terms[i][TERM_A]*cos(terms[i][TERM_B]+terms[i][TERM_C]*jme);
668    
669     return sum;
670     }
671    
672 jpye 2584 static double earth_values(double term_sum[], int count, double jme)
673 jpye 2582 {
674     int i;
675     double sum=0;
676    
677     for (i = 0; i < count; i++)
678     sum += term_sum[i]*pow(jme, i);
679    
680     sum /= 1.0e8;
681    
682     return sum;
683     }
684    
685 jpye 2584 static double earth_heliocentric_longitude(double jme)
686 jpye 2582 {
687     double sum[L_COUNT];
688     int i;
689    
690     for (i = 0; i < L_COUNT; i++)
691     sum[i] = earth_periodic_term_summation(L_TERMS[i], l_subcount[i], jme);
692    
693     return limit_degrees(rad2deg(earth_values(sum, L_COUNT, jme)));
694    
695     }
696    
697 jpye 2584 static double earth_heliocentric_latitude(double jme)
698 jpye 2582 {
699     double sum[B_COUNT];
700     int i;
701    
702     for (i = 0; i < B_COUNT; i++)
703     sum[i] = earth_periodic_term_summation(B_TERMS[i], b_subcount[i], jme);
704    
705     return rad2deg(earth_values(sum, B_COUNT, jme));
706    
707     }
708    
709 jpye 2584 static double earth_radius_vector(double jme)
710 jpye 2582 {
711     double sum[R_COUNT];
712     int i;
713    
714     for (i = 0; i < R_COUNT; i++)
715     sum[i] = earth_periodic_term_summation(R_TERMS[i], r_subcount[i], jme);
716    
717     return earth_values(sum, R_COUNT, jme);
718    
719     }
720    
721 jpye 2584 static double geocentric_longitude(double l)
722 jpye 2582 {
723     double theta = l + 180.0;
724    
725     if (theta >= 360.0) theta -= 360.0;
726    
727     return theta;
728     }
729    
730 jpye 2584 static double geocentric_latitude(double b)
731 jpye 2582 {
732     return -b;
733     }
734    
735 jpye 2584 static double mean_elongation_moon_sun(double jce)
736 jpye 2582 {
737     return third_order_polynomial(1.0/189474.0, -0.0019142, 445267.11148, 297.85036, jce);
738     }
739    
740 jpye 2584 static double mean_anomaly_sun(double jce)
741 jpye 2582 {
742     return third_order_polynomial(-1.0/300000.0, -0.0001603, 35999.05034, 357.52772, jce);
743     }
744    
745 jpye 2584 static double mean_anomaly_moon(double jce)
746 jpye 2582 {
747     return third_order_polynomial(1.0/56250.0, 0.0086972, 477198.867398, 134.96298, jce);
748     }
749    
750 jpye 2584 static double argument_latitude_moon(double jce)
751 jpye 2582 {
752     return third_order_polynomial(1.0/327270.0, -0.0036825, 483202.017538, 93.27191, jce);
753     }
754    
755 jpye 2584 static double ascending_longitude_moon(double jce)
756 jpye 2582 {
757     return third_order_polynomial(1.0/450000.0, 0.0020708, -1934.136261, 125.04452, jce);
758     }
759    
760 jpye 2584 static double xy_term_summation(int i, double x[TERM_X_COUNT])
761 jpye 2582 {
762     int j;
763     double sum=0;
764    
765     for (j = 0; j < TERM_Y_COUNT; j++)
766     sum += x[j]*Y_TERMS[i][j];
767    
768     return sum;
769     }
770    
771 jpye 2584 static void nutation_longitude_and_obliquity(double jce, double x[TERM_X_COUNT], double *del_psi,
772 jpye 2582 double *del_epsilon)
773     {
774     int i;
775     double xy_term_sum, sum_psi=0, sum_epsilon=0;
776    
777     for (i = 0; i < Y_COUNT; i++) {
778     xy_term_sum = deg2rad(xy_term_summation(i, x));
779     sum_psi += (PE_TERMS[i][TERM_PSI_A] + jce*PE_TERMS[i][TERM_PSI_B])*sin(xy_term_sum);
780     sum_epsilon += (PE_TERMS[i][TERM_EPS_C] + jce*PE_TERMS[i][TERM_EPS_D])*cos(xy_term_sum);
781     }
782    
783     *del_psi = sum_psi / 36000000.0;
784     *del_epsilon = sum_epsilon / 36000000.0;
785     }
786    
787 jpye 2584 static double ecliptic_mean_obliquity(double jme)
788 jpye 2582 {
789     double u = jme/10.0;
790    
791     return 84381.448 + u*(-4680.93 + u*(-1.55 + u*(1999.25 + u*(-51.38 + u*(-249.67 +
792     u*( -39.05 + u*( 7.12 + u*( 27.87 + u*( 5.79 + u*2.45)))))))));
793     }
794    
795 jpye 2584 static double ecliptic_true_obliquity(double delta_epsilon, double epsilon0)
796 jpye 2582 {
797     return delta_epsilon + epsilon0/3600.0;
798     }
799    
800 jpye 2584 static double aberration_correction(double r)
801 jpye 2582 {
802     return -20.4898 / (3600.0*r);
803     }
804    
805 jpye 2584 static double apparent_sun_longitude(double theta, double delta_psi, double delta_tau)
806 jpye 2582 {
807     return theta + delta_psi + delta_tau;
808     }
809    
810 jpye 2584 static double greenwich_mean_sidereal_time (double jd, double jc)
811 jpye 2582 {
812     return limit_degrees(280.46061837 + 360.98564736629 * (jd - 2451545.0) +
813     jc*jc*(0.000387933 - jc/38710000.0));
814     }
815    
816 jpye 2584 static double greenwich_sidereal_time (double nu0, double delta_psi, double epsilon)
817 jpye 2582 {
818     return nu0 + delta_psi*cos(deg2rad(epsilon));
819     }
820    
821 jpye 2584 static double geocentric_sun_right_ascension(double lamda, double epsilon, double beta)
822 jpye 2582 {
823     double lamda_rad = deg2rad(lamda);
824     double epsilon_rad = deg2rad(epsilon);
825    
826     return limit_degrees(rad2deg(atan2(sin(lamda_rad)*cos(epsilon_rad) -
827     tan(deg2rad(beta))*sin(epsilon_rad), cos(lamda_rad))));
828     }
829    
830 jpye 2584 static double geocentric_sun_declination(double beta, double epsilon, double lamda)
831 jpye 2582 {
832     double beta_rad = deg2rad(beta);
833     double epsilon_rad = deg2rad(epsilon);
834    
835     return rad2deg(asin(sin(beta_rad)*cos(epsilon_rad) +
836     cos(beta_rad)*sin(epsilon_rad)*sin(deg2rad(lamda))));
837     }
838    
839 jpye 2584 static double observer_hour_angle(double nu, double longitude, double alpha_deg)
840 jpye 2582 {
841     return limit_degrees(nu + longitude - alpha_deg);
842     }
843    
844 jpye 2584 static double sun_equatorial_horizontal_parallax(double r)
845 jpye 2582 {
846     return 8.794 / (3600.0 * r);
847     }
848    
849 jpye 2584 static void sun_right_ascension_parallax_and_topocentric_dec(double latitude, double elevation,
850 jpye 2582 double xi, double h, double delta, double *delta_alpha, double *delta_prime)
851     {
852     double delta_alpha_rad;
853     double lat_rad = deg2rad(latitude);
854     double xi_rad = deg2rad(xi);
855     double h_rad = deg2rad(h);
856     double delta_rad = deg2rad(delta);
857     double u = atan(0.99664719 * tan(lat_rad));
858     double y = 0.99664719 * sin(u) + elevation*sin(lat_rad)/6378140.0;
859     double x = cos(u) + elevation*cos(lat_rad)/6378140.0;
860    
861     delta_alpha_rad = atan2( - x*sin(xi_rad) *sin(h_rad),
862     cos(delta_rad) - x*sin(xi_rad) *cos(h_rad));
863    
864     *delta_prime = rad2deg(atan2((sin(delta_rad) - y*sin(xi_rad))*cos(delta_alpha_rad),
865     cos(delta_rad) - x*sin(xi_rad) *cos(h_rad)));
866    
867     *delta_alpha = rad2deg(delta_alpha_rad);
868     }
869    
870 jpye 2584 static double topocentric_sun_right_ascension(double alpha_deg, double delta_alpha)
871 jpye 2582 {
872     return alpha_deg + delta_alpha;
873     }
874    
875 jpye 2584 static double topocentric_local_hour_angle(double h, double delta_alpha)
876 jpye 2582 {
877     return h - delta_alpha;
878     }
879    
880 jpye 2584 static double topocentric_elevation_angle(double latitude, double delta_prime, double h_prime)
881 jpye 2582 {
882     double lat_rad = deg2rad(latitude);
883     double delta_prime_rad = deg2rad(delta_prime);
884    
885     return rad2deg(asin(sin(lat_rad)*sin(delta_prime_rad) +
886     cos(lat_rad)*cos(delta_prime_rad) * cos(deg2rad(h_prime))));
887     }
888    
889 jpye 2584 static double atmospheric_refraction_correction(double pressure, double temperature,
890 jpye 2582 double atmos_refract, double e0)
891     {
892     double del_e = 0;
893    
894     if (e0 >= -1*(SUN_RADIUS + atmos_refract))
895     del_e = (pressure / 1010.0) * (283.0 / (273.0 + temperature)) *
896     1.02 / (60.0 * tan(deg2rad(e0 + 10.3/(e0 + 5.11))));
897    
898     return del_e;
899     }
900    
901 jpye 2584 static double topocentric_elevation_angle_corrected(double e0, double delta_e)
902 jpye 2582 {
903     return e0 + delta_e;
904     }
905    
906 jpye 2584 static double topocentric_zenith_angle(double e)
907 jpye 2582 {
908     return 90.0 - e;
909     }
910    
911 jpye 2584 static double topocentric_azimuth_angle_neg180_180(double h_prime, double latitude, double delta_prime)
912 jpye 2582 {
913     double h_prime_rad = deg2rad(h_prime);
914     double lat_rad = deg2rad(latitude);
915    
916     return rad2deg(atan2(sin(h_prime_rad),
917     cos(h_prime_rad)*sin(lat_rad) - tan(deg2rad(delta_prime))*cos(lat_rad)));
918     }
919    
920 jpye 2584 static double topocentric_azimuth_angle_zero_360(double azimuth180)
921 jpye 2582 {
922     return azimuth180 + 180.0;
923     }
924    
925 jpye 2584 static double surface_incidence_angle(double zenith, double azimuth180, double azm_rotation,
926 jpye 2582 double slope)
927     {
928     double zenith_rad = deg2rad(zenith);
929     double slope_rad = deg2rad(slope);
930    
931     return rad2deg(acos(cos(zenith_rad)*cos(slope_rad) +
932     sin(slope_rad )*sin(zenith_rad) * cos(deg2rad(azimuth180 - azm_rotation))));
933     }
934    
935 jpye 2584 static double sun_mean_longitude(double jme)
936 jpye 2582 {
937     return limit_degrees(280.4664567 + jme*(360007.6982779 + jme*(0.03032028 +
938     jme*(1/49931.0 + jme*(-1/15300.0 + jme*(-1/2000000.0))))));
939     }
940    
941 jpye 2584 static double eot(double m, double alpha, double del_psi, double epsilon)
942 jpye 2582 {
943     return limit_minutes(4.0*(m - 0.0057183 - alpha + del_psi*cos(deg2rad(epsilon))));
944     }
945    
946 jpye 2584 static double approx_sun_transit_time(double alpha_zero, double longitude, double nu)
947 jpye 2582 {
948     return (alpha_zero - longitude - nu) / 360.0;
949     }
950    
951 jpye 2584 static double sun_hour_angle_at_rise_set(double latitude, double delta_zero, double h0_prime)
952 jpye 2582 {
953     double h0 = -99999;
954     double latitude_rad = deg2rad(latitude);
955     double delta_zero_rad = deg2rad(delta_zero);
956     double argument = (sin(deg2rad(h0_prime)) - sin(latitude_rad)*sin(delta_zero_rad)) /
957     (cos(latitude_rad)*cos(delta_zero_rad));
958    
959     if (fabs(argument) <= 1) h0 = limit_degrees180(rad2deg(acos(argument)));
960    
961     return h0;
962     }
963    
964 jpye 2584 static void approx_sun_rise_and_set(double *m_rts, double h0)
965 jpye 2582 {
966     double h0_dfrac = h0/360.0;
967    
968     m_rts[SUN_RISE] = limit_zero2one(m_rts[SUN_TRANSIT] - h0_dfrac);
969     m_rts[SUN_SET] = limit_zero2one(m_rts[SUN_TRANSIT] + h0_dfrac);
970     m_rts[SUN_TRANSIT] = limit_zero2one(m_rts[SUN_TRANSIT]);
971     }
972    
973 jpye 2584 static double rts_alpha_delta_prime(double *ad, double n)
974 jpye 2582 {
975     double a = ad[JD_ZERO] - ad[JD_MINUS];
976     double b = ad[JD_PLUS] - ad[JD_ZERO];
977    
978     if (fabs(a) >= 2.0) a = limit_zero2one(a);
979     if (fabs(b) >= 2.0) b = limit_zero2one(b);
980    
981     return ad[JD_ZERO] + n * (a + b + (b-a)*n)/2.0;
982     }
983    
984 jpye 2584 static double rts_sun_altitude(double latitude, double delta_prime, double h_prime)
985 jpye 2582 {
986     double latitude_rad = deg2rad(latitude);
987     double delta_prime_rad = deg2rad(delta_prime);
988    
989     return rad2deg(asin(sin(latitude_rad)*sin(delta_prime_rad) +
990     cos(latitude_rad)*cos(delta_prime_rad)*cos(deg2rad(h_prime))));
991     }
992    
993 jpye 2584 static double sun_rise_and_set(double *m_rts, double *h_rts, double *delta_prime, double latitude,
994 jpye 2582 double *h_prime, double h0_prime, int sun)
995     {
996     return m_rts[sun] + (h_rts[sun] - h0_prime) /
997     (360.0*cos(deg2rad(delta_prime[sun]))*cos(deg2rad(latitude))*sin(deg2rad(h_prime[sun])));
998     }
999    
1000     ////////////////////////////////////////////////////////////////////////////////////////////////
1001     // Calculate required SPA parameters to get the right ascension (alpha) and declination (delta)
1002     // Note: JD must be already calculated and in structure
1003     ////////////////////////////////////////////////////////////////////////////////////////////////
1004 jpye 2584 static void calculate_geocentric_sun_right_ascension_and_declination(spa_data *spa)
1005 jpye 2582 {
1006     double x[TERM_X_COUNT];
1007    
1008     spa->jc = julian_century(spa->jd);
1009    
1010     spa->jde = julian_ephemeris_day(spa->jd, spa->delta_t);
1011     spa->jce = julian_ephemeris_century(spa->jde);
1012     spa->jme = julian_ephemeris_millennium(spa->jce);
1013    
1014     spa->l = earth_heliocentric_longitude(spa->jme);
1015     spa->b = earth_heliocentric_latitude(spa->jme);
1016     spa->r = earth_radius_vector(spa->jme);
1017    
1018     spa->theta = geocentric_longitude(spa->l);
1019     spa->beta = geocentric_latitude(spa->b);
1020    
1021     x[TERM_X0] = spa->x0 = mean_elongation_moon_sun(spa->jce);
1022     x[TERM_X1] = spa->x1 = mean_anomaly_sun(spa->jce);
1023     x[TERM_X2] = spa->x2 = mean_anomaly_moon(spa->jce);
1024     x[TERM_X3] = spa->x3 = argument_latitude_moon(spa->jce);
1025     x[TERM_X4] = spa->x4 = ascending_longitude_moon(spa->jce);
1026    
1027     nutation_longitude_and_obliquity(spa->jce, x, &(spa->del_psi), &(spa->del_epsilon));
1028    
1029     spa->epsilon0 = ecliptic_mean_obliquity(spa->jme);
1030     spa->epsilon = ecliptic_true_obliquity(spa->del_epsilon, spa->epsilon0);
1031    
1032     spa->del_tau = aberration_correction(spa->r);
1033     spa->lamda = apparent_sun_longitude(spa->theta, spa->del_psi, spa->del_tau);
1034     spa->nu0 = greenwich_mean_sidereal_time (spa->jd, spa->jc);
1035     spa->nu = greenwich_sidereal_time (spa->nu0, spa->del_psi, spa->epsilon);
1036    
1037     spa->alpha = geocentric_sun_right_ascension(spa->lamda, spa->epsilon, spa->beta);
1038     spa->delta = geocentric_sun_declination(spa->beta, spa->epsilon, spa->lamda);
1039     }
1040    
1041     ////////////////////////////////////////////////////////////////////////
1042     // Calculate Equation of Time (EOT) and Sun Rise, Transit, & Set (RTS)
1043     ////////////////////////////////////////////////////////////////////////
1044    
1045 jpye 2584 static void calculate_eot_and_sun_rise_transit_set(spa_data *spa)
1046 jpye 2582 {
1047     spa_data sun_rts;
1048     double nu, m, h0, n;
1049     double alpha[JD_COUNT], delta[JD_COUNT];
1050     double m_rts[SUN_COUNT], nu_rts[SUN_COUNT], h_rts[SUN_COUNT];
1051     double alpha_prime[SUN_COUNT], delta_prime[SUN_COUNT], h_prime[SUN_COUNT];
1052     double h0_prime = -1*(SUN_RADIUS + spa->atmos_refract);
1053     int i;
1054    
1055     sun_rts = *spa;
1056     m = sun_mean_longitude(spa->jme);
1057     spa->eot = eot(m, spa->alpha, spa->del_psi, spa->epsilon);
1058    
1059     sun_rts.hour = sun_rts.minute = sun_rts.second = 0;
1060     sun_rts.timezone = 0.0;
1061    
1062     sun_rts.jd = julian_day (sun_rts.year, sun_rts.month, sun_rts.day,
1063     sun_rts.hour, sun_rts.minute, sun_rts.second, sun_rts.timezone);
1064    
1065     calculate_geocentric_sun_right_ascension_and_declination(&sun_rts);
1066     nu = sun_rts.nu;
1067    
1068     sun_rts.delta_t = 0;
1069     sun_rts.jd--;
1070     for (i = 0; i < JD_COUNT; i++) {
1071     calculate_geocentric_sun_right_ascension_and_declination(&sun_rts);
1072     alpha[i] = sun_rts.alpha;
1073     delta[i] = sun_rts.delta;
1074     sun_rts.jd++;
1075     }
1076    
1077     m_rts[SUN_TRANSIT] = approx_sun_transit_time(alpha[JD_ZERO], spa->longitude, nu);
1078     h0 = sun_hour_angle_at_rise_set(spa->latitude, delta[JD_ZERO], h0_prime);
1079    
1080     if (h0 >= 0) {
1081    
1082     approx_sun_rise_and_set(m_rts, h0);
1083    
1084     for (i = 0; i < SUN_COUNT; i++) {
1085    
1086     nu_rts[i] = nu + 360.985647*m_rts[i];
1087    
1088     n = m_rts[i] + spa->delta_t/86400.0;
1089     alpha_prime[i] = rts_alpha_delta_prime(alpha, n);
1090     delta_prime[i] = rts_alpha_delta_prime(delta, n);
1091    
1092     h_prime[i] = limit_degrees180pm(nu_rts[i] + spa->longitude - alpha_prime[i]);
1093    
1094     h_rts[i] = rts_sun_altitude(spa->latitude, delta_prime[i], h_prime[i]);
1095     }
1096    
1097     spa->srha = h_prime[SUN_RISE];
1098     spa->ssha = h_prime[SUN_SET];
1099     spa->sta = h_rts[SUN_TRANSIT];
1100    
1101     spa->suntransit = dayfrac_to_local_hr(m_rts[SUN_TRANSIT] - h_prime[SUN_TRANSIT] / 360.0,
1102     spa->timezone);
1103    
1104     spa->sunrise = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime,
1105     spa->latitude, h_prime, h0_prime, SUN_RISE), spa->timezone);
1106    
1107     spa->sunset = dayfrac_to_local_hr(sun_rise_and_set(m_rts, h_rts, delta_prime,
1108     spa->latitude, h_prime, h0_prime, SUN_SET), spa->timezone);
1109    
1110     } else spa->srha= spa->ssha= spa->sta= spa->suntransit= spa->sunrise= spa->sunset= -99999;
1111    
1112     }
1113    
1114     ///////////////////////////////////////////////////////////////////////////////////////////
1115     // Calculate all SPA parameters and put into structure
1116     // Note: All inputs values (listed in header file) must already be in structure
1117     ///////////////////////////////////////////////////////////////////////////////////////////
1118 jpye 2584 #include <stdio.h>
1119    
1120 jpye 2582 int spa_calculate(spa_data *spa)
1121     {
1122     int result;
1123    
1124     result = validate_inputs(spa);
1125    
1126 jpye 2584 if (result == 0){
1127    
1128     //fprintf(stderr,"\n\n\nfunction = %d = %s\n\n\n",spa->function,(spa->function == SPA_ZA_JD?"SPA_ZA_JD":"other"));
1129    
1130     /* only calculate the JD from the datetime field if requested to do so, otherwise
1131     we should be able to use the JD directly as a time input. */
1132     switch(spa->function){
1133     case SPA_ZA:
1134     case SPA_ZA_INC:
1135     case SPA_ZA_RTS:
1136     case SPA_ALL:
1137     case SPA_JD:
1138     spa->jd = julian_day (spa->year, spa->month, spa->day,
1139 jpye 2582 spa->hour, spa->minute, spa->second, spa->timezone);
1140 jpye 2584 if(spa->function == SPA_JD){
1141     return result;
1142     }
1143     break;
1144     default:
1145     break;
1146     }
1147 jpye 2582
1148     calculate_geocentric_sun_right_ascension_and_declination(spa);
1149    
1150     spa->h = observer_hour_angle(spa->nu, spa->longitude, spa->alpha);
1151     spa->xi = sun_equatorial_horizontal_parallax(spa->r);
1152    
1153     sun_right_ascension_parallax_and_topocentric_dec(spa->latitude, spa->elevation, spa->xi,
1154     spa->h, spa->delta, &(spa->del_alpha), &(spa->delta_prime));
1155    
1156     spa->alpha_prime = topocentric_sun_right_ascension(spa->alpha, spa->del_alpha);
1157     spa->h_prime = topocentric_local_hour_angle(spa->h, spa->del_alpha);
1158    
1159     spa->e0 = topocentric_elevation_angle(spa->latitude, spa->delta_prime, spa->h_prime);
1160     spa->del_e = atmospheric_refraction_correction(spa->pressure, spa->temperature,
1161     spa->atmos_refract, spa->e0);
1162     spa->e = topocentric_elevation_angle_corrected(spa->e0, spa->del_e);
1163    
1164     spa->zenith = topocentric_zenith_angle(spa->e);
1165     spa->azimuth180 = topocentric_azimuth_angle_neg180_180(spa->h_prime, spa->latitude,
1166     spa->delta_prime);
1167     spa->azimuth = topocentric_azimuth_angle_zero_360(spa->azimuth180);
1168    
1169 jpye 2584 if((spa->function == SPA_ZA_INC) || (spa->function == SPA_ZA_INC_JD) || (spa->function == SPA_ALL)){
1170 jpye 2582 spa->incidence = surface_incidence_angle(spa->zenith, spa->azimuth180,
1171     spa->azm_rotation, spa->slope);
1172 jpye 2584 }
1173 jpye 2582
1174 jpye 2584 if((spa->function == SPA_ZA_RTS) || (spa->function == SPA_ZA_RTS_JD) || (spa->function == SPA_ALL)){
1175 jpye 2582 calculate_eot_and_sun_rise_transit_set(spa);
1176 jpye 2584 }
1177 jpye 2582 }
1178    
1179     return result;
1180     }
1181     ///////////////////////////////////////////////////////////////////////////////////////////

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