1 /************************************************************************
2 * Copyright (c) 2003-2006, 2009-2010 Arabeyes, Thamer Mahmoud
3 *
4 * A full featured Muslim Prayer Times calculator
5 *
6 * Most of the astronomical values and formulas used in this file are based
7 * upon a subset of the VSOP87 planetary theory developed by Jean Meeus. Page
8 * and formula numbers in-line below are references to his book: Astronomical
9 * Algorithms. Willmann-Bell, second edition, 1998.
10 *
11 * (www.arabeyes.org - under LGPL license - see COPYING file)
12 ************************************************************************/
13
14 #include "astro.h"
15
16 typedef struct
17 {
18 double dec;
19 double ra;
20 double sidtime;
21 double dra;
22 double rsum;
23
24 } AstroDay ;
25
26 enum Type { SUNRISE,
27 SUNSET };
28
29 static void computeAstroDay(double JD, AstroDay* astroday);
30 static void computeTopAstro(const Location* loc, const Astro* astro, Astro* tastro);
31 static double getRiseSet (const Location* loc, const Astro* astro, int type);
32 static double getRefraction(const Location* loc, double sunAlt);
33 static double limitAngle180(double L);
34 static double limitAngle(double L);
35 static double limitAngle1(double L);
36 static double limitAngle180between(double L);
37
getSunrise(const Location * loc,const Astro * tastro)38 double getSunrise (const Location* loc, const Astro* tastro)
39 {
40 return getRiseSet (loc, tastro, SUNRISE);
41 }
42
getSunset(const Location * loc,const Astro * tastro)43 double getSunset (const Location* loc, const Astro* tastro)
44 {
45 return getRiseSet (loc, tastro, SUNSET);
46 }
47
getTransit(double lon,const Astro * tastro)48 double getTransit(double lon, const Astro* tastro)
49 {
50
51 double M, sidG;
52 double ra0=tastro->ra[0], ra2=tastro->ra[2];
53 double A, H;
54
55 M = ((tastro->ra[1] - lon - tastro->sid[1]) / 360.0);
56 M = limitAngle1(M);
57
58 /* Sidereal time at Greenwich (p. 103) */
59 sidG = tastro->sid[1] + 360.985647 * M;
60
61 if (tastro->ra[1] > 350 && tastro->ra[2] < 10)
62 ra2 += 360;
63 if (tastro->ra[0] > 350 && tastro->ra[1] < 10)
64 ra0 = 0;
65
66 /* Interpolation of 1-day intervals (pp. 24-25) */
67 A = tastro->ra[1] +
68 (M * ((tastro->ra[1] - ra0) + ( ra2 - tastro->ra[1]) +
69 (( ra2 - tastro->ra[1]) - (tastro->ra[1] - ra0)) * M) / 2.0 );
70
71 H = limitAngle180between(sidG + lon - A);
72
73 return 24.0 * (M - (H/360.0));
74 }
75
getRiseSet(const Location * loc,const Astro * tastro,int type)76 static double getRiseSet (const Location* loc, const Astro* tastro, int type)
77 {
78 /* p. 101 */
79 double lhour, M, sidG, ra0, ra2;
80 double A, B, H, sunAlt, delM, tH, rDec, rLat, rB;
81 double part1, part2, part3;
82
83 rDec = DEG_TO_RAD(tastro->dec[1]);
84 rLat = DEG_TO_RAD(loc->degreeLat);
85
86 ra0=tastro->ra[0];
87 ra2=tastro->ra[2];
88
89 /* Compute the hour angle */
90 part1 = sin(DEG_TO_RAD(CENTER_OF_SUN_ANGLE)) - (sin (rLat) * sin (rDec));
91 part2 = cos (rLat) * cos (rDec);
92 part3 = part1 / part2;
93
94 if ( part3 < -INVALID_TRIGGER || part3 > INVALID_TRIGGER)
95 return 99;
96
97 lhour = limitAngle180 (( RAD_TO_DEG (acos (part3))));
98
99 /* Eastern Longitudes are positive throughout this file. */
100 M = ((tastro->ra[1] - loc->degreeLong - tastro->sid[1]) / 360.0);
101
102 if (type == SUNRISE)
103 M = M - (lhour/360.0);
104 if (type == SUNSET)
105 M = M + (lhour/360.0);
106
107 M = limitAngle1(M);
108
109 /* Sidereal time at Greenwich (p. 103) */
110 sidG = limitAngle(tastro->sid[1] + 360.985647 * M);
111
112 ra0 = tastro->ra[0];
113 ra2 = tastro->ra[2];
114
115 if (tastro->ra[1] > 350 && tastro->ra[2] < 10)
116 ra2 += 360;
117 if (tastro->ra[0] > 350 && tastro->ra[1] < 10)
118 ra0 = 0;
119
120 /* Interpolation of 1-day intervals (pp. 24-25) */
121 A = tastro->ra[1] + (M * (( tastro->ra[1] - ra0) +
122 (ra2 - tastro->ra[1] ) +
123 (( ra2 - tastro->ra[1] ) -
124 ( tastro->ra[1] - ra0)) * M) / 2.0 );
125
126 B = tastro->dec[1] + (M * ((tastro->dec[1] - tastro->dec[0]) +
127 (tastro->dec[2] - tastro->dec[1]) +
128 ((tastro->dec[2] - tastro->dec[1]) -
129 (tastro->dec[1] - tastro->dec[0])) * M) / 2.0 );
130 rB = DEG_TO_RAD(B);
131
132 H = limitAngle180between(sidG + loc->degreeLong - A);
133
134 tH = DEG_TO_RAD(H) - tastro->dra[1];
135
136 /* Airless Sun's altitude at local horizontal coordinates (p. 93, 13.6) */
137 sunAlt = RAD_TO_DEG(asin ( sin(rLat) * sin(rB)
138 + cos(rLat) * cos(rB)
139 * cos(tH) ));
140
141 sunAlt += getRefraction(loc, sunAlt);
142
143 /* (p. 103) */
144 delM = (sunAlt - CENTER_OF_SUN_ANGLE) / (360.0 * cos(rB) * cos(rLat)
145 * sin(tH));
146
147 return (M + delM) * 24.0;
148
149 }
150
151
getRefraction(const Location * loc,double sunAlt)152 double getRefraction(const Location* loc, double sunAlt)
153 {
154 double part1, part2;
155
156 part1 = (loc->pressure/1010.0) * (283/(273 + loc->temperature));
157 part2 = 1.02 / (RAD_TO_DEG(tan(DEG_TO_RAD(sunAlt + (10.3/(sunAlt + 5.11))))) + 0.0019279);
158
159 return (part1 * part2) / 60.0;
160 }
161
162
163 const double DT2[]={
164 63.4673, 63.8285, 64.0908, 64.2998, 64.4734, /* 1999-2003 */
165 64.5736, 64.7052, 64.8452, 65.1464, 65.4574, /* 2004-2008 */
166 65.7768, /* 2009 */
167 66.5, 67.1, 68, 68, 69, /* 2010-2014 predictions */
168 69, 70, 70 /* 2015-2017 predictions */
169 };
170
computeDeltaT(double year)171 static double computeDeltaT(double year)
172 {
173 int i;
174 double tempdt;
175 /* pp. 78-80 */
176 double t = (year - 2000) / 100.0;
177 if (year < 948)
178 return 2177 + (497 * t) + (44.1 * pow (t, 2));
179 else if (year >= 1620 && year <= 1998)
180 return 0; /* FIXIT: Support historical delta-t values for years before
181 * 1998. In the DT table above, each line represents 5 even
182 * years in the range 1620-1998. We should first complete the
183 * table to include data for both odd and even years. */
184 else if ((year > 1998 && year < 2100) || year < 1620)
185 {
186 /* NOTE: The "2017" found below this comment should be changed to
187 reflect the last year added to the DT2 table. */
188 if (year >= 1999 && year <= 2017) {
189 i = year-1999;
190 return DT2[i];
191 }
192 /* As of 2007, the two formulas given by Meeus seem to be off by
193 many seconds when compared to observed values found at
194 <http://maia.usno.navy.mil>. The DT2 table overrides these values
195 with observed and predicted ones for delta-t (on January, other
196 months not yet supported). Extrapolated (and wrong) values are still
197 used for years after 2017. */
198 else tempdt = 102 + (102 * t) + (25.3 * pow (t, 2));
199
200 if (year >= 2000)
201 return tempdt + (0.37 * (year - 2100));
202 else return tempdt;
203 }
204 return 0;
205 }
206
207
getJulianDay(const Date * date,double gmt)208 double getJulianDay(const Date* date, double gmt)
209 {
210 double jdB=0, jdY, jdM, JD;
211
212 jdY=date->year;
213 jdM=date->month;
214
215 if (date->month <= 2) {
216 jdY--;
217 jdM+=12;
218 }
219
220 if (date->year < 1)
221 jdY++;
222
223 if ((date->year > 1582) || ((date->year == 1582) &&
224 ((date->month > 10) ||
225 ((date->month == 10) && (date->day >= 4)))))
226 jdB = 2 - floor(jdY/100.0) + floor((jdY/100.0)/4.0);
227
228 JD = floor(365.25 * (jdY + 4716.0)) + floor(30.6001 * (jdM + 1))
229 + (date->day + (-gmt)/24.0) + jdB - 1524.5 ;
230
231 JD = JD + (computeDeltaT(date->year) / 86400.0);
232 return JD;
233
234 }
235
236 const double L0[64][3]={
237 {175347046, 0, 0},
238 {3341656, 4.6692568, 6283.07585},
239 {34894, 4.6261, 12566.1517},
240 {3497, 2.7441, 5753.3849},
241 {3418, 2.8289, 3.5231},
242 {3136, 3.6277, 77713.7715},
243 {2676, 4.4181, 7860.4194},
244 {2343, 6.1352, 3930.2097},
245 {1324, 0.7425, 11506.7698},
246 {1273, 2.0371, 529.691},
247 {1199, 1.1096, 1577.3435},
248 {990, 5.233, 5884.927},
249 {902, 2.045, 26.298},
250 {857, 3.508, 398.149},
251 {780, 1.179, 5223.694},
252 {753, 2.533, 5507.553},
253 {505, 4.583, 18849.228},
254 {492, 4.205, 775.523},
255 {357, 2.92, 0.067},
256 {317, 5.849, 11790.629},
257 {284, 1.899, 796.298},
258 {271, 0.315, 10977.079},
259 {243, 0.345, 5486.778},
260 {206, 4.806, 2544.314},
261 {205, 1.869, 5573.143},
262 {202, 2.458, 6069.777},
263 {156, 0.833, 213.299},
264 {132, 3.411, 2942.463},
265 {126, 1.083, 20.775},
266 {115, 0.645, 0.98},
267 {103, 0.636, 4694.003},
268 {102, 0.976, 15720.839},
269 {102, 4.267, 7.114},
270 {99, 6.21, 2146.17},
271 {98, 0.68, 155.42},
272 {86, 5.98, 161000.69},
273 {85, 1.3, 6275.96},
274 {85, 3.67, 71430.7},
275 {80, 1.81, 17260.15},
276 {79, 3.04, 12036.46},
277 {75, 1.76, 5088.63},
278 {74, 3.5, 3154.69},
279 {74, 4.68, 801.82},
280 {70, 0.83, 9437.76},
281 {62, 3.98, 8827.39},
282 {61, 1.82, 7084.9},
283 {57, 2.78, 6286.6},
284 {56, 4.39, 14143.5},
285 {56, 3.47, 6279.55},
286 {52, 0.19, 12139.55},
287 {52, 1.33, 1748.02},
288 {51, 0.28, 5856.48},
289 {49, 0.49, 1194.45},
290 {41, 5.37, 8429.24},
291 {41, 2.4, 19651.05},
292 {39, 6.17, 10447.39},
293 {37, 6.04, 10213.29},
294 {37, 2.57, 1059.38},
295 {36, 1.71, 2352.87},
296 {36, 1.78, 6812.77},
297 {33, 0.59, 17789.85},
298 {30, 0.44, 83996.85},
299 {30, 2.74, 1349.87},
300 {25, 3.16, 4690.48}
301 };
302
303 const double L1[][3]={
304 {628331966747.0, 0, 0},
305 {206059, 2.678235, 6283.07585},
306 {4303, 2.6351, 12566.1517},
307 {425, 1.59, 3.523},
308 {119, 5.796, 26.298},
309 {109, 2.966, 1577.344},
310 {93, 2.59, 18849.23},
311 {72, 1.14, 529.69},
312 {68, 1.87, 398.15},
313 {67, 4.41, 5507.55},
314 {59, 2.89, 5223.69},
315 {56, 2.17, 155.42},
316 {45, 0.4, 796.3},
317 {36, 0.47, 775.52},
318 {29, 2.65, 7.11},
319 {21, 5.34, 0.98},
320 {19, 1.85, 5486.78},
321 {19, 4.97, 213.3},
322 {17, 2.99, 6275.96},
323 {16, 0.03, 2544.31},
324 {16, 1.43, 2146.17},
325 {15, 1.21, 10977.08},
326 {12, 2.83, 1748.02},
327 {12, 3.26, 5088.63},
328 {12, 5.27, 1194.45},
329 {12, 2.08, 4694},
330 {11, 0.77, 553.57},
331 {10, 1.3, 6286.6},
332 {10, 4.24, 1349.87},
333 {9, 2.7, 242.73},
334 {9, 5.64, 951.72},
335 {8, 5.3, 2352.87},
336 {6, 2.65, 9437.76},
337 {6, 4.67, 4690.48}
338 };
339
340 const double L2[][3]={
341 {52919, 0, 0},
342 {8720, 1.0721, 6283.0758},
343 {309, 0.867, 12566.152},
344 {27, 0.05, 3.52},
345 {16, 5.19, 26.3},
346 {16, 3.68, 155.42},
347 {10, 0.76, 18849.23},
348 {9, 2.06, 77713.77},
349 {7, 0.83, 775.52},
350 {5, 4.66, 1577.34},
351 {4, 1.03, 7.11},
352 {4, 3.44, 5573.14},
353 {3, 5.14, 796.3},
354 {3, 6.05, 5507.55},
355 {3, 1.19, 242.73},
356 {3, 6.12, 529.69},
357 {3, 0.31, 398.15},
358 {3, 2.28, 553.57},
359 {2, 4.38, 5223.69},
360 {2, 3.75, 0.98}
361 };
362
363 const double L3[][3]={
364 {289, 5.844, 6283.076},
365 {35, 0, 0},
366 {17, 5.49, 12566.15},
367 {3, 5.2, 155.42},
368 {1, 4.72, 3.52},
369 {1, 5.3, 18849.23},
370 {1, 5.97, 242.73}
371 };
372
373 const double L4[][3]={
374 {114, 3.142, 0},
375 {8, 4.13, 6283.08},
376 {1, 3.84, 12566.15}
377 };
378
379 const double L5[][3]={
380 {1, 3.14, 0}
381 };
382
383 const double B0[][3]={
384
385 {280, 3.199, 84334.662},
386 {102, 5.422, 5507.553},
387 {80, 3.88, 5223.69},
388 {44, 3.7, 2352.87},
389 {32, 4, 1577.34}
390 };
391
392 const double B1[][3]={
393
394 {9, 3.9, 5507.55},
395 {6, 1.73, 5223.69}
396 };
397
398 const double R0[][3]={
399 {100013989, 0, 0},
400 {1670700, 3.0984635, 6283.07585},
401 {13956, 3.05525, 12566.1517},
402 {3084, 5.1985, 77713.7715},
403 {1628, 1.1739, 5753.3849},
404 {1576, 2.8469, 7860.4194},
405 {925, 5.453, 11506.77},
406 {542, 4.564, 3930.21},
407 {472, 3.661, 5884.927},
408 {346, 0.964, 5507.553},
409 {329, 5.9, 5223.694},
410 {307, 0.299, 5573.143},
411 {243, 4.273, 11790.629},
412 {212, 5.847, 1577.344},
413 {186, 5.022, 10977.079},
414 {175, 3.012, 18849.228},
415 {110, 5.055, 5486.778},
416 {98, 0.89, 6069.78},
417 {86, 5.69, 15720.84},
418 {86, 1.27, 161000.69},
419 {65, 0.27, 17260.15},
420 {63, 0.92, 529.69},
421 {57, 2.01, 83996.85},
422 {56, 5.24, 71430.7},
423 {49, 3.25, 2544.31},
424 {47, 2.58, 775.52},
425 {45, 5.54, 9437.76},
426 {43, 6.01, 6275.96},
427 {39, 5.36, 4694},
428 {38, 2.39, 8827.39},
429 {37, 0.83, 19651.05},
430 {37, 4.9, 12139.55},
431 {36, 1.67, 12036.46},
432 {35, 1.84, 2942.46},
433 {33, 0.24, 7084.9},
434 {32, 0.18, 5088.63},
435 {32, 1.78, 398.15},
436 {28, 1.21, 6286.6},
437 {28, 1.9, 6279.55},
438 {26, 4.59, 10447.39}
439 };
440
441 const double R1[][3]={
442
443 {103019, 1.10749, 6283.07585},
444 {1721, 1.0644, 12566.1517},
445 {702, 3.142, 0},
446 {32, 1.02, 18849.23},
447 {31, 2.84, 5507.55},
448 {25, 1.32, 5223.69},
449 {18, 1.42, 1577.34},
450 {10, 5.91, 10977.08},
451 {9, 1.42, 6275.96},
452 {9, 0.27, 5486.78}
453 };
454
455 const double R2[][3]={
456
457 {4359, 5.7846, 6283.0758},
458 {124, 5.579, 12566.152},
459 {12, 3.14, 0},
460 {9, 3.63, 77713.77},
461 {6, 1.87, 5573.14},
462 {3, 5.47, 18849.23}
463
464 };
465
466 const double R3[][3]={
467 {145, 4.273, 6283.076},
468 {7, 3.92, 12566.15}
469 };
470
471 const double R4[][3]={
472 {4, 2.56, 6283.08}
473 };
474
475 const double PN[][4]={
476 {-171996, -174.2, 92025, 8.9},
477 {-13187, -1.6, 5736, -3.1},
478 {-2274, -0.2, 977, -0.5},
479 {2062, 0.2, -895, 0.5},
480 {1426, -3.4, 54, -0.1},
481 {712, 0.1, -7, 0},
482 {-517, 1.2, 224, -0.6},
483 {-386, -0.4, 200, 0},
484 {-301, 0, 129, -0.1},
485 {217, -0.5, -95, 0.3},
486 {-158, 0, 0, 0},
487 {129, 0.1, -70, 0},
488 {123, 0, -53, 0},
489 {63, 0, 0, 0},
490 {63, 0.1, -33, 0},
491 {-59, 0, 26, 0},
492 {-58, -0.1, 32, 0},
493 {-51, 0, 27, 0},
494 {48, 0, 0, 0},
495 {46, 0, -24, 0},
496 {-38, 0, 16, 0},
497 {-31, 0, 13, 0},
498 {29, 0, 0, 0},
499 {29, 0, -12, 0},
500 {26, 0, 0, 0},
501 {-22, 0, 0, 0},
502 {21, 0, -10, 0},
503 {17, -0.1, 0, 0},
504 {16, 0, -8, 0},
505 {-16, 0.1, 7, 0},
506 {-15, 0, 9, 0},
507 {-13, 0, 7, 0},
508 {-12, 0, 6, 0},
509 {11, 0, 0, 0},
510 {-10, 0, 5, 0},
511 {-8, 0, 3, 0},
512 {7, 0, -3, 0},
513 {-7, 0, 0, 0},
514 {-7, 0, 3, 0},
515 {-7, 0, 3, 0},
516 {6, 0, 0, 0},
517 {6, 0, -3, 0},
518 {6, 0, -3, 0},
519 {-6, 0, 3, 0},
520 {-6, 0, 3, 0},
521 {5, 0, 0, 0},
522 {-5, 0, 3, 0},
523 {-5, 0, 3, 0},
524 {-5, 0, 3, 0},
525 {4, 0, 0, 0},
526 {4, 0, 0, 0},
527 {4, 0, 0, 0},
528 {-4, 0, 0, 0},
529 {-4, 0, 0, 0},
530 {-4, 0, 0, 0},
531 {3, 0, 0, 0},
532 {-3, 0, 0, 0},
533 {-3, 0, 0, 0},
534 {-3, 0, 0, 0},
535 {-3, 0, 0, 0},
536 {-3, 0, 0, 0},
537 {-3, 0, 0, 0},
538 {-3, 0, 0, 0}
539 };
540
541 const int COEFF[][5]={
542 {0, 0, 0, 0, 1},
543 {-2, 0, 0, 2, 2},
544 {0, 0, 0, 2, 2},
545 {0, 0, 0, 0, 2},
546 {0, 1, 0, 0, 0},
547 {0, 0, 1, 0, 0},
548 {-2, 1, 0, 2, 2},
549 {0, 0, 0, 2, 1},
550 {0, 0, 1, 2, 2},
551 {-2, -1, 0, 2, 2},
552 {-2, 0, 1, 0, 0},
553 {-2, 0, 0, 2, 1},
554 {0, 0, -1, 2, 2},
555 {2, 0, 0, 0, 0},
556 {0, 0, 1, 0, 1},
557 {2, 0, -1, 2, 2},
558 {0, 0, -1, 0, 1},
559 {0, 0, 1, 2, 1},
560 {-2, 0, 2, 0, 0},
561 {0, 0, -2, 2, 1},
562 {2, 0, 0, 2, 2},
563 {0, 0, 2, 2, 2},
564 {0, 0, 2, 0, 0},
565 {-2, 0, 1, 2, 2},
566 {0, 0, 0, 2, 0},
567 {-2, 0, 0, 2, 0},
568 {0, 0, -1, 2, 1},
569 {0, 2, 0, 0, 0},
570 {2, 0, -1, 0, 1},
571 {-2, 2, 0, 2, 2},
572 {0, 1, 0, 0, 1},
573 {-2, 0, 1, 0, 1},
574 {0, -1, 0, 0, 1},
575 {0, 0, 2, -2, 0},
576 {2, 0, -1, 2, 1},
577 {2, 0, 1, 2, 2},
578 {0, 1, 0, 2, 2},
579 {-2, 1, 1, 0, 0},
580 {0, -1, 0, 2, 2},
581 {2, 0, 0, 2, 1},
582 {2, 0, 1, 0, 0},
583 {-2, 0, 2, 2, 2},
584 {-2, 0, 1, 2, 1},
585 {2, 0, -2, 0, 1},
586 {2, 0, 0, 0, 1},
587 {0, -1, 1, 0, 0},
588 {-2, -1, 0, 2, 1},
589 {-2, 0, 0, 0, 1},
590 {0, 0, 2, 2, 1},
591 {-2, 0, 2, 0, 1},
592 {-2, 1, 0, 2, 1},
593 {0, 0, 1, -2, 0},
594 {-1, 0, 1, 0, 0},
595 {-2, 1, 0, 0, 0},
596 {1, 0, 0, 0, 0},
597 {0, 0, 1, 2, 0},
598 {0, 0, -2, 2, 2},
599 {-1, -1, 1, 0, 0},
600 {0, 1, 1, 0, 0},
601 {0, -1, 1, 2, 2},
602 {2, -1, -1, 2, 2},
603 {0, 0, 3, 2, 2},
604 {2, -1, 0, 2, 2}
605 };
606
getAstroValuesByDay(double julianDay,const Location * loc,Astro * astro,Astro * topAstro)607 void getAstroValuesByDay(double julianDay, const Location* loc, Astro* astro,
608 Astro* topAstro)
609 {
610 AstroDay ad;
611
612 if (astro->jd == julianDay-1)
613 {
614 /* Copy cached values */
615 astro->ra[0] = astro->ra[1];
616 astro->ra[1] = astro->ra[2];
617 astro->dec[0] = astro->dec[1];
618 astro->dec[1] = astro->dec[2];
619 astro->sid[0] = astro->sid[1];
620 astro->sid[1] = astro->sid[2];
621 astro->dra[0] = astro->dra[1];
622 astro->dra[1] = astro->dra[2];
623 astro->rsum[0] = astro->rsum[1];
624 astro->rsum[1] = astro->rsum[2];
625 /* Compute next day values */
626 computeAstroDay(julianDay+1, &ad);
627 astro->ra[2] = ad.ra;
628 astro->dec[2] = ad.dec;
629 astro->sid[2] = ad.sidtime;
630 astro->dra[2] = ad.dra;
631 astro->rsum[2] = ad.rsum;
632 }
633 else if (astro->jd == julianDay + 1)
634 {
635 /* Copy cached values */
636 astro->ra[2] = astro->ra[1];
637 astro->ra[1] = astro->ra[0];
638 astro->dec[2] = astro->dec[1];
639 astro->dec[1] = astro->dec[0];
640 astro->sid[2] = astro->sid[1];
641 astro->sid[1] = astro->sid[0];
642 astro->dra[2] = astro->dra[1];
643 astro->dra[1] = astro->dra[0];
644 astro->rsum[2] = astro->rsum[1];
645 astro->rsum[1] = astro->rsum[0];
646 /* Compute previous day values */
647 computeAstroDay(julianDay-1, &ad);
648 astro->ra[0] = ad.ra;
649 astro->dec[0] = ad.dec;
650 astro->sid[0] = ad.sidtime;
651 astro->dra[0] = ad.dra;
652 astro->rsum[0] = ad.rsum;
653
654
655 } else if (astro->jd != julianDay)
656 {
657 /* Compute 3 day values */
658 computeAstroDay(julianDay-1, &ad);
659 astro->ra[0] = ad.ra;
660 astro->dec[0] = ad.dec;
661 astro->sid[0] = ad.sidtime;
662 astro->dra[0] = ad.dra;
663 astro->rsum[0] = ad.rsum;
664 computeAstroDay(julianDay, &ad);
665 astro->ra[1] = ad.ra;
666 astro->dec[1] = ad.dec;
667 astro->sid[1] = ad.sidtime;
668 astro->dra[1] = ad.dra;
669 astro->rsum[1] = ad.rsum;
670 computeAstroDay(julianDay+1, &ad);
671 astro->ra[2] = ad.ra;
672 astro->dec[2] = ad.dec;
673 astro->sid[2] = ad.sidtime;
674 astro->dra[2] = ad.dra;
675 astro->rsum[2] = ad.rsum;
676
677 }
678
679 astro->jd = julianDay;
680 computeTopAstro(loc, astro, topAstro);
681
682 }
683
684
computeAstroDay(double JD,AstroDay * astroday)685 void computeAstroDay(double JD, AstroDay* astroday)
686 {
687
688 int i =0;
689 double R, Gg, rGg, G;
690
691 double tL, L;
692 double tB, B;
693
694 double D, M, M1, F, O;
695
696 double U, E0, E, rE, lambda, rLambda, V0, V;
697
698 double RAn, RAd, RA, DEC;
699
700 double B0sum=0, B1sum=0;
701 double R0sum=0, R1sum=0, R2sum=0, R3sum=0, R4sum=0;
702 double L0sum=0, L1sum=0, L2sum=0, L3sum=0, L4sum=0, L5sum=0;
703
704 double PNsum=0, psi=0, epsilon=0;
705 double deltaPsi, deltaEps;
706
707 double JC = (JD - 2451545)/36525.0;
708 double JM = JC/10.0;
709 double JM2 = pow (JM, 2);
710 double JM3 = pow (JM, 3);
711 double JM4 = pow (JM, 4);
712 double JM5 = pow (JM, 5);
713
714 /* FIXIT: By default, the getJulianDay function returns JDE rather then JD,
715 * make sure this is accurate, and works in last-day-of-year
716 * circumstances. */
717 double JDE = JD;
718
719 double T = (JDE - 2451545)/36525.0;
720
721 for(i=0; i < 64; i++)
722 L0sum += L0[i][0] * cos(L0[i][1] + L0[i][2] * JM);
723 for(i=0; i < 34; i++)
724 L1sum += L1[i][0] * cos(L1[i][1] + L1[i][2] * JM);
725 for(i=0; i < 20; i++)
726 L2sum += L2[i][0] * cos(L2[i][1] + L2[i][2] * JM);
727 for(i=0; i < 7; i++)
728 L3sum += L3[i][0] * cos(L3[i][1] + L3[i][2] * JM);
729 for(i=0; i < 3; i++)
730 L4sum += L4[i][0] * cos(L4[i][1] + L4[i][2] * JM);
731 L5sum = L5[0][0] * cos(L5[0][1] + L5[0][2] * JM);
732
733
734 tL = (L0sum + (L1sum * JM) + (L2sum * JM2)
735 + (L3sum * JM3) + (L4sum * JM4)
736 + (L5sum * JM5)) / pow (10, 8);
737
738 L = limitAngle(RAD_TO_DEG(tL));
739
740 for(i=0; i<5; i++)
741 B0sum += B0[i][0] * cos(B0[i][1] + B0[i][2] * JM);
742 for(i=0; i<2; i++)
743 B1sum += B1[i][0] * cos(B1[i][1] + B1[i][2] * JM);
744
745
746 tB= (B0sum + (B1sum * JM)) / pow (10, 8);
747 B = RAD_TO_DEG(tB);
748
749
750 for(i=0; i < 40; i++)
751 R0sum += R0[i][0] * cos(R0[i][1] + R0[i][2] * JM);
752 for(i=0; i < 10; i++)
753 R1sum += R1[i][0] * cos(R1[i][1] + R1[i][2] * JM);
754 for(i=0; i < 6; i++)
755 R2sum += R2[i][0] * cos(R2[i][1] + R2[i][2] * JM);
756 for(i=0; i < 2; i++)
757 R3sum += R3[i][0] * cos(R3[i][1] + R3[i][2] * JM);
758 R4sum = R4[0][0] * cos(R4[0][1] + R4[0][2] * JM);
759
760 R = (R0sum + (R1sum * JM) + (R2sum * JM2)
761 + (R3sum * JM3) + (R4sum * JM4)) / pow (10, 8);
762
763 G = limitAngle((L + 180));
764 Gg = -B;
765 rGg = DEG_TO_RAD(Gg);
766 /* Compute the fundamental arguments (p. 144) */
767 D = 297.85036 + (445267.111480 * T) - (0.0019142 * pow (T, 2)) +
768 (pow (T, 3)/189474.0);
769 M = 357.52772 + (35999.050340 * T) - (0.0001603 * pow (T, 2)) -
770 (pow (T, 3)/300000.0);
771 M1 = 134.96298 + (477198.867398 * T) + (0.0086972 * pow (T, 2)) +
772 (pow (T, 3)/56250.0);
773 F = 93.27191 + (483202.017538 * T) - ( 0.0036825 * pow (T, 2)) +
774 (pow (T, 3)/327270.0);
775 O = 125.04452 - (1934.136261 * T) + (0.0020708 * pow (T, 2)) +
776 (pow (T, 3)/450000.0);
777 /* Add the terms (pp. 144-6) */
778 for (i=0; i<63; i++) {
779 PNsum += D * COEFF[i][0];
780 PNsum += M * COEFF[i][1];
781 PNsum += M1 * COEFF[i][2];
782 PNsum += F * COEFF[i][3];
783 PNsum += O * COEFF[i][4];
784 psi += (PN[i][0] + JC*PN[i][1])*sin(DEG_TO_RAD(PNsum));
785 epsilon += (PN[i][2] + JC*PN[i][3])*cos(DEG_TO_RAD(PNsum));
786 PNsum=0;
787 }
788
789 deltaPsi = psi/36000000.0;
790 /* Nutation in obliquity */
791 deltaEps = epsilon/36000000.0;
792
793
794 /* The obliquity of the ecliptic (p. 147, 22.3) */
795 U = JM/10.0;
796 E0 = 84381.448 - 4680.93 * U - 1.55 * pow(U,2) + 1999.25 * pow(U,3)
797 - 51.38 * pow(U,4) - 249.67 * pow(U,5) - 39.05 * pow(U,6) + 7.12
798 * pow(U,7) + 27.87 * pow(U,8) + 5.79 * pow(U,9) + 2.45 * pow(U,10);
799 /* Real/true obliquity (p. 147) */
800 E = E0/3600.0 + deltaEps;
801 rE = DEG_TO_RAD(E);
802
803 lambda = G + deltaPsi + (-20.4898/(3600.0 * R));
804 rLambda = DEG_TO_RAD(lambda);
805
806 /* Mean Sidereal time (p. 88) */
807 V0 = 280.46061837 + 360.98564736629 * ( JD - 2451545) +
808 0.000387933 * pow(JC,2) - pow(JC,3)/ 38710000.0;
809 /* Apparent sidereal time */
810 V = limitAngle(V0) + deltaPsi * cos(rE);
811
812 RAn = sin(rLambda) * cos(rE) - tan(rGg) * sin(rE);
813 RAd = cos(rLambda);
814 RA = limitAngle(RAD_TO_DEG(atan2(RAn,RAd)));
815
816 DEC = asin( sin(rGg) * cos(rE) + cos(rGg) * sin(rE) *
817 sin(rLambda));
818
819 astroday->ra = RA;
820 astroday->dec = DEC;
821 astroday->sidtime = V;
822 astroday->dra = 0;
823 astroday->rsum = R;
824
825 }
826
computeTopAstro(const Location * loc,const Astro * astro,Astro * topAstro)827 void computeTopAstro(const Location* loc, const Astro* astro, Astro* topAstro)
828 {
829 int i;
830 double lHour, SP, rlHour, rLat;
831 double tU, tpCos, tpSin, tRA0 ,tRA ,tDEC;
832
833 rLat = DEG_TO_RAD(loc->degreeLat);
834
835 for (i=0; i<3; i++)
836 {
837 lHour = limitAngle(astro->sid[i] + loc->degreeLong - astro->ra[i]);
838 rlHour = DEG_TO_RAD(lHour);
839
840 SP = DEG_TO_RAD (8.794/(3600 * astro->rsum[i]));
841
842 /* (p. 82, with b/a = 0.99664719) */
843 tU = atan (0.99664719 * tan(rLat));
844
845 tpSin = 0.99664719 * sin(tU) + (loc->seaLevel/EARTH_RADIUS) *
846 sin(rLat);
847
848 tpCos = cos(tU) + (loc->seaLevel/EARTH_RADIUS) * cos(rLat);
849
850 /* (p. 297, 40.2) */
851 tRA0 = (((-tpCos) * sin(SP) * sin(rlHour)) / (cos(astro->dec[i]) -
852 tpCos * sin(SP) * cos(rlHour)));
853
854 tRA = astro->ra[i] + RAD_TO_DEG(tRA0);
855
856 /* (p. 297, 40.3) */
857 tDEC = RAD_TO_DEG(atan2((sin(astro->dec[i]) - tpSin * sin(SP)) * cos(tRA0),
858 cos(astro->dec[i]) - tpCos * sin(SP) *
859 cos(rlHour)));
860
861 topAstro->ra[i] = tRA;
862 topAstro->dec[i] = tDEC;
863 topAstro->sid[i] = astro->sid[i];
864 topAstro->dra[i] = tRA0;
865 topAstro->rsum[i] = astro->rsum[i];
866
867 }
868
869 }
870
limitAngle(double L)871 static double limitAngle(double L)
872 {
873 double F;
874 L /= 360.0;
875 F = L - floor(L);
876 if (F > 0)
877 return 360 * F;
878 else if (F < 0)
879 return 360 - 360 * F;
880 else return L;
881 }
882
883
limitAngle180(double L)884 static double limitAngle180(double L)
885 {
886 double F;
887 L /= 180.0;
888 F = L - floor(L);
889 if (F > 0)
890 return 180 * F;
891 else if (F < 0)
892 return 180 - 180 * F;
893 else return L;
894 }
895
896 /* Limit between 0 and 1 (fraction of day)*/
limitAngle1(double L)897 static double limitAngle1(double L)
898 {
899 double F;
900 F = L - floor(L);
901 if (F < 0)
902 return F += 1;
903 return F;
904 }
905
limitAngle180between(double L)906 static double limitAngle180between(double L)
907 {
908 double F;
909 L /= 360.0;
910 F = (L - floor(L)) * 360.0;
911 if (F < -180)
912 F += 360;
913 else if (F > 180)
914 F -= 360;
915 return F;
916 }
917
918