/*******************************************************************************/
/*                 Calculation of the Precession                               */
/*              (Astronomical algorithms 2nd p.134)                            */
/*******************************************************************************/

#define pai (3.1415926535897932/180.0)

void Precession (double year, double year0, double RA0,double DEC0,double V1,double V2,double *RA,double *DEC)
{
	double jd,jd0,t,T, Zeta,Zeta1,Zeta2,Zeta3,Theta,Theta1,Theta2,Theta3,Zed,Zed1,Zed2,Zed3;
	double CosTheta,SinTheta ,CosDec0,SinDec0,CosRaZeta ,SinRaZeta;

// year0: Epoch of the Star Catalogue
// year : year to calculate the position

	jd0=JD( year0,1,1,12,0,0);
	jd=JD( year,1,1,12,0,0);

	T=(jd0-2451545.0)/36525.0;
	t=(jd-jd0)/36525.0;

	Zeta1  = 2306.2181 + 1.39656*T - 0.000139*T*T;
	Zeta2  = 0.30188 - 0.000344*T;
	Zeta3  = 0.017998;
	Theta1 = 2004.3109 - 0.85330*T - 0.000217*T*T;
	Theta2 = - ( 0.42665 + 0.000217*T);
	Theta3 = - 0.041833;
	Zed1   = 2306.2181 + 1.39656*T - 0.000139*T*T;
	Zed2   = 1.09468 + 0.000066*T;
	Zed3   = 0.018203;

	Zeta  = (  Zeta1*t +  Zeta2*t*t + Zeta3*t*t*t )/3600.0;
	Zed   = (   Zed1*t +   Zed2*t*t +  Zed3*t*t*t )/3600.0;
	Theta = ( Theta1*t + Theta2*t*t +Theta3*t*t*t )/3600.0;

	Zeta = Zeta*pai;
	Zed  = Zed*pai;
	Theta= Theta*pai;
//	RA0  = (RA0 + V1*t*100/3600.0/cos(DEC0*pai))*pai; //BSC parameter
	RA0  = (RA0 + V1*15.0*t*100/3600.0)*pai;          //Sky parameter
	DEC0 = (DEC0+ V2*t*100/3600.0)*pai;

	CosTheta  = cos (Theta);
	SinTheta  = sin (Theta);
	CosDec0   = cos (DEC0);
	SinDec0   = sin (DEC0);
	CosRaZeta = cos(RA0+Zeta);
	SinRaZeta = sin(RA0+Zeta);

	*DEC   = asin (CosTheta*SinDec0 + SinTheta*CosDec0*CosRaZeta)/pai;
	*RA =atan2(CosDec0*SinRaZeta ,CosTheta*CosDec0*CosRaZeta-SinTheta*SinDec0)/pai;
	*RA =*RA + Zed/pai;
	if (*RA<0) *RA += 360.0;
}