# include "astro.h"

void precess (Ra, Dec, in_epoch, out_epoch)
double *Ra, *Dec;
double in_epoch, out_epoch;
{

  double T, zeta, theta, z, A, D;
  double SD, CD, SA, CA;

  T = (out_epoch - in_epoch) / 100.0;
  
  zeta  = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T);
  theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T);
  z     =          0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T;

  A = *Ra;
  D = *Dec;

  SD =  cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D);
  CD = sqrt (1 - SD*SD);
  SA =  sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD;
  CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD;
  
  *Dec = asin(SD) / RAD_DEG;
  *Ra  = atan2(SA, CA) / RAD_DEG + z;
  
  if (*Ra < 0)
    *Ra += 360;
}  
   
