# include <dvo.h>

static Coords *mosaic = NULL;

void RegisterMosaic (Coords *coords) {
  mosaic = coords;
}

int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) {

  int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl;
  char *type;
  double L, M, X, Y, T, Z, Z2;
  double R, sphi, cphi, stht, ctht;
  double alpha, delta, salp, calp, sdel, sdp, cdp;
  
  *ra  = 0;
  *dec = 0;
  stht = ctht = 1;
  type = &coords[0].ctype[4];
  
  /* PLY is equiv to LIN with higher order terms
     ZPL is equiv to ZEA with higher order terms
     DIS is equiv to TAN with higher order terms
     WRP is equiv to PLY, with implied mosaic */

  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP") || !strcmp(type, "-ZPL");
  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
  Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
  Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
  Zenithal   = !strcmp(type, "-ZEA") || !strcmp(type, "-ZPL") || Zenith1 || Zenith2;
  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);

  /** convert pixel coordinates to cartesian system **/
  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);

  /** extra polynomial terms **/
  if (coords[0].Npolyterms > 1) {
    L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
    M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
  }
  if (coords[0].Npolyterms > 2) {
    L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
    M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
  }

  /**** Locally Cartesian Projections ****/
  if (Cartesian) {
    *ra  = L + coords[0].crval1;
    *dec = M + coords[0].crval2;

    /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */
    if (!strcmp(type, "-WRP")) {
      if (mosaic == NULL) return (FALSE);
      XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic);
    }
    return (TRUE);
  }
  
  /**** Zenithal Projections ****/
  if (Zenithal) {
    R = hypot (L,M);
    if ((L == 0) && (M == 0)) {
      sphi = 0;
      cphi = 1;
    } else {
      sphi =  L / R;
      cphi = -M / R;
    }

    /* this is wrong : STG is not TAN - need to put in correct relationships.  but is a close approx */
    if (Zenith1) {
      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 (Zenith2) {
      ctht = RAD_DEG * R;
      stht = sqrt (1 - ctht*ctht);
    }
    if (!strcmp(type, "-ZEA") || !strcmp(type, "-ZPL")) {
      if (R > 2*DEG_RAD) {
	*ra = L;
	*dec = M;
	return (FALSE);
      }
      stht = 1 - 0.5*SQ(R*RAD_DEG);
      ctht = sqrt (1 - stht*stht);
    }

    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 (salp, calp);
    delta = asin (sdel);
    
    *ra  = DEG_RAD*alpha + coords[0].crval1;
    *dec = DEG_RAD*delta;

    /* rationalize ra range 0 - 360.0 */
    while (*ra < 0.0)   *ra += 360.0;
    while (*ra > 360.0) *ra -= 360.0;

    return (TRUE);
  }
  
  /**** Other Conventional Projections ****/
  if (PseudoCyl) {
    if (!strcmp(type, "-AIT")) {
      Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
      if (Z2 < 0) return (FALSE);
      Z = sqrt (Z2);
      alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0);
      delta = DEG_RAD * asin (RAD_DEG*M*Z);
      *ra  = alpha + coords[0].crval1;
      *dec = delta + coords[0].crval2;
    }
    if (!strcmp(type, "-GLS")) {
      /* L,M in degrees, alpha,delta in degrees */
      alpha = L / cos (RAD_DEG * M);
      delta = M;
      *ra  = alpha + coords[0].crval1;
      *dec = delta + coords[0].crval2;
    }
    if (!strcmp(type, "-PAR")) {
      /* L,M in degrees, alpha,delta in degrees */
      alpha = L / (1.0 - SQ(2.0*M/180));
      delta = 3 * DEG_RAD * asin (M/180.0);
      *ra  = alpha + coords[0].crval1;
      *dec = delta + coords[0].crval2;
    }

    /* rationalize ra range 0 - 360.0 */
    while (*ra < 0.0)   *ra += 360.0;
    while (*ra > 360.0) *ra -= 360.0;

    return (TRUE);
  }
  return (FALSE);
}

int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) {

  char *type;
  int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl;
  double phi, theta;
  double determ;
  double X, Y, L, M, Lo, Mo, dL, dM;
  double sphi, cphi, stht;
  double salp, calp, sdel, cdel, sdp, cdp;
  double P, A, Rc;

  status = TRUE;
  *x = 0;
  *y = 0;
  type = &coords[0].ctype[4];
  L = M = 0;

  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP") || !strcmp(type, "-ZPL");
  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
  Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
  Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
  Zenithal   = !strcmp(type, "-ZEA") || !strcmp(type, "-ZPL") || Zenith1 || Zenith2;
  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);

  /**** Locally Cartesian Projections ****/
  if (Cartesian) {
    if (!strcmp(type, "-WRP")) {
      if (mosaic == NULL) return (FALSE);
      RD_to_XY (&Lo, &Mo, ra, dec, mosaic);
      L = (Lo - coords[0].crval1);
      M = (Mo - coords[0].crval2);
    } else {
      L = (ra  - coords[0].crval1);
      M = (dec - coords[0].crval2);
    }
  }
  
  /**** Zenithal Projections ****/
  if (Zenithal)  {
    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;    /* sin(theta) */
    sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
    cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
    if (stht < 0) status = FALSE;

    if (Zenith1) {
      L =  DEG_RAD * sphi / stht;
      M = -DEG_RAD * cphi / stht;
    }
    if (Zenith2) {
      L =  DEG_RAD * sphi;
      M = -DEG_RAD * cphi;
    }
    if (!strcmp(type, "-ZEA") || !strcmp(type, "-ZPL")) {
      Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
      L =  Rc * sphi;
      M = -Rc * cphi;
      status = TRUE;
    }
  }

  /**** Other Standard Projections ****/
  if (PseudoCyl) {
    if (!strcmp(type, "-AIT")) {
      phi = RAD_DEG*(ra - coords[0].crval1);
      theta = RAD_DEG*(dec - coords[0].crval2);
      P = 1.0 + cos (theta) * cos (0.5*phi);
      if (P != 0.0) {
	A =  DEG_RAD * sqrt (2.0 / P);
	L =  2.0 * A * cos (theta) * sin (0.5*phi);
	M =  A * sin (theta);
      } else { 
	L =  0.0;
	M =  0.0;
      }	
    }
    if (!strcmp(type, "-GLS")) {
      phi = ra - coords[0].crval1;
      theta = dec - coords[0].crval2;
      L = phi * cos(RAD_DEG * theta);
      M = theta;
    }
    if (!strcmp(type, "-PAR")) {
      phi = ra - coords[0].crval1;
      theta = dec - coords[0].crval2;
      L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
      M = 180.0 * sin (RAD_DEG*theta/3.0);
    }
  }

  /* convert L,M to X,Y */
  determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
  X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M);
  Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L);

  /** extra polynomial terms **/
  if (coords[0].Npolyterms > 1) {
    for (i = 0; i < 10; i++) {
      Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
      Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
      if (coords[0].Npolyterms > 1) {
	Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
	Mo += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
      }
      if (coords[0].Npolyterms > 2) {
	Lo += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
	Mo += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
      }
      dL = (L - Lo);
      dM = (M - Mo);
      // fprintf (stderr, "%d: %f,%f : %f,%f : %f,%f : %f,%f\n", i, L, M, X, Y, Lo, Mo, dL, dM);

      X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM);
      Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dM);
    }
  }
  /* check for correct size (iterate?) */

  *x = X / coords[0].cdelt1 + coords[0].crpix1;
  *y = Y / coords[0].cdelt2 + coords[0].crpix2;

  return (status);
  
}

int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords) {

  int status;
  double tmpx, tmpy;

  status = RD_to_XY (&tmpx, &tmpy, ra, dec, coords);
  *x = tmpx;
  *y = tmpy;
  
  return (status);

}

int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords) {

  int status;
  double tmpr, tmpd;

  status = XY_to_RD (&tmpr, &tmpd, x, y, coords);
  *ra = tmpr;
  *dec = tmpd;
  
  return (status);

}

int GetCoords (Coords *coords, Header *header) {
  
  int status, itmp, Polynomial, Polyterm;
  double Lambda, rotate, scale;
  double equinox;
  char *ctype;
  
  rotate = 0.0;
  coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0;
  coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0;
  coords[0].pc1_1 = coords[0].pc2_2 = 1.0;
  coords[0].pc2_1 = coords[0].pc1_2 = 0.0;
  coords[0].Npolyterms = 1;
  strcpy (coords[0].ctype, "NONE");
  
  status = FALSE; 
  if (gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) {
    status  = gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1);
    status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1);
    status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2);  
    status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2);

    if (gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) {
      status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2);
      if (gfits_scan (header, "CROTA2", "%lf", 1, &rotate)) {
	Lambda = coords[0].cdelt2 / coords[0].cdelt1;
	coords[0].pc1_1 =  cos(rotate*RAD_DEG);
	coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
	coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
	coords[0].pc2_2 =  cos(rotate*RAD_DEG);
      }
      if (gfits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) {
	status &= gfits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2);
	status &= gfits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1);
	status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2);
      }

      /* set NPLYTERM based on header.  if NPLYTERM is missing, it should have a 
	 value of 0, unless the projection type is one of PLY, DIS, WRP, in which
	 case it should be set to 3 */
      ctype = &coords[0].ctype[4];
      Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
      Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp);

      coords[0].Npolyterms = 0;
      if (Polynomial && !Polyterm) coords[0].Npolyterms = 3;
      if (Polyterm) coords[0].Npolyterms = itmp;

      switch (coords[0].Npolyterms) {
	case 3:
	  status &= gfits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]);
	  status &= gfits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]);
	  status &= gfits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]);
	  status &= gfits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]);
	  status &= gfits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]);
	  status &= gfits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]);
	  status &= gfits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]);
	  status &= gfits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]);
	case 2:
	  status &= gfits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]);
	  status &= gfits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]);
	  status &= gfits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]);
	  status &= gfits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]);
	  status &= gfits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]);
	  status &= gfits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]);
	case 0:
	case 1:
	  break;
      }
    } else {
      if (gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) {
	status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2);
	status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1);
	status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2);
	/* renormalize */
	scale = hypot (coords[0].pc1_1, coords[0].pc1_2);
	coords[0].cdelt1 = coords[0].cdelt2 = scale;
	coords[0].pc1_1 /= scale;
	coords[0].pc1_2 /= scale;
	coords[0].pc2_1 /= scale;
	coords[0].pc2_2 /= scale;
      } else {
	status = FALSE;
      }
    }
  } else {
    /* some of my thesis data uses this simple linear model - convert on read? */
    if (gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) {
      status  = gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1);
      status &= gfits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2);
      status &= gfits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2);  
      status &= gfits_scan (header, "DEC_X", "%f", 1, &coords[0].pc2_1);
      status &= gfits_scan (header, "DEC_Y", "%f", 1, &coords[0].pc2_2);
      coords[0].crpix1 = coords[0].crpix2 = 0.0;
      coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
      strcpy (coords[0].ctype, "GENE");
    }
  }
  if (status) {
    if (!gfits_scan (header, "EQUINOX", "%lf", 1, &equinox)) {
      if (!gfits_scan (header, "EPOCH", "%lf", 1, &equinox)) {
	equinox = 2000.0;
      }
    }
    if (fabs (equinox - 2000.0) > 0.1) {
      coords_precess (&coords[0].crval1, &coords[0].crval2, equinox, 2000.0);
    } 
  }
  if (!status) {
    fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype);
    coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0;
    coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0;
    coords[0].pc1_1 = coords[0].pc2_2 = 1.0;
    coords[0].pc2_1 = coords[0].pc1_2 = 0.0;
    strcpy (coords[0].ctype, "NONE");
  }
  return (status);
}

int PutCoords (Coords *coords, Header *header) {
  
  int OldAIPS;
  char csys[16], ctype[16];
  double rotate, Lambda;

  /* modifications to the ctype? */
  OldAIPS = FALSE;
  gfits_modify (header, "CTYPE2",   "%s",  1, coords[0].ctype);
  if (!strcmp(coords[0].ctype, "MM")) {
    gfits_modify (header, "CTYPE1",   "%s",  1, "LL");
    OldAIPS = TRUE;
  } else {
    strcpy (csys, "NONE");
    if (!strncmp (coords[0].ctype, "DEC-", 4)) strcpy (csys, "RA--");
    if (!strncmp (coords[0].ctype, "GLAT", 4)) strcpy (csys, "GLON");
    if (!strncmp (coords[0].ctype, "ELAT", 4)) strcpy (csys, "ELON");
    if (!strncmp (coords[0].ctype, "HLAT", 4)) strcpy (csys, "HLON");
    if (!strncmp (coords[0].ctype, "SLAT", 4)) strcpy (csys, "SLON");
    if (!strcmp (csys, "NONE")) return (FALSE);
    sprintf (ctype, "%s-%s", csys, &coords[0].ctype[5]);
    gfits_modify (header, "CTYPE1",   "%s",  1, ctype);
  }    

  gfits_modify (header, "CDELT1",   "%le", 1, coords[0].cdelt1); 
  gfits_modify (header, "CDELT2",   "%le", 1, coords[0].cdelt2);
  gfits_modify (header, "CRVAL1",   "%lf", 1, coords[0].crval1);
  gfits_modify (header, "CRVAL2",   "%lf", 1, coords[0].crval2);  
  gfits_modify (header, "CRPIX1",   "%lf", 1, coords[0].crpix1);
  gfits_modify (header, "CRPIX2",   "%lf", 1, coords[0].crpix2);

  if (OldAIPS) {
    Lambda = coords[0].cdelt2 / coords[0].cdelt1;
    rotate = DEG_RAD*atan2 (coords[0].pc2_1*Lambda, coords[0].pc1_1);
    gfits_modify (header, "CROTA1", "%f", 1, rotate);
    gfits_modify (header, "CROTA2", "%f", 1, rotate);
    return (TRUE);
  } 

  gfits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
  gfits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
  gfits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
  gfits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
  gfits_modify (header, "NPLYTERM", "%d",  1, coords[0].Npolyterms);

  /* RA Terms */
  if (coords[0].Npolyterms > 1) {
    gfits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]);   /* polyterms[0]); */
    gfits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]);   /* polyterms[1]); */
    gfits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]);   /* polyterms[2]); */
  }
  if (coords[0].Npolyterms > 2) {
    gfits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]);   /* polyterms[3]); */
    gfits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]);   /* polyterms[4]); */
    gfits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]);   /* polyterms[5]); */
    gfits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]);   /* polyterms[6]); */
  }

  /* Dec Terms */
  if (coords[0].Npolyterms > 1) {
    gfits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]);   /* polyterms[7]); */
    gfits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]);   /* polyterms[8]); */
    gfits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]);   /* polyterms[9]); */
  }
  if (coords[0].Npolyterms > 2) {
    gfits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]);   /* polyterms[10]); */
    gfits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]);   /* polyterms[11]); */
    gfits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]);   /* polyterms[12]); */
    gfits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]);   /* polyterms[13]); */
  }
  return (TRUE);
}

void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch) {

  double T;
  double A, D, RA, DEC, zeta, z, theta;
  double SA, CA, SD, CD;
  
  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 = DEG_RAD*asin(SD);
  RA  = DEG_RAD*atan2(SA, CA) + z;
  
  if (RA < 0)
    RA += 360;
  
  *ra = RA;
  *dec = DEC; 
}

/* -PLY projection is an extrapolation of the -TAN projection. 
   In addition to the usual linear terms of CRPIXi, PC00i00j, there are 
   higher order polynomial terms (up to 3rd order):
   Axis 1 terms:
   PCA1X2Y0 = coords.polyterm[0][0] = x^2                                             
   PCA1X1Y1 = coords.polyterm[1][0] = xy                                          
   PCA1X0Y2 = coords.polyterm[2][0] = y^2                                             
   PCA1X3Y0 = coords.polyterm[3][0] = x^3                                             
   PCA1X2Y1 = coords.polyterm[4][0] = x^2 y                                             
   PCA1X1Y2 = coords.polyterm[5][0] = x y^2                                             
   PCA1X0Y3 = coords.polyterm[6][0] = y^2                                              
   Axis 2 terms:
   PCA2X2Y0 = coords.polyterm[0][1] = x^2                                               
   PCA2X1Y1 = coords.polyterm[1][1] = xy                                                
   PCA2X0Y2 = coords.polyterm[2][1] = y^2                                               
   PCA2X3Y0 = coords.polyterm[3][1] = x^3                                               
   PCA2X2Y1 = coords.polyterm[4][1] = x^2 y                                             
   PCA2X1Y2 = coords.polyterm[5][1] = x y^2                                             
   PCA2X0Y3 = coords.polyterm[6][1] = y^2                                               
*/

# if (0)

  /** convert pixel coordinates to cartesian system **/
  X = coords[0].cdelt1*(x - coords[0].crpix1);
  Y = coords[0].cdelt2*(y - coords[0].crpix2);
  if (Polynomi) {
    if (coords[0].Npolyterms > 2) {
      X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
      Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]);
    }
    if (coords[0].Npolyterms > 2) {
      X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]);
      Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]);
    }
  }

  L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
  M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
/** this code is the old method used for higher order terms.  they
    are essentially 6th order, with weird coupled terms.
    I don't think any real data used these terms, but they should 
    be re-calculated, I would think 
**/

# endif

/*

Projections and Transformations

The Coords structure is used to represent the standard FITS header representations of the WCS terms.
The FITS WCS encapsulates several concepts in one system: coordinate frame, projection, and
transformations.  The ctype variable defines both the frame (ie, RA/DEC, GLON/GLAT, etc) and the
projection (ie, AIT, SIN, TAN, etc).  The associated terms (PC00i00j) define transformations from
the projection system to another linear coordinate system.  I have extended the basic WCS
transformation terms to include higher-order polynomial terms.  The presence of higher-order terms
is indicated in the header by the NPLYTERM keyword, with a value greater than 1 (a value of 0 or 1
implies no higher-order coefficients).  The coefficients have keywords of the form PCAnXiYj where
'n' is the axis corresponding to the dependent variable, 'i' is the order of the X component and 'j'
is the order of the Y component.  A value of 2 indicates that the second-order coefficients are
defined (PCAnX2Y0, PCAnX1Y1, PCAnX0Y2); A value of 3 indicates that the third-order coefficients are
also defined (PCAnX3Y0, PCAnX2Y1, PCAnX1Y2, PCAnX0Y3).  Some headers in the past were generated
without the NPLYTERM keyword.  a value of PLY, DIS, or WRP for the projection without a
corresponding value for NPLYTERM implies that the value should be interpreted as '3'.

*/
