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

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