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