# include "../include/astro.h"

main (argc, argv)
int argc;
char **argv;
{
  
  Coords coords;

  coords.crpix1 = 1024;
  coords.crpix2 = 512;
  coords.crval1 = 10.0;
  coords.crval2 = 89.0;

  coords.cdelt1 = 2.6;
  coords.cdelt2 = 2.6;

  coords.pc1_1 = 1.0;    coords.pc1_2 = 1.0;
  coords.pc2_1 = 1.0;    coords.pc2_2 = 1.0;


  while (fscanf (stdin, "%lf%lf", &ra, &dec)!=EOF) {
    RD_to_XY (&x, &y, ra, dec, &coords);
    XY_to_RD (&Ra, &Dec, x, y, &coords);
    fprintf (stderr, "%f %f --> %f %f --> %f %f\n", ra, dec, x, y, Ra, Dec);
  }

}
 
XY_to_RD (ra, dec, x, y, coords)
double *ra, *dec;
double  x, y;
Coords coords[];
{

  double L, M, X, Y, T;
  double R, sphi, cphi, stht, ctht;
  double alpha, delta, salp, calp, sdel, cdel, sdp, cdp;
  
  *ra  = 0;
  *dec = 0;

  X = coords[0].cdelt1*(x - coords[0].crpix1);
  Y = coords[0].cdelt2*(y - coords[0].crpix2);
  L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
  M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);

  fprintf (stderr, "XY LM: %f %f %f %f\n", X, Y, L, M);
  
  if (!strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
    R = hypot (L,M);
    if ((L == 0) && (M == 0)) {
      sphi = 0;
      cphi = 1;
    }
    else {
      sphi =  L / R;
      cphi = -M / R;
    }

    if (!strcmp(&coords[0].ctype[4], "-TAN") ) {
      if (R == 0) {
	stht = 1.0;
	ctht = 0.0;
      }
      else {
	T = DEG_RAD / R;
	stht =   T / sqrt ( 1.0 + T*T);
	ctht = 1.0 / sqrt ( 1.0 + T*T);
      }
    }
    if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
      ctht = RAD_DEG * R;
      stht = sqrt (1 - ctht*ctht);
    }

    sdp  = sin(RAD_DEG*coords[0].crval2);
    cdp  = cos(RAD_DEG*coords[0].crval2);
    
    sdel = stht*sdp - ctht*cphi*cdp;
    salp = ctht*sphi;
    calp = stht*cdp + ctht*cphi*sdp;
    alpha = atan2 (calp, salp);
    delta = asin (sdel);

    *ra  = DEG_RAD*alpha + coords[0].crval1;
    *dec = DEG_RAD*delta;
  }
  if (!strcmp(&coords[0].ctype[4], "-LIN") || !strcmp(&coords[0].ctype[0], "GENE")) {
    *ra  = L + coords[0].crval1;
    *dec = M + coords[0].crval2;
  }

}


RD_to_XY (x, y, ra, dec, coords)
double *x, *y;
double  ra, dec;
Coords coords[];
{

  double tmp_d;
  double X, Y, R, sphi, cphi, stht, ctht;
  double alpha, delta, salp, calp, sdel, cdel, sdp, cdp;

  *x = 0;
  *y = 0;
  /* 
  if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) {
    X = (ra  - coords[0].crval1);
    Y = (dec - coords[0].crval2);
  }
  */
  
  if (!strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
    sdp  = sin(RAD_DEG*coords[0].crval2);
    cdp  = cos(RAD_DEG*coords[0].crval2);
    salp = sin(RAD_DEG*(ra - coords[0].crval1));
    calp = cos(RAD_DEG*(ra - coords[0].crval1));
    sdel = sin(RAD_DEG*dec);
    cdel = cos(RAD_DEG*dec);

    stht = sdel*sdp + cdel*cdp*calp;  
    sphi = cdel*salp; 
    cphi = cdel*sdp*calp - sdel*cdp; 

    if (!strcmp(&coords[0].ctype[4], "-TAN") ) {
      X =  DEG_RAD * sphi / stht;
      Y = -DEG_RAD * cphi / stht;
    }
    if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
      X =  DEG_RAD * sphi;
      Y = -DEG_RAD * cphi;
    }
  }
  
  tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
  *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1;
  *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2;
 
}


/*** without sin/cos relationship simplifications **
  if (!strcmp(&ctype[4], "-TAN") || !strcmp(&ctype[4], "-SIN") || !strcmp(&ctype[0], "MM")) {
    X = (x - CRPIX1)*CDELT1;
    Y = (y - CRPIX2)*CDELT2;
  
    R = hypot (X,Y);
    if ((X == 0) && (tmp_Y == 0)) 
      phi = 0;
    else 
      phi = atan2 (X,-Y);  *** phi in radians **
    }
    sphi = sin(phi);
    cphi = cos(phi); 

    if (!strcmp(&ctype[4], "-TAN") ) {
      if (R == 0) 
	tht = 0.0;
      else 
	tht = atan (DEG_RAD / R); *** R in degrees, tht in radians **
    }
    if (!strcmp(&ctype[4], "-SIN") || !strcmp(&ctype[0], "MM")) {
      tht = acos (RAD_DEG * R);   *** R in degrees, tht in radians **
    }

    sdp  = sin(RAD_DEG*CRVAL2);
    cdp  = cos(RAD_DEG*CRVAL2);
    
    sdel = stht*sdp - ctht*cphi*cdp;
    salp = ctht*sphi;
    calp = stht*cdp + ctht*cphi*sdp;
    alpha = atan2 (salp, calp);
    delta = asin (sdel);
    
    *ra = DEG_RAD*alpha + CRVAL1;
    *dec = DEG_RAD*delta;
  }
*/
