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