unit PTime; {* praytime.c Version 2.1 * * Program to compute Islamic prayer hour schedules. * This program is a modified version of Kamal Abdali's 'praytimer' * which produces TeX code for the schedule on the standard output. * 'praytimer' includes the following notice: * * Copyright (c) 1987--1992 Kamal Abdali * * Permission for nonprofit use of this software and its documentation * is hereby granted without fee, provided that the above copyright notice * appear in all copies and that both that copyright notice and this * permission notice appear in supporting documentation, and that the name * of Kamal Abdali not be used in advertising or publicity pertaining to * distribution of the software without specific, written prior permission. * * * Modified by Waleed A. Muhanna (wmuhanna@magnus.acs.ohio-state.edu) to * improve input checking, remove all TeX related stuff, allow a default * location to be easily configured in, and permit the user to print a * schedule for the current date or a specified month and year. * * Please send any comments/suggestions/fixes/additions to * wmuhanna@magnus.acs.ohio-state.edu * *} {* * The program generates Islamic prayer time schedules for any * location. If no argument is given, the prayer time schedule * for the current date at the current (configured) location is * produced. If a year is specified, the program outputs a schedule * for the entire year. The year can be between 1900 and 2200, or 0 * for a 'perpetual' schedule. If a month is also specified, a schedule * just for that month is printed. * * The following command line arguments control the calculation method: * * -a angle Sun's angle of depression at Fajr in degrees (usually 18) * -t time Time interval from Fajr to sunrise in minutes (usually 90) * -f fiqh Value should be S(hafii) or H(anafi) * -r ratio Shadow ratio at Asr. (Usually 1 for SHafii, 2 for others) * * NOTE: It is an error to specify both -a and -t or both -f and -r. * DEFAULT is as if called with '-a 18 -r 1' * * To generate schedules for locations other than the default (configured) * location, the following command line option should be indicated. * * -i reads name and geographical data for the location from * the standard input, instead of using the default (configured) * location. If standard input is a terminal, the program * prompts the user for the values. * * Data on standard input must contain (in given order) * Name of location (upto 40 caharacters) * Latitude degrees and minutes, and N or S to specify north or south * Longitude degrees and minutes, and E or W to specify east or west * Time Zone in hours (Decimal for fractional hour zones, negative if * West of Greenwich) * Y or 1 if DayLight Saving Time adjustment needed. N or 0, otherwise. * * Data items should be separated by whitespace, but the name must be on * a separate line by itself because it may contain spaces or punctuation. * * Input may contain data for more than one location. * An interactive session may be ended by typing the End-of-File character * (e.g. CTRL-D in UNIX) when prompted for Name. *} interface uses hConv, Math, SysUtils, Classes; type rLoc = record Name: string; LatD, LatM: integer; LatO: char; LongD, LongM: integer; LongO: char; TimeZ: byte; DayLight: boolean; end; rPTime = record QiblahD, QiblahM: integer; QiblahO: char; Day, DoW: byte; Mark: char; Fajr, Shuruq, Dhuhr, Asr, Maghrib, Isha: string[5]; end; var myrLoc: rLoc; myrPTime: rPTime; procedure GetPrayTime(f: char; d, m, y: integer); overload; procedure GetPrayTime(f: char; dt: TDateTime); overload; procedure header(day,month,year: integer); procedure Display(first, last, startdate, cd, cm, header: integer); procedure getDatas; procedure makeSchedule(day,month,year: integer); procedure computeHours(year, first, last: integer); procedure computeConstants(year: integer); function Qibla: double; function NoonTime(nday: integer; var coaltn: double): double; function tempus(nday: integer; coalt, time0: double): double; procedure DayLight(year, hasDayLt: integer; var leap, beg, finish: integer); function deg2rad(degree: double): double; function dm2deg(degree, minute: integer): double; function dms2deg(degree: longint; min: integer; sec: double): double; procedure deg2dm(degree: double; var deg, min: integer); function hms2h(hour, min: integer; sec: double): double; function FABS(a: variant): variant; function FMOD(a,b: variant): variant; implementation const PTLOC = 'Makkah Al-Mukarramah 21 25 N 39 49 E 3 N'; MAXNAMEL = 40; DPR = (57.29577951308230876799);{* degree per radian (180/pi) *} RPD = (0.01745329251994329577); {* radians per degree (pi/180) *} HPR = (3.81971863420548805845); {* hours per radian (12/pi) *} MAXENV = 9; dowl: array[0..6] of string = ('Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'); dir: array[0..3] of string = ('E', 'W', 'N', 'S'); var leap, beginDayLight, endDayLight: integer; latd, latm, longd, longm, latIsS, longIsW, fajrByInterval, fajrByMaxInterval: integer; tim: array[0..366,0..8] of real; asrShadowRatio, fajrInterval, fajrDepr: double; cosobl, sinobl, perigee, dmanom, anom0: double; c1, c2, delsid, sidtm0: double; hd: sdate; name: string[MAXNAMEL+2]; latitude, longitude, timeZone: double; hasDayLt: integer; hflag: integer; stdow: integer; function FABS(a: variant): variant; begin if ( a > 0 ) then Result := a else Result := -a; end; function FMOD(a,b: variant): variant; begin Result := (a-(floor(a/b)*b)); end; procedure GetPrayTime(f: char; d, m, y: integer); var day, month, year: integer; fiqh: string; ratio, intvl: double; begin fiqh := UpperCase(f); intvl := 90.0; fajrByMaxInterval := 0; fajrByInterval := 0; fajrInterval := 1.5; fajrDepr := 18.0; asrShadowRatio := 1.0; if ( fiqh = 'S' ) then ratio := 1.0 else ratio := 2.0; asrShadowRatio := ratio; day := d; month := m; year := y; getDatas; makeSchedule(day,month,year); end; procedure GetPrayTime(f: char; dt: TDateTime); var day, month, year: word; fiqh: string; ratio, intvl: double; begin fiqh := UpperCase(f); intvl := 90.0; fajrByMaxInterval := 0; fajrByInterval := 0; fajrInterval := 1.5; fajrDepr := 18.0; asrShadowRatio := 1.0; if ( fiqh = 'S' ) then ratio := 1.0 else ratio := 2.0; asrShadowRatio := ratio; DecodeDate(dt,year,month,day); getDatas; makeSchedule(day,month,year); end; {* * Obtain name and geographical data for the location for which * the schedule is desired. If INTERACTIVE is true, then * the user is prompted for the info at the terminal. * Otherwise, the info is obtained from the standard input. * Note: time data and values to control the calculation method are * to be given as command line arguments, not as input. *} procedure getDatas; begin name := myrLoc.Name; latd := myrLoc.LatD; latm := myrLoc.LatM; if ( UpperCase(myrLoc.LatO) = 'S' ) then latIsS := 1 else latIsS := 0; longd := myrLoc.LongD; longm := myrLoc.LongM; if ( UpperCase(myrLoc.LongO) = 'W' ) then longIsW := 1 else longIsW := 0; timeZone := myrLoc.TimeZ; if ( myrLoc.DayLight ) then hasDayLt := 1 else hasDayLt := 0; latitude := deg2rad(dm2deg(latd,latm)); if ( latIsS = 1 ) then latitude := - latitude; longitude := deg2rad(dm2deg(longd,longm)); if ( longIsW = 0 ) then longitude := -longitude; end; {* * Computes day number (first is 0) in the year. *} function daynum(d,m,y: integer): integer; var i, dn: integer; begin i := 1; Dec(d); dn := d; while ( i < m ) do begin dn := dn+gDays(i,y); Inc(i); end; Result := dn; end; procedure setdow(y,m,d: integer); begin stdow := (Trunc(julianday(y,m,d, 0.0)+1.5)) mod 7; end; procedure hSchedule(m,y: integer); var fd, ld: sdate; cd, cm, ndays0, ndays1, sd, dh: integer; begin sd := 1; dh := 1; // mnthheader(m,y,sl); fd := gdate(y,m,1); ld := gdate(y, m+1, 1); ld := caldate(julianday(ld.Year, ld.Mon, ld.Day, 0.0)-1.0); ndays0 := daynum(fd.Day, fd.Mon, fd.Year); setdow(fd.Year,fd.Mon,fd.Day); cd := fd.Day; cm := fd.Mon; if ( fd.mon <> ld.mon ) then begin ndays1 := daynum(gDays(fd.Mon,fd.Year), fd.Mon, fd.Year)+1; computeHours(fd.Year,ndays0,ndays1); Display(ndays0, ndays1,sd,cd,cm, dh); Dec(dh); sd := ndays1-ndays0+1; ndays0 := daynum(1, ld.Mon, ld.Year); cd := 1; cm := ld.Mon; end; ndays1 := daynum(ld.Day, ld.Mon, ld.Year)+1; computeHours(ld.Year,ndays0,ndays1); Display(ndays0,ndays1, sd, cd, cm, dh); end; {* * Schedule computation section. * Prayer hours are computed basically following the algorithms given in * 'Prayer Schedules for North America', American Trust Publications, * Indianapolis, Indiana, 1978, Appendices A and B. * *} procedure makeSchedule(day,month,year: integer); var i: integer; ndays0, ndays1: integer; begin header(day,month,year); if ( day <> 0 ) then begin ndays0 := daynum(day,month,year); ndays1 := ndays0+1; computeHours(year,ndays0,ndays1); setdow(year,month,day); if ( hflag = 1 ) then Display(ndays0,ndays1,hd.day,day,month,1) else Display(ndays0,ndays1,day,day,month,1); end else if ( month <> 0 ) then begin if ( hflag = 0 ) then begin // mnthHeader(month,year); nDays0 := daynum(1,month,year); nDays1 := ndays0+gDays(month,year); computeHours(year,ndays0,ndays1); setDoW(year,month,1); Display(ndays0, ndays1, 1,0, 0, 1); end else hSchedule(month,year); end else begin if ( hflag = 0 ) then begin computeHours(year,0,365+leap); ndays0 := 0; setdow(year,1,1); for i:=1 to 12 do begin // mnthHeader(i,year); nDays1 := nDays0+gDays(i,year); Display(ndays0,ndays1,1,0, 0, 1); ndays0 := ndays1; end; end else for i :=1 to 12 do hSchedule(i,year); end; end; {* * Computes times for range of days first..last-1. *} procedure computeHours(year, first, last: integer); var coaltn: double; time0: array[0..6] of double; coalt: array[0..6] of double; t: double; i,k,l: integer; begin {* For perpetual, use 1994. (Need middle year of 4-year leap cycle *} if ( year = 0 ) then computeConstants(1994) else computeConstants(year); {* find beginning and ending days for DayLight saving time *} DayLight(year, hasDayLt, leap, beginDayLight, endDayLight); {* Approximate times of fajr,shuruq,asr,maghrib,isha *} time0[0] := 4.0; time0[1] := 6.0; time0[3] := 15.0; time0[4] := 18.0; time0[5] := 20.0; {* Coaltitudes of sun at fajr,shuruq,maghrib,isha *} /// coalt[0] := deg2rad((double)(90+fajrDepr)); coalt[0] := deg2rad(90+fajrDepr); coalt[1] := deg2rad(90.83); coalt[4] := coalt[1]; coalt[5] := coalt[0]; {* Get approximate times for the first day specified. *} {* Later on, each day's times used as approximate times for next day *} NoonTime(first,coaltn); coalt[3] := arctan(asrShadowRatio+tan(coaltn)); t := tempus(first,coalt[1],time0[1]); if ( t < 24.0 ) then time0[1] := t else time0[1] := 6.0; t := tempus(first,coalt[3],time0[3]); if ( t < 24.0 ) then time0[3] := t else time0[3] := 15.0; t := tempus(first,coalt[4],time0[4]); if ( t < 24.0 ) then time0[4] := t else time0[4] := 18.0; if ( fajrByInterval = 1 ) then begin time0[0] := time0[1]-fajrInterval; time0[5] := time0[4]+fajrInterval; end else begin t := tempus(first,coalt[0],time0[0]); if ( t < 24.0 ) then time0[0] := t else time0[0] := 4.0; t := tempus(first,coalt[5],time0[5]); if ( t < 24.0 ) then time0[5] := t else time0[5] := 20.0; end; {* compute times for the whole range of days *} for l := first to last-1 do begin {* for perpetual calendar, february 29 and march 1 have same times *} k := l; if ( (l > 59) and (year = 0) ) then k := l-1; tim[l,2] := NoonTime(k+1,coaltn); coalt[3] := ArcTan(asrShadowRatio+tan(coaltn)); t := tempus(k+1,coalt[1],time0[1]); tim[l][1] := t; if ( t < 24.0 ) then time0[1] := t else time0[1] := 6.0; t := tempus(k+1,coalt[3],time0[3]); tim[l,3] := t; if ( t < 24.0 ) then time0[3] := t else time0[3] := 15.0; t := tempus(k+1,coalt[4],time0[4]); tim[l,4] := t; if ( t < 24.0 ) then time0[4] := t else time0[4] := 18.0; if ( fajrByInterval = 1 ) then begin tim[l,0] := time0[1]-fajrInterval; time0[0] := time0[1]-fajrInterval; tim[l,5] := time0[4]+fajrInterval; time0[5] := time0[4]+fajrInterval; end else begin t := tempus(k+1,coalt[0],time0[0]); tim[l,0] := t; if ( t < 24.0 ) then time0[0] := t else time0[0] := 4.0; t := tempus(k+1,coalt[5],time0[5]); tim[l,5] := t; if ( t < 24.0 ) then time0[5] := t else time0[5] := 20.0; end; end; {* correct for DayLight saving time (if necessary) *} if ( (endDayLight <> 0) and (first < endDayLight) ) then begin i := beginDayLight-1; if ( first > i ) then i := first; while ( (i < endDayLight) and (i < last) ) do begin for k := 0 to 5 do tim[i,k] := tim[i,k]+1.0; Inc(i); end; end; end; {* * Computes astro constants for Jan 0 of given year *} procedure computeConstants(year: integer); {* ndays = time from 12 hr(noon), Jan 0, 1900 to 0 hr, Jan 0 of year *} {* t = same in julian centuries (units of 36525 days) *} {* obl = obliquity of ecliptic *} {* perigee = sun's longitude at perigee *} {* eccy = earth's eccentricity *} {* dmanom,delsid = daily motion (change) in *} {* sun's anomaly, sidereal time *} {* anom0,sidtm0 = sun's mean anomaly, *} {* sidereal time, all at 0 hr, jan 0 of year year *} {* c1,c2 = coefficients in equation of center *} var t: double; ndays: longint; obl, eccy: double; begin ndays := Trunc((Trunc(year-1900))*365+(year-1901)/4); t := (ndays-0.5)/36525.0; obl := deg2rad(dms2deg(23,27,8.26)-dms2deg(0,0,46.845)*t); cosobl := cos(obl); sinobl := sin(obl); eccy := 0.01675104-4.180e-5*t-1.26e-7*t*t; perigee := deg2rad(FMOD(dms2deg(281,13,15.0)+dms2deg(1,43,9.03)*t+ dms2deg(0,0,1.63)*t*t,360.0)); dmanom := deg2rad(dms2deg(35999,2,59.10)/36525.0); anom0 := deg2rad(FMOD(dms2deg(358,28,33.0)-dms2deg(0,0,0.54)*t*t+ FMOD(dms2deg(35999,2,59.10)*t,360.0),360.0)); delsid := hms2h(2400,3,4.542)/36525.0; sidtm0 := FMOD(hms2h(6,38,45.836)+FMOD(hms2h(2400,3,4.542)*t,24.0),24.0); c1 := eccy*(2-eccy*eccy/4); c2 := 5*eccy*eccy/4; end; {* * Double-duty function for leap year and DayLight Saving dates info. * Finds whether year is leap (sets leap = 1 if yes, 0 if no). * If hasDayLt is non-zero, then also computes the day numbers of the * start and end of DayLight Savings Time. * Sets begin = Day no. of the first Sunday of April, and * finish = Day no. of the Saturday before the last Sunday of October. *} procedure DayLight(year, hasDayLt: integer; var leap, beg, finish: integer); var m4,m1,jan0,napr1,noct31,apr1,oct31: integer; begin m4 := year mod 400; m1 := year mod 100; if ( (((year mod 4) = 0) and (m1 <> 0)) or (m4 = 0) ) then leap := 1 else leap := 0; if ( hasDayLt = 0 ) then begin {* No adjustment for DayLight Saving Time (year zero for perpetual) *} beg := 367; finish := 0; Exit; end; if ( year = 0 ) then begin {* DayLight Saving Time in perpetual calendar. April 1 thru Oct 31 *} beg := 92; {* April 1, 31+29+31+1 *} finish := beg+213; {* Oct 31, -1+30+31+30+31+31+30+31 *} Exit; end; {* Non-zero year. for annual calendar *} {* jan0,apr1,oct31 = day of week on those dates (fri=0,sat=1,sun=2,...) *} {* napr1,noct31 = Day no. in year on those dates *} jan0 := Trunc(m4/100*124+1+m1+m1/4-leap) mod 7; napr1 := 91+leap; {* 31+28+*leap+31+1 *} noct31 := 304+leap; {* 365+*leap-31-30 *} apr1 := (napr1+jan0) mod 7; oct31 := (noct31+jan0) mod 7; beg := napr1+2-apr1; if ( beg < napr1 ) then beg := beg+7; finish := noct31+2-oct31; if ( finish > noct31 ) then finish := finish-7; finish := finish-1; end; {* * Place sun's coaltitude at noon in coaltn, * and return time of noon for day no. nday of year *} function NoonTime(nday: integer; var coaltn: double): double; { {* slong = sun's true longitude at noon *} {* ra = sun's right ascension, decl = sun's declination *} {* ha = sun's hour angle west *} {* locmt = local mean time of phenomenon *} var t: double; longh,days,anomaly,slong,sinslong,ra,decl,locmt: double; begin longh := longitude*HPR; days := nday+(12.0+longh)/24.0; anomaly := anom0+dmanom*days; slong := perigee+anomaly+c1*sin(anomaly)+c2*sin(anomaly*2); sinslong := sin(slong); ra := ArcTan2(cosobl*sinslong,cos(slong))*HPR; if ( ra < 0.0 ) then ra := ra+24.0; decl := ArcSin(sinobl*sinslong); locmt := ra-delsid*days-sidtm0; t := locmt+longh+timeZone; if ( t < 0.0 ) then t := t+24.0; if ( t > 24.0 ) then t := t-24.0; coaltn := FABS(latitude-decl); Result := t; end; {* * Returns time on day no. nday of year when sun's coaltitude is coalt. * If no such time, then returns a large number. * time0 is approximate time of phenomenon *} function tempus(nday: integer; coalt, time0: double): double; { {* slong = true longitude *} {* ra = sun's right ascension, sindcl = sin(sun's declination) *} {* ha = sun's hour angle west *} {* locmt = local mean time of phenomenon *} var longh,days,anomaly,slong,sinslong,ra,sindcl,cosha,ha,locmt: double; t: double; begin longh := longitude*HPR; days := nday+(time0+longh)/24.0; anomaly := anom0+dmanom*days; slong := perigee+anomaly+c1*sin(anomaly)+c2*sin(anomaly*2); sinslong := sin(slong); ra := ArcTan2(cosobl*sinslong,cos(slong))*HPR; if ( ra < 0.0 ) then ra := ra+24.0; sindcl := sinobl*sinslong; cosha := (cos(coalt)-sindcl*sin(latitude))/(sqrt(1.0-sindcl*sindcl)*cos(latitude)); {* if cos(ha)>1, then time cannot be evaluated *} if ( FABS(cosha) > 1.0 ) then begin Result := 1.0e7; Exit; end; ha := ArcCos(cosha)*HPR; if ( time0 < 12.0 ) then ha := 24.0-ha; locmt := ha+ra-delsid*days-sidtm0; t := locmt+longh+timeZone; if ( t < 0.0 ) then t := t+24.0; if ( t> 24.0 ) then t := t-24.0; Result := t; end; {* * Returns the direction of Qibla in radians. Eastward from north is positive. *} function Qibla: double; var lat0, long0, dflong: double; begin {* lat0, long0 are Makkah's latitude and longitude in radians *} lat0 := 0.3739077; long0 := -0.69504828; dflong := longitude-long0; Result := ArcTan2(sin(dflong), cos(latitude)*tan(lat0)-sin(latitude)*cos(dflong)); {* cos(lat0)*tan(latitude)-sin(latitude)*cos(dflong)) ); *} end; {* * Print title material. *} procedure Header(day,month,year: integer); var direc: double; qibd, qibm: integer; sgnqib: string; begin { if ( latIsS = 1 ) then sgnlat := dir[3] else sgnlat := dir[2]; if ( longIsW = 1 ) then sgnlng := dir[1] else sgnlng := dir[0]; deg2dm(timeZone,zoneH,zoneM); if ( timeZone < 0 ) then sgnzon := '-' else sgnzon := '+'; } direc := Qibla*DPR; deg2dm(direc,qibd,qibm); if ( direc < 0 ) then sgnqib := dir[1] else sgnqib := dir[0]; myrPTime.QiblahD := qibd; myrPTime.QiblahM := qibm; myrPTime.QiblahO := sgnqib[1]; end; {* * Print times for range of days (in the year) first..last-1. * Days in month range from startdate to startdate+last-first *} procedure Display(first, last, startdate, cd, cm, header: integer); var hour, minute: integer; t: double; i, j, l: integer; s: string; begin i := startdate; for l := first to last-1 do begin myrPTime.Day := i; myrPTime.DoW := stdow; Inc(stdow); stdow := stdow mod 7; for j := 0 to 5 do begin t := tim[l,j]; if ( fajrByMaxInterval = 1 ) then begin if ( (j = 0) and ((t > 360.0) or (t < 0.0) or ((tim[l,1]-t) > fajrInterval)) ) then t := tim[l,1]-fajrInterval; if ( (j = 5) and ((t > 360.0) or (t < 0.0) or ((t-tim[l,4]) > fajrInterval)) ) then t := tim[l][4] + fajrInterval; end; if ( t > 360.0 ) then begin myrPTime.Mark := '*'; end else if ( t < 0.0 ) then begin myrPTime.Mark := ' '; end else begin {* time conversion to am and pm hours and rounded minutes *} hour := Trunc(t); minute := Trunc(60.0*(t-hour)+0.5); if ( minute >= 60 ) then begin minute := 0; hour := hour+1; end; // if ( hour>12 ) then // hour := hour-12; s := Format('%2d:%2d',[hour,minute]); if ( s[1] = ' ' ) then s[1] := '0'; if ( s[2] = ' ' ) then s[2] := '0'; if ( s[4] = ' ' ) then s[4] := '0'; if ( s[5] = ' ' ) then s[5] := '0'; case ( j ) of 0 : myrPTime.Fajr := s; 1 : myrPTime.Shuruq := s; 2 : myrPTime.Dhuhr := s; 3 : myrPTime.Asr := s; 4 : myrPTime.Maghrib := s; 5 : myrPTime.Isha := s; end; end; end; Inc(i); end; end; {* * Utility functions * *} function deg2rad(degree: double): double; begin Result := degree*RPD; end; function dm2deg(degree, minute: integer): double; begin Result := degree+(minute/60.0); end; procedure deg2dm(degree: double; var deg, min: integer); var dabs: double; begin dabs := FABS(degree); deg := Trunc(dabs); min := Trunc(60.0*(dabs-deg)+0.5); if ( min >= 60 ) then begin min := 0; deg := deg+1; end; end; function dms2deg(degree: longint; min: integer; sec: double): double; begin Result := degree+(min/60.0)+(sec/3600.0); end; function hms2h(hour, min: integer; sec: double): double; begin Result := hour+(min/60.0)+(sec/3600.0); end; end.