1 #include "prayer.h"
2
3 calc_method_t calc_methods [] = {
4 { MWL, "Muslim World League (MWL)", ANGLE, 18, 17 },
5 { ISNA, "Islamic Society of North America (ISNA)", ANGLE, 15, 15 },
6 { EGAS, "Egyptian General Authority of Survey", ANGLE, 19.5, 17.5 },
7 { UMAQ, "Umm Al-Qura University, Makkah", OFFSET, 18.5, 90 },
8 { UIS, "University of Islamic Sciences, Karachi", ANGLE, 18, 18 }
9 };
10
11 /* Helper geometric and mathematical funtions */
12 static double to_degrees(const double x);
13 static double to_radians(const double x);
14 static double arccot(const double x);
15 /* round is available only in C99 */
16 static double custom_round(const double x);
17 /* Normalizes the given value x within the range [0,N] */
18 static double normalize(const double x, const double N);
19
20 /* Get the Julian Day number of a given date */
21 static unsigned long get_julian_day_number(const struct tm *date);
22
23 /* Computes the approximate Sun coordinates for
24 the given Julian Day Number jdn */
25 static void get_approx_sun_coord(const unsigned long jdn,
26 struct approx_sun_coord *coord);
27
28 /* T function used to compute Sunrise, Sunset, Fajr and Isha
29 Taken from: http://praytimes.org/calculation/
30 */
31 static double T(const double alpha,
32 const double latitude,
33 const double D);
34
35 /* A function used to compute Asr
36 Taken from: http://praytimes.org/calculation/
37 */
38 static double A(const double t,
39 const double latitude,
40 const double D);
41
42 /* Methods for individual times and prayers */
43 static double get_dhuhr(const struct location *loc,
44 const struct approx_sun_coord *coord);
45
46 static double get_sunrise(const double dhuhr,
47 const struct location *loc,
48 const struct approx_sun_coord *coord);
49
50 static double get_fajr(const double dhuhr,
51 const double sunset,
52 const double sunrise,
53 const struct location *loc,
54 const struct approx_sun_coord *coord);
55
56 static double get_isha(const double dhuhr,
57 const double sunset,
58 const double sunrise,
59 const struct location *loc,
60 const struct approx_sun_coord *coord);
61
62 static void conv_time_to_event(const unsigned long julian_day,
63 const double decimal_time,
64 const round_t round,
65 event_t *event);
66
67
68
69 /**
70 * Convert a given angle specified in degrees into radians
71 */
to_radians(const double x)72 static double to_radians(const double x)
73 {
74 return ((x) * M_PI) / 180.0;
75 }
76
77
78 /**
79 * Convert a given angle specified in radians into degrees
80 */
to_degrees(const double x)81 static double to_degrees(const double x)
82 {
83 return ((x) * 180.0) / M_PI;
84 }
85
86
87 /**
88 * Compute the arccotangent of a given value
89 */
arccot(const double x)90 static double arccot(const double x)
91 {
92 return atan2(1.0, (x));
93 }
94
95
96 /**
97 * Normalizes the given value x to be within the range [0,N]
98 */
normalize(const double x,const double N)99 static double normalize(const double x, const double N)
100 {
101 double n;
102 assert(N > 0.0);
103 n = x - (N * floor( x / N ));
104 assert(n <= N);
105 return n;
106 }
107
108 /**
109 * round is avalilable only in C99. Therefore, we use a custom one
110 * @return The rounded value of x
111 */
custom_round(const double x)112 static double custom_round(const double x)
113 {
114 return floor(x + 0.5);
115 }
116
117 /**
118 * Compute the Julian Day Number for a given date
119 * Source: https://en.wikipedia.org/wiki/Julian_day
120 * @return The Julian day number
121 */
get_julian_day_number(const struct tm * date)122 static unsigned long get_julian_day_number(const struct tm *date)
123 {
124 double a, y, m, jdn;
125
126 assert (date != NULL);
127
128 a = floor((14.0 - (double)(date->tm_mon + 1)) / 12.0);
129 y = (double)(1900 + date->tm_year) + 4800.0 - a;
130 m = (double)(date->tm_mon + 1) + 12.0 * a - 3.0;
131 jdn = (double)date->tm_mday + \
132 floor((153.0 * m + 2.0) / 5.0) + \
133 365.0 * y + \
134 floor(y / 4.0) - \
135 floor(y / 100.0) + \
136 floor(y / 400.0) - \
137 32045.0;
138
139 return (unsigned long)jdn;
140 }
141
142
143 /**
144 * Compute the approximate Sun coordinates from the given Julian Day
145 * Number jdn according to the method given the US Naval Observatory (USNO)
146 * on: http://aa.usno.navy.mil/faq/docs/SunApprox.php
147 * The computed values are stored in struct coord.
148 */
get_approx_sun_coord(const unsigned long jdn,struct approx_sun_coord * coord)149 static void get_approx_sun_coord(const unsigned long jdn,
150 struct approx_sun_coord *coord)
151 {
152 double d; /* Julian Day Number from 1 January 2000 */
153 double g; /* Mean anomaly of the Sun */
154 double q; /* Mean longitude of the Sun */
155 double L; /* Geocentric apparent ecliptic longitude of the Sun
156 (adjusted for aberration) */
157 double R; /* The approximated distance of the Sun from the Earth
158 in astronomical units (AU) */
159 double e; /* The mean obliquity of the ecliptic */
160 double RA; /* The Sun's right ascension */
161 double D; /* The Sun's Declination */
162 double EqT; /* The Equation of Time */
163 double SD; /* The angular semidiameter of the Sun in degrees */
164
165 assert (coord != NULL);
166 assert (jdn > 2451545);
167
168 d = (double)(jdn - 2451545); /* Remove the offset from jdn */
169 g = 357.529 + 0.98560028 * d;
170 q = 280.459 + 0.98564736 * d;
171 L = q + 1.915 * sin(to_radians(g)) + \
172 0.020 * sin(to_radians(2.0*g));
173 R = 1.00014 - 0.01671 * cos(to_radians(g)) - \
174 0.00014 * cos(to_radians(2.0*g));
175 e = 23.439 - 0.00000036 * d;
176 RA = to_degrees(atan2(cos(to_radians(e)) * sin(to_radians(L)), \
177 cos(to_radians(L)))/15.0);
178 D = to_degrees(asin(sin(to_radians(e)) * sin(to_radians(L))));
179 EqT = q/15.0 - RA;
180
181 /* Resulting EqT Can be larger than 360.
182 Therefore, it needs normalization */
183 EqT = normalize(EqT, 360.0);
184 SD = 0.2666 / R;
185
186 coord->D = D;
187 coord->EqT = EqT;
188 coord->R = R;
189 coord->SD = SD;
190 }
191
192
193 /**
194 * T function used to compute Sunrise, Sunset, Fajr, and Isha
195 * Taken from: http://praytimes.org/calculation/
196 */
T(const double alpha,const double latitude,const double D)197 static double T(const double alpha,
198 const double latitude,
199 const double D)
200 {
201 double p1 = 1.0/15.0;
202 double p2 = cos(to_radians(latitude)) * cos(to_radians(D));
203 double p3 = sin(to_radians(latitude)) * sin(to_radians(D));
204 double p4 = -1.0 * sin(to_radians(alpha));
205 double p5 = to_degrees(acos((p4 - p3) / p2));
206 double r = p1 * p5;
207 return r;
208 }
209
210
211 /**
212 * A function used to compute Asr
213 * Taken from: http://praytimes.org/calculation/
214 */
A(const double t,const double latitude,const double D)215 static double A(const double t,
216 const double latitude,
217 const double D)
218 {
219 double p1 = 1.0/15.0;
220 double p2 = cos(to_radians(latitude)) * cos(to_radians(D));
221 double p3 = sin(to_radians(latitude)) * sin(to_radians(D));
222 double p4 = tan(to_radians((latitude - D)));
223 double p5 = arccot(t + p4);
224 double p6 = sin(p5);
225 double p7 = acos((p6 - p3) / p2);
226 double r = p1 * to_degrees(p7);
227 return r;
228 }
229
230
231 /**
232 * Compute the Dhuhr prayer time
233 */
get_dhuhr(const struct location * loc,const struct approx_sun_coord * coord)234 static double get_dhuhr(const struct location *loc,
235 const struct approx_sun_coord *coord)
236 {
237 double dhuhr = 0.0;
238
239 assert(loc != NULL);
240 assert(coord != NULL);
241
242 dhuhr = 12.0 + loc->timezone - loc->longitude/15.0 - coord->EqT;
243 dhuhr = normalize(dhuhr, 24.0);
244 return dhuhr;
245 }
246
247
248 /**
249 * Compute the Asr prayer time
250 */
get_asr(const double dhuhr,const struct location * loc,const struct approx_sun_coord * coord)251 static double get_asr(const double dhuhr,
252 const struct location *loc,
253 const struct approx_sun_coord *coord)
254 {
255 double asr = 0.0;
256
257 assert(dhuhr > 0.0);
258 assert(loc != NULL);
259 assert(coord != NULL);
260 assert(loc->asr_method == SHAFII || loc->asr_method == HANAFI);
261
262 switch(loc->asr_method) {
263 case SHAFII:
264 asr = dhuhr + A(1.0, loc->latitude, coord->D);
265 break;
266 case HANAFI:
267 asr = dhuhr + A(2.0, loc->latitude, coord->D);
268 break;
269 }
270 asr = normalize(asr, 24.0);
271 return asr;
272 }
273
274
275 /**
276 * Compute the Sunrise time
277 * TODO: Check if this is valid for extreme altitudes
278 */
get_sunrise(const double dhuhr,const struct location * loc,const struct approx_sun_coord * coord)279 static double get_sunrise(const double dhuhr,
280 const struct location *loc,
281 const struct approx_sun_coord *coord)
282 {
283 assert(loc != NULL);
284 assert(coord != NULL);
285
286 return (dhuhr - T(0.8333 + 0.0347 * sqrt(loc->altitude),\
287 loc->latitude, coord->D));
288 }
289
290
291 /**
292 * Compute the Sunset time
293 * TODO: Check if this is valid for extreme altitudes
294 */
get_sunset(const double dhuhr,const struct location * loc,const struct approx_sun_coord * coord)295 static double get_sunset(const double dhuhr,
296 const struct location *loc,
297 const struct approx_sun_coord *coord)
298 {
299 assert(loc != NULL);
300 assert(coord != NULL);
301
302 return (dhuhr + T(0.8333 + 0.0347 * sqrt(loc->altitude),\
303 loc->latitude, coord->D));
304 }
305
306
307 /**
308 * Compute the Fajr prayer time
309 */
get_fajr(const double dhuhr,const double sunset,const double sunrise,const struct location * loc,const struct approx_sun_coord * coord)310 static double get_fajr(const double dhuhr,
311 const double sunset,
312 const double sunrise,
313 const struct location *loc,
314 const struct approx_sun_coord *coord)
315 {
316 double fajr = 0.0;
317
318 (void)sunset;
319 (void)sunrise;
320
321 assert(dhuhr > 0.0);
322 assert(loc != NULL);
323 assert(coord != NULL);
324
325 if (loc->extr_method == NONE) {
326 /* Normal latitude. Use the classical methods */
327 fajr = dhuhr - T(loc->calc_method.fajr,\
328 loc->latitude,\
329 coord->D);
330 } else {
331 /* TODO: Extreme methods */
332 fprintf(stderr, "Extreme methods are not implemented yet\n");
333 exit(EXIT_FAILURE);
334 }
335
336 return fajr;
337 }
338
339
340 /**
341 * Compute the Isha prayer time
342 */
get_isha(const double dhuhr,const double sunset,const double sunrise,const struct location * loc,const struct approx_sun_coord * coord)343 static double get_isha(const double dhuhr,
344 const double sunset,
345 const double sunrise,
346 const struct location *loc,
347 const struct approx_sun_coord *coord)
348 {
349 double isha = 0.0;
350
351 (void)sunset;
352 (void)sunrise;
353
354 assert(dhuhr > 0.0);
355 assert(loc != NULL);
356 assert(coord != NULL);
357
358 if (loc->extr_method == NONE) {
359 if (loc->calc_method.isha_type == OFFSET) {
360 /* Umm Al-Qura uses a fixed offset */
361 isha = sunset + loc->calc_method.isha * ONE_MINUTE;
362 } else {
363 isha = dhuhr + T(loc->calc_method.isha,\
364 loc->latitude,\
365 coord->D);
366 }
367 } else {
368 /* TODO: Extreme latitude */
369 fprintf(stderr, "Extreme altitudes are not implemented yet!\n");
370 exit(EXIT_FAILURE);
371 }
372
373 isha = normalize(isha, 24.0);
374 return isha;
375 }
376
377
378 /**
379 * Convert a given time in decimal format to
380 * hh:mm format and store it in the given event struct.
381 * The minutes part can be rounded up or down based on the
382 * round flag
383 */
conv_time_to_event(const unsigned long julian_day,const double decimal_time,const round_t rounding,event_t * t)384 static void conv_time_to_event(const unsigned long julian_day,
385 const double decimal_time,
386 const round_t rounding,
387 event_t *t)
388 {
389 double r = 0.0, f = 0.0;
390
391 assert(julian_day > 0);
392 assert(decimal_time >= 0.0 && decimal_time <= 24.0);
393 assert(t != NULL);
394 assert(rounding == UP ||
395 rounding == DOWN ||
396 rounding == NEAREST);
397
398 t->julian_day = julian_day;
399 f = floor(decimal_time);
400 t->hour = (unsigned int)f;
401 r = (decimal_time - f) * 60.0;
402 switch (rounding) {
403 case UP:
404 t->minute = (unsigned int)(ceil(r));
405 break;
406 case DOWN:
407 t->minute = (unsigned int)(floor(r));
408 break;
409 case NEAREST:
410 t->minute = (unsigned int)(custom_round(r));
411 break;
412 }
413 }
414
415
416 /**
417 * Compute the Qibla direction from North clockwise
418 * using the Equation (1) given in references/qibla.pdf
419 */
get_qibla_direction(const struct location * loc)420 double get_qibla_direction(const struct location *loc)
421 {
422 double p1, p2, p3, qibla;
423
424 assert(loc != NULL);
425
426 p1 = sin(to_radians(KAABA_LONGITUDE - loc->longitude));
427 p2 = cos(to_radians(loc->latitude)) * \
428 tan(to_radians(KAABA_LATITUDE));
429 p3 = sin(to_radians(loc->latitude)) * \
430 cos(to_radians(KAABA_LONGITUDE - loc->longitude));
431 qibla = to_degrees(atan2(p1, (p2 - p3)));
432
433 return qibla;
434 }
435
436
get_prayer_times(const struct tm * date,const struct location * loc,struct prayer_times * pt)437 void get_prayer_times(const struct tm *date,
438 const struct location *loc,
439 struct prayer_times *pt)
440 {
441 unsigned long jdn, jdn_next, jdn_prev;
442 double true_noon, sunrise, sunset;
443 double noon_next, sunrise_next;
444 double noon_prev, sunset_prev;
445 double fajr, dhuhr, asr, maghrib, isha;
446 struct approx_sun_coord coord, coord_next, coord_prev;
447
448 assert(date != NULL);
449 assert(loc != NULL);
450 assert(pt != NULL);
451
452 jdn = get_julian_day_number(date);
453 jdn_next = jdn + 1;
454 jdn_prev = jdn - 1;
455
456 get_approx_sun_coord(jdn, &coord);
457 get_approx_sun_coord(jdn_next, &coord_next);
458 get_approx_sun_coord(jdn_prev, &coord_prev);
459
460 true_noon = get_dhuhr(loc, &coord);
461 noon_next = get_dhuhr(loc, &coord_next);
462 noon_prev = get_dhuhr(loc, &coord_prev);
463
464 sunrise = get_sunrise(true_noon, loc, &coord);
465 sunset = get_sunset(true_noon, loc, &coord);
466 sunrise_next = get_sunrise(noon_next, loc, &coord_next);
467 sunset_prev = get_sunset(noon_prev, loc, &coord_prev);
468
469 fajr = get_fajr(true_noon, sunset_prev, sunrise, loc, &coord);
470 dhuhr = true_noon;
471 asr = get_asr(true_noon, loc, &coord);
472 maghrib = sunset;
473 isha = get_isha(true_noon, sunset, sunrise_next, loc, &coord);
474
475 /* TODO: Find what Fiqh says regarding the rounding... */
476 conv_time_to_event(jdn, fajr, NEAREST, &(pt->fajr));
477 conv_time_to_event(jdn, sunrise, NEAREST, &(pt->sunrise));
478 conv_time_to_event(jdn, dhuhr, NEAREST, &(pt->dhuhr));
479 conv_time_to_event(jdn, asr, NEAREST, &(pt->asr));
480 conv_time_to_event(jdn, maghrib, NEAREST, &(pt->maghrib));
481 conv_time_to_event(jdn, isha, NEAREST, &(pt->isha));
482 }
483