# include "astro.h"
    
CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output) {

  CoordTransform *transform;
  struct timeval now;
  struct tm *local;
  double T;

  ALLOCATE (transform, CoordTransform, 1);
  transform->isIdentity = FALSE;
  
  switch (input) {
    case COORD_CELESTIAL:
      switch (output) {
	case COORD_CELESTIAL:
	  transform->phi = 0.0;
	  transform->Xo  = 0.0;
	  transform->xo  = 0.0;
	  transform->isIdentity = TRUE;
	  break;
	case COORD_GALACTIC:
	  transform->phi = -62.60*RAD_DEG;
	  transform->Xo  = 282.25*RAD_DEG;
	  transform->xo  =  33.00;
	  break;
	case COORD_ECLIPTIC:
	  gettimeofday (&now, NULL);
	  local = localtime (&now.tv_sec);
	  T = local[0].tm_year / 100.0;
	  transform->phi = -1.0*RAD_DEG*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
	  transform->Xo  = 0.0;
	  transform->xo  = 0.0;
	  break;
	default:
	  abort();
      }
      break;
    case COORD_ECLIPTIC:
      switch (output) {
	case COORD_CELESTIAL:
	  gettimeofday (&now, (struct timezone *) NULL);
	  local = localtime (&now.tv_sec);
	  T = local[0].tm_year / 100.0;
	  transform->phi = RAD_DEG*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
	  transform->Xo = 0.0;
	  transform->xo = 0.0;
	  break;
	case COORD_GALACTIC:
	  return NULL;
	  break;
	case COORD_ECLIPTIC:
	  transform->phi = 0.0;
	  transform->Xo  = 0.0;
	  transform->xo  = 0.0;
	  transform->isIdentity = TRUE;
	  break;
	default:
	  abort();
      }
    case COORD_GALACTIC:
      switch (output) {
	case COORD_CELESTIAL:
	  transform->phi =  62.60*RAD_DEG;
	  transform->Xo  =  33.00*RAD_DEG;
	  transform->xo  = 282.25;
	  break;
	case COORD_GALACTIC:
	  transform->phi = 0.0;
	  transform->Xo  = 0.0;
	  transform->xo  = 0.0;
	  transform->isIdentity = TRUE;
	  break;
	case COORD_ECLIPTIC:
	  return NULL;
	default:
	  abort();
      }
    default:
      abort();
  }
 
  // pre-calculated constants:
  transform->sin_phi_cos_Xo = sin(transform->phi)*cos(transform->Xo);
  transform->sin_phi_sin_Xo = sin(transform->phi)*sin(transform->Xo);
  transform->cos_phi        = cos(transform->phi);
  
  transform->cos_phi_cos_Xo = cos(transform->phi)*cos(transform->Xo);
  transform->cos_phi_sin_Xo = cos(transform->phi)*sin(transform->Xo);
  transform->sin_phi        = sin(transform->phi);
  transform->cos_Xo 	    = cos(transform->Xo);
  transform->sin_Xo 	    = sin(transform->Xo);

  return transform;
}

// input and output coordinates are in degrees
int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform) {

  double sin_x, sin_y, cos_x, cos_y;


  if (transform == NULL) return (FALSE);

  if (transform->isIdentity) {
    *x = X;
    *y = Y;
    return (TRUE);
  }

  X *= RAD_DEG;
  Y *= RAD_DEG;

  // recast with constants extracted:
  sin_y = cos(Y)*sin(X)*transform->sin_phi_cos_Xo - cos(Y)*cos(X)*transform->sin_phi_sin_Xo + sin(Y)*transform->cos_phi;
  cos_y = sqrt (1 - sin_y*sin_y);

  sin_x = cos(Y)*sin(X)*transform->cos_phi_cos_Xo - cos(Y)*cos(X)*transform->cos_phi_sin_Xo - sin(Y)*transform->sin_phi;
  cos_x = cos(Y)*cos(X)*transform->cos_Xo + cos(Y)*sin(X)*transform->sin_Xo;
      
  // atan2 returns -pi : +pi
  *x = DEG_RAD * atan2 (sin_x, cos_x) + transform->xo;
  if ((*x) < 0.0) (*x) += 360;

  // should be in range -pi/2 : +pi/2
  *y = DEG_RAD * atan2 (sin_y, cos_y);

  return (TRUE);
}

// sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);

// sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
// cos_x = cos(Y)*cos(X - Xo) / cos_y;

// multiplying both sides by cos_y:
// sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi));
// cos_x = cos(Y)*cos(X - Xo);
      
// sin(a - b) = sin(a)*cos(b) - sin(b)*cos(a);
// cos(a - b) = sin(a)*sin(b) + cos(a)*cos(b);
