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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2599 - (show 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 /*
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
6 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 /////////////////////////////////////////////
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 // NOTICE
41 //
42 //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 //
47 //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 //
51 //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 //
58 ///////////////////////////////////////////////////////////////////////////////////////////////
59
60
61 ///////////////////////////////////////////////////////////////////////////////////////////////
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 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
129 ///////////////////////////////////////////////////
130 /// Earth Periodic Terms
131 ///////////////////////////////////////////////////
132 static 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 static 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 static 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 static 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 static 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 static double rad2deg(double radians)
505 {
506 return (180.0/PI)*radians;
507 }
508
509 static double deg2rad(double degrees)
510 {
511 return (PI/180.0)*degrees;
512 }
513
514 static 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 static 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 static 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 static 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 static 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 static double dayfrac_to_local_hr(double dayfrac, double timezone)
569 {
570 return 24.0*limit_zero2one(dayfrac + timezone/24.0);
571 }
572
573 static 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 static int validate_inputs(spa_data *spa)
580 {
581 switch(spa->function){
582 case SPA_ZA_JD:
583 case SPA_ZA_INC_JD:
584 case SPA_ZA_RTS_JD:
585 if ((spa->jd < 990557.5) || (spa->jd > 3912879.500000)) return 16;
586 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 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 double julian_day (int year, int month, int day, int hour, int minute, int second, double tz)
618 {
619 double day_decimal, julian_day, a;
620
621 //fprintf(stderr,"Julian day for %d/%d/%d %d:%d:%d (%+0.2f) is ",day,month,year,hour,minute,second,tz);
622
623 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 //fprintf(stderr,"%f\n",julian_day);
638 return julian_day;
639 }
640
641 static double julian_century(double jd)
642 {
643 return (jd-2451545.0)/36525.0;
644 }
645
646 static double julian_ephemeris_day(double jd, double delta_t)
647 {
648 return jd+delta_t/86400.0;
649 }
650
651 static double julian_ephemeris_century(double jde)
652 {
653 return (jde - 2451545.0)/36525.0;
654 }
655
656 static double julian_ephemeris_millennium(double jce)
657 {
658 return (jce/10.0);
659 }
660
661 static double earth_periodic_term_summation(const double terms[][TERM_COUNT], int count, double jme)
662 {
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 static double earth_values(double term_sum[], int count, double jme)
673 {
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 static double earth_heliocentric_longitude(double jme)
686 {
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 static double earth_heliocentric_latitude(double jme)
698 {
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 static double earth_radius_vector(double jme)
710 {
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 static double geocentric_longitude(double l)
722 {
723 double theta = l + 180.0;
724
725 if (theta >= 360.0) theta -= 360.0;
726
727 return theta;
728 }
729
730 static double geocentric_latitude(double b)
731 {
732 return -b;
733 }
734
735 static double mean_elongation_moon_sun(double jce)
736 {
737 return third_order_polynomial(1.0/189474.0, -0.0019142, 445267.11148, 297.85036, jce);
738 }
739
740 static double mean_anomaly_sun(double jce)
741 {
742 return third_order_polynomial(-1.0/300000.0, -0.0001603, 35999.05034, 357.52772, jce);
743 }
744
745 static double mean_anomaly_moon(double jce)
746 {
747 return third_order_polynomial(1.0/56250.0, 0.0086972, 477198.867398, 134.96298, jce);
748 }
749
750 static double argument_latitude_moon(double jce)
751 {
752 return third_order_polynomial(1.0/327270.0, -0.0036825, 483202.017538, 93.27191, jce);
753 }
754
755 static double ascending_longitude_moon(double jce)
756 {
757 return third_order_polynomial(1.0/450000.0, 0.0020708, -1934.136261, 125.04452, jce);
758 }
759
760 static double xy_term_summation(int i, double x[TERM_X_COUNT])
761 {
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 static void nutation_longitude_and_obliquity(double jce, double x[TERM_X_COUNT], double *del_psi,
772 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 static double ecliptic_mean_obliquity(double jme)
788 {
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 static double ecliptic_true_obliquity(double delta_epsilon, double epsilon0)
796 {
797 return delta_epsilon + epsilon0/3600.0;
798 }
799
800 static double aberration_correction(double r)
801 {
802 return -20.4898 / (3600.0*r);
803 }
804
805 static double apparent_sun_longitude(double theta, double delta_psi, double delta_tau)
806 {
807 return theta + delta_psi + delta_tau;
808 }
809
810 static double greenwich_mean_sidereal_time (double jd, double jc)
811 {
812 return limit_degrees(280.46061837 + 360.98564736629 * (jd - 2451545.0) +
813 jc*jc*(0.000387933 - jc/38710000.0));
814 }
815
816 static double greenwich_sidereal_time (double nu0, double delta_psi, double epsilon)
817 {
818 return nu0 + delta_psi*cos(deg2rad(epsilon));
819 }
820
821 static double geocentric_sun_right_ascension(double lamda, double epsilon, double beta)
822 {
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 static double geocentric_sun_declination(double beta, double epsilon, double lamda)
831 {
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 static double observer_hour_angle(double nu, double longitude, double alpha_deg)
840 {
841 return limit_degrees(nu + longitude - alpha_deg);
842 }
843
844 static double sun_equatorial_horizontal_parallax(double r)
845 {
846 return 8.794 / (3600.0 * r);
847 }
848
849 static void sun_right_ascension_parallax_and_topocentric_dec(double latitude, double elevation,
850 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 static double topocentric_sun_right_ascension(double alpha_deg, double delta_alpha)
871 {
872 return alpha_deg + delta_alpha;
873 }
874
875 static double topocentric_local_hour_angle(double h, double delta_alpha)
876 {
877 return h - delta_alpha;
878 }
879
880 static double topocentric_elevation_angle(double latitude, double delta_prime, double h_prime)
881 {
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 static double atmospheric_refraction_correction(double pressure, double temperature,
890 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 static double topocentric_elevation_angle_corrected(double e0, double delta_e)
902 {
903 return e0 + delta_e;
904 }
905
906 static double topocentric_zenith_angle(double e)
907 {
908 return 90.0 - e;
909 }
910
911 static double topocentric_azimuth_angle_neg180_180(double h_prime, double latitude, double delta_prime)
912 {
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 static double topocentric_azimuth_angle_zero_360(double azimuth180)
921 {
922 return azimuth180 + 180.0;
923 }
924
925 static double surface_incidence_angle(double zenith, double azimuth180, double azm_rotation,
926 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 static double sun_mean_longitude(double jme)
936 {
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 static double eot(double m, double alpha, double del_psi, double epsilon)
942 {
943 return limit_minutes(4.0*(m - 0.0057183 - alpha + del_psi*cos(deg2rad(epsilon))));
944 }
945
946 static double approx_sun_transit_time(double alpha_zero, double longitude, double nu)
947 {
948 return (alpha_zero - longitude - nu) / 360.0;
949 }
950
951 static double sun_hour_angle_at_rise_set(double latitude, double delta_zero, double h0_prime)
952 {
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 static void approx_sun_rise_and_set(double *m_rts, double h0)
965 {
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 static double rts_alpha_delta_prime(double *ad, double n)
974 {
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 static double rts_sun_altitude(double latitude, double delta_prime, double h_prime)
985 {
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 static double sun_rise_and_set(double *m_rts, double *h_rts, double *delta_prime, double latitude,
994 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 static void calculate_geocentric_sun_right_ascension_and_declination(spa_data *spa)
1005 {
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 static void calculate_eot_and_sun_rise_transit_set(spa_data *spa)
1046 {
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 #include <stdio.h>
1119
1120 int spa_calculate(spa_data *spa)
1121 {
1122 int result;
1123
1124 result = validate_inputs(spa);
1125
1126 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 spa->hour, spa->minute, spa->second, spa->timezone);
1140 if(spa->function == SPA_JD){
1141 return result;
1142 }
1143 break;
1144 default:
1145 break;
1146 }
1147
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 if((spa->function == SPA_ZA_INC) || (spa->function == SPA_ZA_INC_JD) || (spa->function == SPA_ALL)){
1170 spa->incidence = surface_incidence_angle(spa->zenith, spa->azimuth180,
1171 spa->azm_rotation, spa->slope);
1172 }
1173
1174 if((spa->function == SPA_ZA_RTS) || (spa->function == SPA_ZA_RTS_JD) || (spa->function == SPA_ALL)){
1175 calculate_eot_and_sun_rise_transit_set(spa);
1176 }
1177 }
1178
1179 return result;
1180 }
1181 ///////////////////////////////////////////////////////////////////////////////////////////

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