Index: unk/Ohana/src/libdvo/src/coordops.update.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/coordops.update.c	(revision 14289)
+++ 	(revision )
@@ -1,662 +1,0 @@
-# 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;
-  
-  proj = GetProjection (coords[0].ctype);
-  mode = GetProjectionMode (proj);
-  if (proj == PROJ_NONE) return (FALSE);
-  if (proj == PROJ_MODE_NONE) 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 (mode == PROJ_MODE_CARTESIAN) {
-    *ra  = L + coords[0].crval1;
-    *dec = M + coords[0].crval2;
-
-    /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */
-    if (proj == PROJ_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 (mode == PROJ_MODE_ZENITHAL) {
-    R = hypot (L,M);
-    if ((L == 0) && (M == 0)) {
-      sphi = 0;
-      cphi = 1;
-    } else {
-      sphi =  L / R;
-      cphi = -M / R;
-    }
-
-    switch (proj) {
-      case PROJ_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);
-	}
-	break;
-      case PROJ_STG:
-	stht = (4 - RAD_DEG*R) / (4 + RAD_DEC*R);
-	ctht = sqrt (1 - stht*stht);
-	break;
-      case PROJ_SIN:
-	ctht = RAD_DEG * R;
-	stht = sqrt (1 - ctht*ctht);
-	break;
-      case PROJ_ZEA:
-      case PROJ_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);
-	break;
-    }
-    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 (mode == PROJ_MODE_PSEUDOCLY) {
-    switch (proj) {
-      case PROJ_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);
-      break;
-      
-      case PROJ_GLS:
-	/* L,M in degrees, alpha,delta in degrees */
-	alpha = L / cos (RAD_DEG * M);
-	delta = M;
-	break;
-      case PROJ_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);
-	break;
-    }
-    *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;
-
-  proj = GetProjection (coords[0].ctype);
-  mode = GetProjectionMode (proj);
-  if (proj == PROJ_NONE) return (FALSE);
-  if (proj == PROJ_MODE_NONE) return (FALSE);
-
-  /**** Locally Cartesian Projections ****/
-  if (mode == PROJ_MODE_CARTESIAN) {
-    if (proj == PROJ_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 (mode == PROJ_MODE_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;
-
-    switch (proj) {
-      case PROJ_TAN:
-      case PROJ_DIS:
-	L =  DEG_RAD * sphi / stht;
-	M = -DEG_RAD * cphi / stht;
-	break;
-      case PROJ_SIN:
-	L =  DEG_RAD * sphi;
-	M = -DEG_RAD * cphi;
-	break;
-      case PROJ_ZEA:
-      case PROJ_ZPL:
-	Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
-	L =  Rc * sphi;
-	M = -Rc * cphi;
-	status = TRUE;
-	break;
-    }
-  }
-
-  /**** Other Standard Projections ****/
-  if (mode == PROJ_MODE_PSEUDOCYL) {
-    switch (proj) {
-      case PROJ_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;
-	}	
-	break;
-      case PROJ_GLS:
-	phi = ra - coords[0].crval1;
-	theta = dec - coords[0].crval2;
-	L = phi * cos(RAD_DEG * theta);
-	M = theta;
-	break;
-      case PROJ_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);
-	break;
-    }
-  }
-
-  /* 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);
-
-      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'.
-
-*/
-
-  /* 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 */
-
-int GetProjection (char *ctype) {
-  if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA;
-  if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL;
-  if (!strcmp(&ctype[4], "-ARC")) return PROJ_ARC;
-  if (!strcmp(&ctype[4], "-STG")) return PROJ_STG;
-  if (!strcmp(&ctype[4], "-SIN")) return PROJ_SIN;
-  if (!strcmp(&ctype[0], "MM"))   return PROJ_SIN; // note ctype[0]
-  if (!strcmp(&ctype[4], "-TAN")) return PROJ_TAN;
-  if (!strcmp(&ctype[4], "-DIS")) return PROJ_DIS;
-  if (!strcmp(&ctype[4], "-LIN")) return PROJ_LIN;
-  if (!strcmp(&ctype[0], "GENE")) return PROJ_LIN; // note ctype[0]
-  if (!strcmp(&ctype[4], "-PLY")) return PROJ_PLY;
-  if (!strcmp(&ctype[4], "-WRP")) return PROJ_WRP;
-  if (!strcmp(&ctype[4], "-AIT")) return PROJ_AIT;
-  if (!strcmp(&ctype[4], "-GLS")) return PROJ_GLS;
-  if (!strcmp(&ctype[4], "-PAR")) return PROJ_PAR;
-  return PROJ_NONE;
-}
-  
-int SetProjection (char *ctype, int proj) {
-  switch (proj) {
-    case PROJ_ZEA: strcpy(&ctype[4], "-ZEA") return TRUE;
-    case PROJ_ZPL: strcpy(&ctype[4], "-ZPL") return TRUE;
-    case PROJ_ARC: strcpy(&ctype[4], "-ARC") return TRUE;
-    case PROJ_STG: strcpy(&ctype[4], "-STG") return TRUE;
-    case PROJ_SIN: strcpy(&ctype[4], "-SIN") return TRUE;
-    case PROJ_TAN: strcpy(&ctype[4], "-TAN") return TRUE;
-    case PROJ_DIS: strcpy(&ctype[4], "-DIS") return TRUE;
-    case PROJ_LIN: strcpy(&ctype[4], "-LIN") return TRUE;
-    case PROJ_PLY: strcpy(&ctype[4], "-PLY") return TRUE;
-    case PROJ_WRP: strcpy(&ctype[4], "-WRP") return TRUE;
-    case PROJ_AIT: strcpy(&ctype[4], "-AIT") return TRUE;
-    case PROJ_GLS: strcpy(&ctype[4], "-GLS") return TRUE;
-    case PROJ_PAR: strcpy(&ctype[4], "-PAR") return TRUE;
-  }
-  return FALSE;
-}  
-
-int GetProjectionMode (int proj) {
-  switch (proj) {
-    case PROJ_ZEA:
-    case PROJ_ZPL:
-    case PROJ_ARC:
-    case PROJ_STG:
-    case PROJ_SIN:
-    case PROJ_TAN:
-    case PROJ_DIS:
-      return PROJ_MODE_ZENITHAL;
-    case PROJ_LIN: 
-    case PROJ_PLY: 
-    case PROJ_WRP: 
-      return PROJ_MODE_CARTESIAN;
-    case PROJ_AIT:
-    case PROJ_GLS:
-    case PROJ_PAR:
-      return PROJ_MODE_PSEUDOCLY;
-  }
-  return PROJ_MODE_NONE;
-}
-