alqrah.cpp

VC+6.0 编写的关于计算穆斯林祷告的程序. 比较使用的应用小程序

/*
  CopyRight by Fayez Alhargan, 2001
  King Abdulaziz City for Science and Technology
  Computer and Electronics Research Institute
  Riyadh, Saudi Arabia
  [e-mail]
  Tel:[phone] Fax:[phone]

  This is a program that computes prayer times and sunrise.

  version: opn1.2
  last modified 20-5-2003

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser Library General Public License for more details.


*/


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
//#include <iostream.h>
#include <fstream.h>

#define pi 3.1415926535897932
#define DToR (pi / 180)
#define RToH (12 / pi)
#define EarthRadius 6378.14


      void GetRatior(int yg,int mg,int dg,double param[],double *IshRt,double *FajrRt);

/*------------------------- Sun ------------------------------------------*/
     double atanxy(double x,double y);
     void EclipToEquator(double lmdr,double betar,double &alph, double &dltr);
     double RoutinR2(double M,double e);
     double GCalendarToJD(int yg,int mg, double dg );
     double SunParamr(int yg,int mg, int dg,double ObsLon,double ObsLat, double TimeZone,
		double *Rise,double *Transit,double *Setting,double *RA,double *Decl,int *RiseSetFlags);


/*=====================================================================*/
/*             Computation for the Sun                                 */
/*=====================================================================*/
double atanxy(double x,double y)
{
     double argm;

     if(x==0)  argm=0.5*pi; else argm=atan(y/x);

     if(x>0 && y<0) argm=2.0*pi+argm;
     if(x<0) argm=pi+argm;

     return argm;
}
void EclipToEquator(double lmdr,double betar,double &alph, double &dltr)
{
 /*

   Convert Ecliptic to Equatorial Coordinate
   p.40 No.27, Peter Duffett-Smith book
   input: lmdr,betar  in radians
   output: alph,dltr in radians
 */

 double eps=23.441884;  // (in degrees) this changes with time
 double sdlt,epsr;
 double x,y,alpr;
 double rad=0.017453292;  // =pi/180.0

 epsr=eps*rad;  // convert to radians
 sdlt=sin(betar)*cos(epsr)+cos(betar)*sin(epsr)*sin(lmdr);
 dltr=asin(sdlt);
 y=sin(lmdr)*cos(epsr)-tan(betar)*sin(epsr);
 x=cos(lmdr);
 alph=atanxy(x,y);

}

double RoutinR2(double M,double e)
{
   /*
    Routine R2:
    Calculate the value of E
    p.91, Peter Duffett-Smith book
  */
  double dt=1,dE,Ec;
  Ec=M;
  while(fabs(dt)>1e-9)
  {
    dt=Ec-e*sin(Ec)-M;
    dE=dt/(1-e*cos(Ec));
    Ec=Ec-dE;
  }
 return Ec;
}

double SunParamr(int yg,int mg, int dg,double ObsLon,double ObsLat, double TimeZone,
		double *Rise,double *Transit,double *Setting,double *RA,double *Decl,int *RiseSetFlags)
{
  /*
    p.99 of the Peter Duffett-Smith book

  */

  double UT,ET,y,L,e,M,omg;
  double eps,T,JD,Ec;
  double tnv,v,tht;
  double K,angl,T1,T2,H,cH;
  *RiseSetFlags=0;

  JD=GCalendarToJD(yg,mg,dg);
  T=(JD+ TimeZone/24.0 - 2451545.0) / 36525.0;

  L=279.6966778+36000.76892*T+0.0003025*T*T;  // in degrees
  while(L>360) L=L-360;
  while(L<0) L=L+360;
  L=L*pi/180.0;  // radians

  M=358.47583+35999.04975*T-0.00015*T*T-0.0000033*T*T*T;
  while(M>360) M=M-360;
  while(M<0) M=M+360;
  M=M*pi/180.0;

  e=0.01675104-0.0000418*T-0.000000126*T*T;
  Ec=23.452294-0.0130125*T-0.00000164*T*T+0.000000503*T*T*T;
  Ec=Ec*pi/180.0;

  y=tan(0.5*Ec);
  y=y*y;
  ET=y*sin(2*L)-2*e*sin(M)+4*e*y*sin(M)*cos(2*L)-0.5*y*y*sin(4*L)-5*0.25*e*e*sin(2*M);
  UT=ET*180.0/(15.0*pi);   // from radians to hours

  Ec=RoutinR2(M,e);
  tnv=sqrt((1+e)/(1-e))*tan(0.5*Ec);
  v=2.0*atan(tnv);
  tht=L+v-M;
  EclipToEquator(tht,0,*RA,*Decl);

  K=12-UT-TimeZone+ObsLon*12.0/pi;  // (Noon)
  *Transit=K;
      /*  Sunrise and Sunset*/

   angl=(-0.833333)*DToR;  // Meeus p.98
   T1=(sin(angl)-sin(*Decl)*sin(ObsLat));
   T2=(cos(*Decl)*cos(ObsLat));  // p.38  Hour angle for the Sun
   cH=T1/T2;
   if(cH>1)  {*RiseSetFlags=16;cH=1;}  /*At this day and place the sun does not rise or set  */
   H=acos(cH);
   H=H*12.0/pi;
   *Rise=K-H; 	       // Sunrise
   *Setting=K+H; // SunSet

   return JD;
}

/*
  For international prayer times see Islamic Fiqah Council of the Muslim
  World League:  Saturday 12 Rajeb 1406H, concerning prayer times and fasting
  times for countries of high latitudes.
  This program is based on the above.
*/
/*****************************************************************************/
/* Name:    OmAlQrah                                                         */
/* Type:    Procedure                                                        */
/* Purpose: Compute prayer times and sunrise                                 */
/* Arguments:                                                                */
/*   yg,mg,dg : Date in Greg                                                 */
/*   param[0]: Safety time  in hours should be 0.016383h                     */
/*   longtud,latud: param[1],[2] : The place longtude and latitude in radians*/
/*   HeightdifW : param[3]: The place western herizon height difference in meters */
/*   HeightdifE : param[4]: The place eastern herizon height difference in meters */
/*   Zonh :param[5]: The place zone time dif. from GMT  West neg and East pos*/
/*          in decimal hours                                                 */
/*  fjrangl: param[6]: The angle (radian) used to compute                    */
/*            Fajer prayer time (OmAlqrah  -19 deg.)                         */
/*  ashangl: param[7]: The angle (radian) used to compute Isha  prayer time  */
/*          ashangl=0 then use  (OmAlqrah: ash=SunSet+1.5h)                  */
/*  asr  : param[8]: The Henfy (asr=2) Shafi (asr=1, Omalqrah asr=1)         */
/*  param[9]: latude (radian) that should be used for places above -+65.5    */
/*            should be 45deg as suggested by Rabita                         */
/*   param[10]: The Isha fixed time from Sunset                              */
/*  Output:                                                                  */
/*  lst[]: lst[n], 1:Fajer 2:Sunrise 3:Zohar 4:Aser  5:Magreb  6:Ishe        */
/*                 7:Fajer using exact Rabita method for places >48          */
/*                 8:Ash   using exact Rabita method for places >48          */
/*                 9: Eid Prayer Time                                        */
/*          for places above 48 lst[1] and lst[6] use a modified version of  */
/*          Rabita method that tries to eliminate the discontinuity          */
/*         all in 24 decimal hours                                           */
/*         returns flag:0 if there are problems, flag:1 no problems          */
/*****************************************************************************/
int OmAlQrahr(int yg,int mg,int dg, double param[], double lst[])
{
    int flag=1,flagrs,problm=0;
    double RA,Decl;
    double Rise,Transit,Setting;
    double SINd,COSd;
    double act,H,angl,K,cH;
    double X,MaxLat;
    double H0,Night,IshRt,FajrRt;
    double HightCorWest=0,HightCorEast=0;
    double IshFix,FajrFix;
  /*
    Main Local variables:
    RA= Sun's right ascension
    Decl= Sun's declination
    H= Hour Angle for the Sun
    K= Noon time
    angl= The Sun altitude for the required time
    flagrs: sunrise sunset flags
	    0:no problem
	    16: Sun always above horizon (at the ploes for some days in the year)
	    32: Sun always below horizon
  */

       /* Compute the Sun various Parameters */
  SunParamr(yg,mg,dg,-param[1],param[2],-param[5],
		&Rise,&Transit,&Setting,&RA,&Decl,&flagrs);

    /* Compute General Values */
  SINd=sin(Decl)*sin(param[2]);
  COSd=cos(Decl)*cos(param[2]);