Index: /trunk/Ohana/src/libohana/src/coordops.c
===================================================================
--- /trunk/Ohana/src/libohana/src/coordops.c	(revision 3325)
+++ /trunk/Ohana/src/libohana/src/coordops.c	(revision 3326)
@@ -15,6 +15,7 @@
   type = &coords[0].ctype[4];
   
+  /* DIS and PLY are identical & apply projection; WRP is cartesian */
   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
-  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
   Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
@@ -58,5 +59,5 @@
     }
 
-    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
+    if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN")) {
       if (R == 0) {
 	stht = 1.0;
@@ -137,5 +138,5 @@
 
   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
-  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
   Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
@@ -162,5 +163,5 @@
     if (stht < 0) status = FALSE;
 
-    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
+    if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN")) {
       L =  DEG_RAD * sphi / stht;
       M = -DEG_RAD * cphi / stht;
@@ -269,7 +270,8 @@
 int GetCoords (Coords *coords, Header *header) {
   
-  int status, itmp;
-  double rotate, scale;
+  int status, itmp, Polynomial;
+  double Lambda, rotate, scale;
   double equinox;
+  char *ctype;
   
   rotate = 0.0;
@@ -290,7 +292,8 @@
       status &= fits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2);
       if (fits_scan (header, "CROTA2", "%f", 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);
-	coords[0].pc2_1 =  sin(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);
       }
@@ -300,5 +303,9 @@
 	status &= fits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2);
       }
-      if (!strcmp (coords[0].ctype, "DEC--PLY")) {
+
+      ctype = &coords[0].ctype[4];
+      Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
+
+      if (Polynomial) {
 	if (fits_scan (header, "NPLYTERM", "%d", 1, &itmp)) {
 	  coords[0].Npolyterms = itmp;
@@ -345,4 +352,5 @@
     }
   } else {
+    /* some of my thesis data uses this simple linear model - convert on read? */
     if (fits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) {
       status  = fits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1);
@@ -375,4 +383,77 @@
   }
   return (status);
+}
+
+int PutCoords (Coords *coords, Header *header) {
+  
+  int OldAIPS;
+  char csys[16], ctype[16];
+  double rotate, Lambda;
+
+  /* modifications to the ctype? */
+  OldAIPS = FALSE;
+  fits_modify (header, "CTYPE2",   "%s",  1, coords[0].ctype);
+  if (!strcmp(coords[0].ctype, "MM")) {
+    fits_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]);
+    fits_modify (header, "CTYPE1",   "%s",  1, ctype);
+  }    
+
+  fits_modify (header, "CDELT1",   "%le", 1, coords[0].cdelt1); 
+  fits_modify (header, "CDELT2",   "%le", 1, coords[0].cdelt2);
+  fits_modify (header, "CRVAL1",   "%lf", 1, coords[0].crval1);
+  fits_modify (header, "CRVAL2",   "%lf", 1, coords[0].crval2);  
+  fits_modify (header, "CRPIX1",   "%lf", 1, coords[0].crpix1);
+  fits_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);
+    fits_modify (header, "CROTA1", "%f", 1, rotate);
+    fits_modify (header, "CROTA2", "%f", 1, rotate);
+    return (TRUE);
+  } 
+
+  fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
+  fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
+  fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
+  fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
+  fits_modify (header, "NPLYTERM", "%d",  1, coords[0].Npolyterms);
+
+  /* RA Terms */
+  if (coords[0].Npolyterms > 1) {
+    fits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]);   /* polyterms[0]); */
+    fits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]);   /* polyterms[1]); */
+    fits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]);   /* polyterms[2]); */
+  }
+  if (coords[0].Npolyterms > 2) {
+    fits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]);   /* polyterms[3]); */
+    fits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]);   /* polyterms[4]); */
+    fits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]);   /* polyterms[5]); */
+    fits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]);   /* polyterms[6]); */
+  }
+
+  /* Dec Terms */
+  if (coords[0].Npolyterms > 1) {
+    fits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]);   /* polyterms[7]); */
+    fits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]);   /* polyterms[8]); */
+    fits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]);   /* polyterms[9]); */
+  }
+  if (coords[0].Npolyterms > 2) {
+    fits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]);   /* polyterms[10]); */
+    fits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]);   /* polyterms[11]); */
+    fits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]);   /* polyterms[12]); */
+    fits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]);   /* polyterms[13]); */
+  }
+  return (TRUE);
 }