Index: /trunk/Ohana/src/gastro2/Makefile
===================================================================
--- /trunk/Ohana/src/gastro2/Makefile	(revision 8240)
+++ /trunk/Ohana/src/gastro2/Makefile	(revision 8241)
@@ -23,4 +23,10 @@
 $(SRC)/gregions2.$(ARCH).o
 
+COORDTEST = \
+$(SRC)/coordtest.$(ARCH).o \
+$(SRC)/gaussj.$(ARCH).o \
+$(SRC)/gpairs.$(ARCH).o \
+$(SRC)/polyfit.$(ARCH).o
+
 GASTRO = \
 $(SRC)/getptolemy.$(ARCH).o \
@@ -29,4 +35,6 @@
 $(SRC)/gheader2.$(ARCH).o \
 $(SRC)/gfit2.$(ARCH).o \
+$(SRC)/gpairs.$(ARCH).o \
+$(SRC)/polyfit.$(ARCH).o \
 $(SRC)/plotstuff.$(ARCH).o \
 $(SRC)/rotate2.$(ARCH).o \
@@ -62,4 +70,7 @@
 	$(CC) $(2MASS) -o $(BIN)/extr2mass.linux $(CCFLAGS)
 
+coordtest: $(COORDTEST)
+	$(CC) $(COORDTEST) -o $(BIN)/coordtest.linux $(CCFLAGS)
+
 # dependancy rules for binary code ##########################
 $(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH)
Index: /trunk/Ohana/src/gastro2/include/gastro2.h
===================================================================
--- /trunk/Ohana/src/gastro2/include/gastro2.h	(revision 8240)
+++ /trunk/Ohana/src/gastro2/include/gastro2.h	(revision 8241)
@@ -3,29 +3,10 @@
 
 typedef struct {
-  double R;
-  double D;
-  double X;
-  double Y;
-  double M;
-  double dM;
+  double R, D;
+  double P, Q;
+  double X, Y;
+  double M, dM;
   int type;
 } StarData;
-
-# if (0)
-typedef struct {
-  float X;
-  float Y;
-  float M;
-  float dM;
-  float Mgal;
-  float Map;
-  float sky;
-  float fx;
-  float fy;
-  float df;
-  char  dophot;
-  char  dummy[3];
-} Stars;
-# endif
 
 typedef struct {
@@ -129,5 +110,7 @@
 double MAGLIM_MIN;
 double MAGLIM_MAX;
-
+int NGRID_PIX;
+
+int    ASCA;
 int    FORCE;
 double F_RA;
@@ -173,5 +156,5 @@
 void 	  fit_lum_bin (double *x, double *y, int N, double *C0, double *C1);
 void 	  fit_norm (); 
-double    fit_scat (StarData *st, StarData *sr);
+double    fit_scat (StarData *st, StarData *sr, Coords *coords);
 int 	  gaussj (double **a, int n, double **b, int m);
 void 	  gcenter (CmpCatalog *Target, RefCatalog *Ref);
@@ -209,5 +192,5 @@
 void 	  plot_init_gridplot ();
 void 	  plot_lumfunc (CmpCatalog *Target, RefCatalog *Ref);
-void 	  plot_resid (StarData *st, StarData *sr);
+void 	  plot_resid (StarData *st, StarData *sr, Coords *coords);
 void 	  plot_resid_init (int version, double xmax);
 void 	  plot_resid_plot (int version, float *xvect, float *yvect, int Nvect);
Index: /trunk/Ohana/src/gastro2/src/ConfigInit.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/ConfigInit.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/ConfigInit.c	(revision 8241)
@@ -27,4 +27,5 @@
   ScanConfig (config, "ASEC_PIX",          "%lf", 0, &ASEC_PIX);
   ScanConfig (config, "NFIELD",            "%lf", 0, &NFIELD);
+  ScanConfig (config, "NGRID_PIX",         "%d",  0, &NGRID_PIX);
   ScanConfig (config, "NPOLYTERMS",        "%d",  0, &NPOLYTERMS);
   ScanConfig (config, "ROT_ZERO",          "%lf", 0, &ROT_ZERO);
@@ -58,4 +59,5 @@
       exit (1);
   }
+  if (NGRID_PIX == 0) NGRID_PIX = 50.0;
 
   if (strcasecmp (ROUGH_ASTROMETRY, "header") && 
Index: /trunk/Ohana/src/gastro2/src/coordtest.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/coordtest.c	(revision 8241)
+++ /trunk/Ohana/src/gastro2/src/coordtest.c	(revision 8241)
@@ -0,0 +1,143 @@
+# include "gastro2.h"
+
+# define A00 +500.00
+// # define A00 +0.00
+# define A10 +4.00
+# define A01 +1.00
+# define A20 -0.001
+# define A11 +0.000
+# define A02 -0.001
+# define A30 +0.002
+# define A21 +0.001
+# define A12 -0.001
+# define A03 -0.002
+
+# define B00 400.00
+// # define B00 +0.00
+# define B10  -1.00
+# define B01  4.00
+# define B20 -0.001
+# define B11  0.000
+# define B02 -0.001
+# define B30 -0.002
+# define B21 -0.001
+# define B12 +0.001
+# define B03 +0.002
+
+int main (int argc, char **argv) {
+
+  /* generate a set of fake data and send to fitter */
+  int i, N, order;
+  double x, y, z;
+  double L[50000], M[50000], X[50000], Y[50000];
+  Coords coords;
+  FILE *f;
+
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: coordtest (order)\n");
+    exit (2);
+  }
+
+  order = atoi (argv[1]);
+  fit_init (order);
+  N = 0;
+
+  switch (order) {
+    case 1:
+      for (x = -10.0; x < 10.1; x += 0.1) {
+	for (y = -10.0; y < 10.1; y += 0.1) {
+      
+	  L[N] = A00 + A10*x + A01*y;
+	  M[N] = B00 + B10*x + B01*y;
+	  X[N] = x;
+	  Y[N] = y;
+	  fit_add (x, y, L[N], M[N], 1.0);
+	  N++;
+	}
+      }
+      break;
+    case 2:
+      for (x = -10.0; x < 10.1; x += 0.1) {
+	for (y = -10.0; y < 10.1; y += 0.1) {
+      
+	  L[N] = A00 + A10*x + A01*y + A20*x*x + A11*x*y + A02*y*y;
+	  M[N] = B00 + B10*x + B01*y + B20*x*x + B11*x*y + B02*y*y;
+	  X[N] = x;
+	  Y[N] = y;
+	  fit_add (x, y, L[N], M[N], 1.0);
+	  N++;
+	}
+      }
+      break;
+    case 3:
+      for (x = -10.0; x < 10.1; x += 0.1) {
+	for (y = -10.0; y < 10.1; y += 0.1) {
+      
+	  L[N] = A00 + A10*x + A01*y + A20*x*x + A11*x*y + A02*y*y + A30*x*x*x + A21*x*x*y + A12*x*y*y + A03*y*y*y;
+	  M[N] = B00 + B10*x + B01*y + B20*x*x + B11*x*y + B02*y*y + B30*x*x*x + B21*x*x*y + B12*x*y*y + B03*y*y*y;
+	  X[N] = x;
+	  Y[N] = y;
+	  fit_add (x, y, L[N], M[N], 1.0);
+	  N++;
+	}
+      }
+      break;
+  }
+
+  fit_eval ();
+
+  strcpy (coords.ctype, "RA---PLY");
+  coords.crval1 = 0.0;
+  coords.crval2 = 0.0;
+
+  fit_adjust (&coords);
+
+  fprintf (stderr, "CTYPE: %s\n", coords.ctype);
+  fprintf (stderr, "CRVAL: %f %f\n", coords.crval1, coords.crval2);
+  fprintf (stderr, "CRPIX: %f %f\n", coords.crpix1, coords.crpix2);
+  fprintf (stderr, "CDELT: %f %f\n", coords.cdelt1, coords.cdelt2);
+  fprintf (stderr, "PC1_j: %f %f\n", coords.pc1_1,  coords.pc1_2);
+  fprintf (stderr, "PC2_j: %f %f\n", coords.pc2_1,  coords.pc2_2);
+
+  fprintf (stderr, "Npolyterms: %d\n", coords.Npolyterms);
+  for (i = 0; i < 7; i++) {
+    fprintf (stderr, "%f %f\n", coords.polyterms[i][0], coords.polyterms[i][1]);
+  }
+
+  f = fopen ("test.dat", "w");
+  {
+    double Lo, Mo, dL, dL2, dM, dM2;
+    double Xo, Yo, dX, dX2, dY, dY2;
+
+    dL = dL2 = dM = dM2 = 0.0;
+    dX = dX2 = dY = dY2 = 0.0;
+    for (i = 0; i < N; i++) {
+      XY_to_RD (&Lo, &Mo, X[i], Y[i], &coords);
+      dL  += (L[i] - Lo);
+      dM  += (M[i] - Mo);
+      dL2 += SQ(L[i] - Lo);
+      dM2 += SQ(M[i] - Mo);
+
+      RD_to_XY (&Xo, &Yo, L[i], M[i], &coords);
+      dX  += (X[i] - Xo);
+      dY  += (Y[i] - Yo);
+      dX2 += SQ(X[i] - Xo);
+      dY2 += SQ(Y[i] - Yo);
+
+      // fit_apply (&Lx, &Mx, X[i], Y[i]);
+      // fprintf (stderr, "%f,%f -> %f,%f | %f,%f : %f, %f\n", X[i], Y[i], Lx, Mx, L[i], M[i], Lo, Mo);
+
+      fprintf (f, "%f %f : %f %f :: %f %f : %f %f\n",
+	       X[i], Y[i], L[i], M[i], Lo, Mo, Xo, Yo);
+    }
+    fclose (f);
+
+    fprintf (stderr, "dL: %f\n", sqrt(fabs(dL2/N - SQ(dL/N))));
+    fprintf (stderr, "dM: %f\n", sqrt(fabs(dM2/N - SQ(dM/N))));
+    fprintf (stderr, "dX: %f\n", sqrt(fabs(dX2/N - SQ(dX/N))));
+    fprintf (stderr, "dY: %f\n", sqrt(fabs(dY2/N - SQ(dY/N))));
+  }
+
+  exit (0);
+}
+
Index: /trunk/Ohana/src/gastro2/src/gargs.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/gargs.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/gargs.c	(revision 8241)
@@ -122,4 +122,12 @@
   }
 
+  /** XXX temporary trick to deal with the very wide-field ASCA images 
+      this alters the definition of the reference field boundaries */
+  ASCA = FALSE;
+  if ((N = get_argument (*argc, argv, "-asca"))) {
+    ASCA = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
   CATDUMP = FALSE;
   if ((N = get_argument (*argc, argv, "-dump"))) {
Index: /trunk/Ohana/src/gastro2/src/gastro2.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/gastro2.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/gastro2.c	(revision 8241)
@@ -30,5 +30,5 @@
   }
 
-  for (i = 0; i < 3; i++) {
+  for (i = 0; i < 5; i++) {
     gfit (&Target, &Ref, MIN (MAX (1, NPOLYTERMS), 3));
     fprintf (stderr, "precision: %f\n", Target.answer.dR / sqrt(Target.answer.N));
Index: /trunk/Ohana/src/gastro2/src/gfit2.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/gfit2.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/gfit2.c	(revision 8241)
@@ -1,10 +1,8 @@
 # include "gastro2.h"
-
-static int NPAIR, Npair;
-static int *idx1, *idx2;
 
 void gfit (CmpCatalog *Target, RefCatalog *Ref, int order) {
 
   int i, j, j0;
+  int Npair, *idx1, *idx2;
   double Radius, Radius2;
   double dX, dY, dR;
@@ -12,4 +10,5 @@
   StarData *st, *sr;
 
+  /* XXX why is this hardwired here? */
   NFIELD = 0.1;
   gproject (Target, Ref, &Subset);
@@ -25,5 +24,5 @@
   sort_stars_X (Subset.stars, Subset.N);
 
-  Radius = MAX (2.0 * Target[0].answer.dR, 0.5);
+  Radius = MAX (2.5 * Target[0].answer.dR, 0.5);
   Radius2 = Radius*Radius;
 
@@ -63,384 +62,21 @@
   }
   
-# if (0)
-  pair_init ();
-  for (i = 0; i < Target[0].N; i++) {
-    pair_add (i, i);
-  }
-# endif
-
-  fit_init (order);
-  fit_norm ();
-
-  plot_resid_init (0, (double) Target[0].header.Naxis[0]);
-  plot_resid_init (1, (double) Target[0].header.Naxis[1]);
-  if (PLOTSTUFF) plot_resid (st, sr);
-
+  Npair = pair_lists (&idx1, &idx2);
   /* find fit for matched pairs */
   fit_init (order);
   for (i = 0; i < Npair; i++) {
-    fit_add (st[idx1[i]].X, st[idx1[i]].Y, sr[idx2[i]].X, sr[idx2[i]].Y, 1.0);
+    fit_add (st[idx1[i]].X, st[idx1[i]].Y, sr[idx2[i]].P, sr[idx2[i]].Q, 1.0);
   }
   fit_eval ();
 
-  Target[0].answer.dR = fit_scat (st, sr);
+  /* XXX this is weak: the fit_scat call requires the coords from the fit_adjust call */
   Target[0].answer.N  = fit_adjust (&Target[0].coords);
+  Target[0].answer.dR = fit_scat (st, sr, &Target[0].coords);
 
+  if (PLOTSTUFF) fprintf (stderr, "ploting resid (2)\n");
   plot_resid_init (0, (double) Target[0].header.Naxis[0]);
   plot_resid_init (1, (double) Target[0].header.Naxis[1]);
-  if (PLOTSTUFF) plot_resid (st, sr);
+  if (PLOTSTUFF) plot_resid (st, sr, &Target[0].coords);
   free (idx1);
   free (idx2);
 }
-
-static int NTERM, NPOWR, NPARS, NORDER, Npts;
-static double **sum, **xsum, **ysum;
-static double **matrix, **vector;
-
-void fit_init (int order) {
-
-  int i;
-
-  NORDER = order;
-  NPOWR = NORDER + 1;
-  NTERM = 2*NORDER + 1;
-  NPARS = (NORDER + 1)*(NORDER + 2) / 2;
-  Npts  = 0;
-
-  /* allocate arrays for fit solution */
-  ALLOCATE (sum, double *, NTERM);
-  ALLOCATE (xsum, double *, NTERM);
-  ALLOCATE (ysum, double *, NTERM);
-  for (i = 0; i < NTERM; i++) {
-    ALLOCATE (sum[i], double, NTERM);
-    bzero (sum[i], NTERM*sizeof(double));
-    ALLOCATE (xsum[i], double, NTERM);
-    bzero (xsum[i], NTERM*sizeof(double));
-    ALLOCATE (ysum[i], double, NTERM);
-    bzero (ysum[i], NTERM*sizeof(double));
-  }
-  ALLOCATE (matrix, double *, NPARS);
-  ALLOCATE (vector, double *, NPARS);
-  for (i = 0; i < NPARS; i++) {
-    ALLOCATE (matrix[i], double, NPARS);
-    ALLOCATE (vector[i], double, 2);
-    bzero (vector[i], 2*sizeof(double));
-    bzero (matrix[i], NPARS*sizeof(double));
-  }
-
-}
-
-void fit_add (double x1, double y1, double x2, double y2, double wt) {
-
-  int n, m;
-  double xterm, yterm, term;
-
-  xterm = 1;
-  for (n = 0; n < NTERM; n++) {
-    yterm = 1;
-    for (m = 0; m < NTERM; m++) {
-      term = xterm*yterm;
-      if (n+m < NTERM) {
-	sum[n][m] += term;
-      }
-      if (n+m < NPOWR) {
-	xsum[n][m] += x2*term;
-	ysum[n][m] += y2*term;
-      }
-      yterm *= y1;
-    }
-    xterm *= x1;
-  }
-  Npts ++;
-
-}
-
-void fit_eval () {
-
-  int i, j, n, m, M, N;
-  double max;
-
-  fprintf (stderr, "npts: %d\n", Npts);
-
-  i = 0;
-  for (m = 0; m < NPOWR; m++) {
-    for (n = 0; n < NPOWR - m; n++, i++) {
-      vector[i][0] = xsum[n][m];
-      vector[i][1] = ysum[n][m];
-    }	
-  }
-  j = 0;
-  for (M = 0; M < NPOWR; M++) {
-    for (N = 0; N < NPOWR - M; N++, j++) {
-      i = 0;
-      for (m = 0; m < NPOWR; m++) {
-	for (n = 0; n < NPOWR - m; n++, i++) {
-	  matrix[i][j] = sum[n+N][m+M];
-	}	
-      }
-    }
-  }       
-  max = 0.0;
-  for (i = 0; i < NPARS; i++) {
-    for (j = 0; j < NPARS; j++) {
-      max = MAX (max, fabs(matrix[i][j]));
-    }
-    max = MAX (max, fabs(vector[i][0]));
-    max = MAX (max, fabs(vector[i][1]));
-  }
-  for (i = 0; i < NPARS; i++) {
-    for (j = 0; j < NPARS; j++) {
-      matrix[i][j] /= max;
-    }
-    vector[i][0] /= max;
-    vector[i][1] /= max;
-  }
-# if (0)
-  for (i = 0; i < NPARS; i++) {
-    for (j = 0; j < NPARS; j++) {
-      fprintf (stderr, "%10.4e  ", matrix[i][j]);
-    }
-    fprintf (stderr, "  %10.4e  ", vector[i][0]);
-    fprintf (stderr, "  %10.4e  \n", vector[i][1]);
-  }
-# endif
-  gaussj (matrix, NPARS, vector, 2); 
-# if (0)
-  fprintf (stderr, "\n\n");
-  for (i = 0; i < NPARS; i++) {
-    for (j = 0; j < NPARS; j++) {
-      fprintf (stderr, "%11.4e  ", matrix[i][j]);
-    }
-    fprintf (stderr, "  %11.4e  ", vector[i][0]);
-    fprintf (stderr, "  %11.4e  \n", vector[i][1]);
-  }
-# endif
-  i = 0;
-  for (m = 0; m < NPOWR; m++) {
-    for (n = 0; n < NPOWR - m; n++, i++) {
-      xsum[n][m] = vector[i][0];
-      ysum[n][m] = vector[i][1];
-    }	
-  }
-  i = 0;
-  for (m = 0; m < NPOWR; m++) {
-    for (n = 0; n < NPOWR - m; n++, i++) {
-      fprintf (stderr, "RA x^%dy^%d: %10.4g    DEC x^%dy^%d: %10.4g \n", 
-	       n, m, vector[i][0], n, m, vector[i][1]);
-    }	
-  }
-}
-
-void fit_norm () { 
-
-  xsum[0][0] = 0;
-  xsum[1][0] = 1;
-  xsum[0][1] = 0;
-
-  ysum[0][0] = 0;
-  ysum[1][0] = 0;
-  ysum[0][1] = 1;
-}
-
-void fit_apply (double *x, double *y, double X, double Y) {
-
-  int m, n;
-  double xterm, yterm;
-  double Xo, Yo;
-
-  Xo = Yo = 0;
-  yterm = 1;
-  for (m = 0; m < NPOWR; m++) { 
-    xterm = 1;
-    for (n = 0; n < NPOWR - m; n++) {
-      Xo += xterm*yterm*xsum[n][m];
-      Yo += xterm*yterm*ysum[n][m];
-      xterm *= X;
-    }	
-    yterm *= Y;
-  }
-  
-  *x = Xo;
-  *y = Yo;
-}
-
-
-void pair_init () {
-
-  Npair = 0;
-  NPAIR = 100;
-  ALLOCATE (idx1, int, NPAIR);
-  ALLOCATE (idx2, int, NPAIR);
-
-}
-
-void pair_add (int i1, int i2) {
-
-  idx1[Npair] = i1;
-  idx2[Npair] = i2;
-
-  Npair ++;
-
-  if (Npair == NPAIR) {
-    NPAIR += 100;
-    REALLOCATE (idx1, int, NPAIR);
-    REALLOCATE (idx2, int, NPAIR);
-  }
-
-}
-
-double fit_scat (StarData *st, StarData *sr) {
-
-  int i;
-  double x, y, dx, dy, dX, dY, dX2, dY2, dR;
-  
-  dX = dY = dX2 = dY2 = 0;
-  for (i = 0; i < Npair; i++) {
-
-    fit_apply (&x, &y, sr[idx2[i]].X, sr[idx2[i]].Y);
-    
-    dx = x - st[idx1[i]].X;
-    dy = y - st[idx1[i]].Y;
-    
-    dX += dx;
-    dY += dy;
-    dX2 += dx*dx;
-    dY2 += dy*dy;
-    
-  }
-
-  dX = dX / Npair;
-  dY = dY / Npair;
-  fprintf (stderr, "scatter: %f, %f\n", sqrt(dX2/Npair - dX*dX), sqrt(dY2/Npair - dY*dY));
-  fprintf (stderr, "precise: %f, %f\n", sqrt((dX2/Npair - dX*dX) / Npair), sqrt((dY2/Npair - dY*dY)/Npair));
-
-  dR = 0.5 * sqrt(fabs(dX2/Npair - dX*dX)) + 0.5 * sqrt (fabs(dY2/Npair - dY*dY));
-  return (dR);
-}
-
-/* convert new terms to adjustments in coords and to polyterms */
-int fit_adjust (Coords *coords) {
-
-  int i, j, N;
-  double S1, S2, p11, p12, p21, p22;
-  double a0, a1, a2, b0, b1, b2, det;
-  double X, Y;
-  int Np, Nv;
-    
-  S1 = coords[0].cdelt1;
-  S2 = coords[0].cdelt2;
-  p11 = coords[0].pc1_1;    p12 = coords[0].pc1_2;
-  p21 = coords[0].pc2_1;    p22 = coords[0].pc2_2;
-    
-  /* get the correct vector entries for the linear terms */
-  N = mk_vector (0, 0, NORDER);
-  a0 = vector[N][0];  b0 = vector[N][1];
-  N = mk_vector (1, 0, NORDER);
-  a1 = vector[N][0];  b1 = vector[N][1];
-  N = mk_vector (0, 1, NORDER);
-  a2 = vector[N][0];  b2 = vector[N][1];
-
-  det = 1.0 / (a1*b2 - a2*b1);
-
-  coords[0].pc1_1 = p11*a1 + p12*b1*(S2/S1);
-  coords[0].pc2_1 = p21*a1 + p22*b1*(S2/S1);
-    
-  coords[0].pc1_2 = p12*b2 + p11*a2*(S1/S2);
-  coords[0].pc2_2 = p22*b2 + p21*a2*(S1/S2);
-    
-  X = (coords[0].crpix1 - a0);
-  Y = (coords[0].crpix2 - b0);
-  coords[0].crpix1 = det*(X*b2 - Y*a2);
-  coords[0].crpix2 = det*(Y*a1 - X*b1);
-
-  coords[0].Npolyterms = NORDER;
-  if (NORDER > 1) strcpy (coords[0].ctype, "DEC--PLY");
-
-  /* generate higher order terms from vector */
-  for (i = 0; i < NORDER + 1; i++) {
-    for (j = 0; j < (NORDER - i + 1); j++) {
-      if (i + j < 2) continue;
-      Np = mk_polyterm (i, j, NORDER);
-      Nv = mk_vector (i, j, NORDER);
-      coords[0].polyterms[Np][0] = det*(vector[Nv][0]*b2  - vector[Nv][1]*a2);  /* x2 y0 */
-      coords[0].polyterms[Np][1] = det*(vector[Nv][1]*a1  - vector[Nv][0]*b1);  /* x2 y0 */
-    }
-  }
-  while (coords[0].crval1 < 0) coords[0].crval1 += 360.0;
-  while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
-
-  /* test for valid solution */
-  return (Npts);
-
-}
-
-int mk_polyterm (int n, int m, int norder) {
-  
-  int i, nt, N;
-  
-  N = 0;
-  nt = n + m;
-  for (i = 2; i < nt; i++) {
-    N += i + 1;
-  }
-  N += m;
-  return (N);
-}
-
-int mk_vector (int n, int m, int norder) {
-  
-  int i, N;
-  
-  N = 0;
-  for (i = 0; i < m; i++) {
-    N += (norder - i + 1);
-  }
-  N += n;
-  return (N);
-}
-
-void plot_resid (StarData *st, StarData *sr) {
-
-  int i;
-  double x, y, dx, dy;
-  float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2;
-  int Nvect;
-
-  Nvect = Npair;
-  ALLOCATE (xvect0, float, Nvect);
-  ALLOCATE (yvect0, float, Nvect);
-  ALLOCATE (xvect1, float, Nvect);
-  ALLOCATE (yvect1, float, Nvect);
-
-  ALLOCATE (xvect2, float, 2*Nvect);
-  ALLOCATE (yvect2, float, 2*Nvect);
-  
-  for (i = 0; i < Npair; i++) {
-    fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y);
-    
-    dx = x - sr[idx2[i]].X;
-    dy = y - sr[idx2[i]].Y;
-    
-    xvect0[i] = sr[idx2[i]].X;
-    xvect1[i] = sr[idx2[i]].Y;
-    yvect0[i] = dx;
-    yvect1[i] = dy;
-
-    xvect2[2*i+0] = sr[idx2[i]].X;
-    yvect2[2*i+0] = sr[idx2[i]].Y;
-    xvect2[2*i+1] = x;
-    yvect2[2*i+1] = y;
-  }
-
-  plot_resid_plot (0, xvect0, yvect0, Nvect);
-  plot_resid_plot (1, xvect1, yvect1, Nvect);
-  plot_fullfield_pairs (xvect2, yvect2, 2*Nvect);
-
-  free (xvect2);
-  free (yvect2);
-  free (xvect0);
-  free (yvect0);
-  free (xvect1);
-  free (yvect1);
-
-}
Index: /trunk/Ohana/src/gastro2/src/gpairs.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/gpairs.c	(revision 8241)
+++ /trunk/Ohana/src/gastro2/src/gpairs.c	(revision 8241)
@@ -0,0 +1,35 @@
+# include "gastro2.h"
+
+static int NPAIR, Npair;
+static int *idx1, *idx2;
+
+void pair_init () {
+
+  Npair = 0;
+  NPAIR = 100;
+  ALLOCATE (idx1, int, NPAIR);
+  ALLOCATE (idx2, int, NPAIR);
+
+}
+
+void pair_add (int i1, int i2) {
+
+  idx1[Npair] = i1;
+  idx2[Npair] = i2;
+
+  Npair ++;
+
+  if (Npair == NPAIR) {
+    NPAIR += 100;
+    REALLOCATE (idx1, int, NPAIR);
+    REALLOCATE (idx2, int, NPAIR);
+  }
+
+}
+
+int pair_lists (int **index1, int **index2) {
+
+  *index1 = idx1;
+  *index2 = idx2;
+  return (Npair);
+}
Index: /trunk/Ohana/src/gastro2/src/gproject2.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/gproject2.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/gproject2.c	(revision 8241)
@@ -4,7 +4,7 @@
 
   int i, N;
-  double X, Y, M, Moff;
+  double X, Y, P, Q, M, Moff;
   double XMIN, XMAX, YMIN, YMAX, MMIN, MMAX;
-  Coords *coords;
+  Coords TPtoSky, FPtoTP, *coords;
   StarData *in, *out;
 
@@ -15,5 +15,40 @@
   in     = Ref[0].stars;
   out    = Subset[0].stars;
+
   coords = &Target[0].coords;
+
+  /* create Tangent Plane to Sky transformation from input coords */
+  strcpy (TPtoSky.ctype, coords[0].ctype);
+  TPtoSky.crval1 = coords[0].crval1;
+  TPtoSky.crval2 = coords[0].crval2;
+
+  TPtoSky.cdelt1 = TPtoSky.cdelt2 = 1;
+  TPtoSky.crpix1 = TPtoSky.crpix2 = 0;
+  TPtoSky.pc1_1 = TPtoSky.pc2_2 = 1;
+  TPtoSky.pc1_2 = TPtoSky.pc2_1 = 0;
+  TPtoSky.Npolyterms = 0;
+  for (i = 0; i < 7; i++) {
+    TPtoSky.polyterms[i][0] = 0;
+    TPtoSky.polyterms[i][1] = 0;
+  }
+
+  /* create Focal Plane to Tangent Plane transformation from input coords */
+  strcpy (FPtoTP.ctype, "FP---LIN");
+  FPtoTP.crval1 = FPtoTP.crval2 = 0;
+
+  FPtoTP.cdelt1 = coords[0].cdelt1;
+  FPtoTP.cdelt2 = coords[0].cdelt2;
+  FPtoTP.crpix1 = coords[0].crpix1;
+  FPtoTP.crpix2 = coords[0].crpix2;
+  FPtoTP.pc1_1  = coords[0].pc1_1;
+  FPtoTP.pc1_2  = coords[0].pc1_2;
+  FPtoTP.pc2_1  = coords[0].pc2_1;
+  FPtoTP.pc2_2  = coords[0].pc2_2;
+
+  FPtoTP.Npolyterms = coords[0].Npolyterms;
+  for (i = 0; i < 7; i++) {
+    FPtoTP.polyterms[i][0] = coords[0].polyterms[i][0];
+    FPtoTP.polyterms[i][1] = coords[0].polyterms[i][1];
+  }
 
   Moff = Ref[0].Moff;
@@ -55,7 +90,12 @@
     }
 
+    /* get tangent-plane coordinates as well */
+    RD_to_XY (&P, &Q, in[i].R, in[i].D, &TPtoSky);
+
     out[N] = in[i];
     out[N].X = X;
     out[N].Y = Y;
+    out[N].P = P;
+    out[N].Q = Q;
     out[N].M = M;
     N++;
@@ -77,4 +117,11 @@
   REALLOCATE (Subset[0].stars, StarData, MAX (1, Subset[0].N));
 
+  /* get tangent-plane coords for target as well */
+  for (i = 0; i < Target[0].N; i++) {
+    XY_to_RD (&P, &Q, Target[0].stars[i].X, Target[0].stars[i].Y, &FPtoTP);
+    Target[0].stars[i].P = P;
+    Target[0].stars[i].Q = Q;
+  }
+
 }
 
Index: /trunk/Ohana/src/gastro2/src/greference2.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/greference2.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/greference2.c	(revision 8241)
@@ -68,8 +68,14 @@
   Yo = -0.5*NFIELD*NY;
 
-  catstats[0].RA[0] = +360.0;
-  catstats[0].RA[1] = -360.0;
-  catstats[0].DEC[0] = +90.0;
-  catstats[0].DEC[1] = -90.0;
+  if (ASCA) {
+      XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &Target[0].coords);
+      catstats[0].RA[0]  = R - 90.0;
+      catstats[0].RA[1]  = R + 90.0;
+      catstats[0].DEC[0] = MAX (-90.0, D - 90.0);
+      catstats[0].DEC[1] = MIN (+90.0, D + 90.0);
+      if (VERBOSE) fprintf (stderr, "asca region: %f - %f, %f - %f\n", 
+			    catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]);
+      return;
+  }
 
   for (x = 0; x <= 1.0; x += 0.5) {
@@ -89,4 +95,12 @@
     }
   }
+
+  /* is a pole in the image?  if so, include it... */
+  status = RD_to_XY (&X, &Y, 0.0, 90.0, &Target[0].coords);
+  if (status) catstats[0].DEC[1] = 90.0;
+
+  status = RD_to_XY (&X, &Y, 0.0, -90.0, &Target[0].coords);
+  if (status) catstats[0].DEC[0] = -90.0;
+
   if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", 
 	   catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]);
Index: /trunk/Ohana/src/gastro2/src/grid.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/grid.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/grid.c	(revision 8241)
@@ -1,5 +1,5 @@
 # include "gastro2.h"
 
-# define NPIX 50
+/* this needs to be user-configured */
 static int NX, NY, Nbin;
 static double C0x, C1x, C0y, C1y;
@@ -98,8 +98,8 @@
     if (PLOTSTUFF) plot_done_gridplot (dX, dY);
 
-    XMIN -= 0.5*NPIX;
-    XMAX -= 0.5*NPIX;
-    YMIN -= 0.5*NPIX;
-    YMAX -= 0.5*NPIX;
+    XMIN -= 0.5*NGRID_PIX;
+    XMAX -= 0.5*NGRID_PIX;
+    YMIN -= 0.5*NGRID_PIX;
+    YMAX -= 0.5*NGRID_PIX;
     gridfree ();
   }
@@ -111,6 +111,6 @@
 void gridinit (double XMIN, double XMAX, double YMIN, double YMAX, int Nr, int Nt) {
 
-  NX = (XMAX - XMIN) / NPIX;
-  NY = (YMAX - YMIN) / NPIX;
+  NX = (XMAX - XMIN) / NGRID_PIX;
+  NY = (YMAX - YMIN) / NGRID_PIX;
 
   C1x =          NX / (XMAX - XMIN); 
Index: /trunk/Ohana/src/gastro2/src/plots.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/plots.c	(revision 8240)
+++ /trunk/Ohana/src/gastro2/src/plots.c	(revision 8241)
@@ -169,5 +169,5 @@
 
   int i, Nvect;
-  float *xvect, *yvect, *zvect, M, dM;
+  float *xvect, *yvect, *zvect, M, dM, mRefMin, mRefMax;
   char c;
   
@@ -216,8 +216,17 @@
   ALLOCATE (yvect, float, Nvect);
   ALLOCATE (zvect, float, Nvect);
+
+  mRefMin = 30.0;
+  mRefMax = -5.0;
+  for (i = 0; i < Nvect; i++) {
+      mRefMin = MIN (mRefMin, Ref[0].stars[i].M);
+      mRefMax = MAX (mRefMax, Ref[0].stars[i].M);
+  }
+  dM = mRefMax - mRefMin;
+
   for (i = 0; i < Nvect; i++) {
     xvect[i] = Ref[0].stars[i].X;
     yvect[i] = Ref[0].stars[i].Y;
-    M = (Target[0].lum.Mmax - Ref[0].stars[i].M) / dM;
+    M = (mRefMax - Ref[0].stars[i].M) / dM;
     zvect[i] = MIN (1.0, MAX (0.01, M));
   }
@@ -330,4 +339,51 @@
 }
 
+
+void plot_resid (StarData *st, StarData *sr, Coords *coords) {
+
+  int i;
+  double x, y, dx, dy;
+  float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2;
+  int Npair, *idx1, *idx2;
+
+  Npair = pair_lists (&idx1, &idx2);
+  ALLOCATE (xvect0, float, Npair);
+  ALLOCATE (yvect0, float, Npair);
+  ALLOCATE (xvect1, float, Npair);
+  ALLOCATE (yvect1, float, Npair);
+
+  ALLOCATE (xvect2, float, 2*Npair);
+  ALLOCATE (yvect2, float, 2*Npair);
+  
+  for (i = 0; i < Npair; i++) {
+    RD_to_XY (&x, &y, sr[idx2[i]].R, sr[idx2[i]].D, coords);
+    
+    dx = x - st[idx1[i]].X;
+    dy = y - st[idx1[i]].Y;
+    
+    xvect0[i] = st[idx1[i]].X;
+    xvect1[i] = st[idx1[i]].Y;
+    yvect0[i] = dx;
+    yvect1[i] = dy;
+
+    xvect2[2*i+0] = st[idx1[i]].X;
+    yvect2[2*i+0] = st[idx1[i]].Y;
+    xvect2[2*i+1] = x;
+    yvect2[2*i+1] = y;
+  }
+
+  plot_resid_plot (0, xvect0, yvect0, Npair);
+  plot_resid_plot (1, xvect1, yvect1, Npair);
+  plot_fullfield_pairs (xvect2, yvect2, 2*Npair);
+
+  free (xvect2);
+  free (yvect2);
+  free (xvect0);
+  free (yvect0);
+  free (xvect1);
+  free (yvect1);
+
+}
+
 void fill_lumfunc (StarData *stars, int N, float *lbin, float *bin, int *nb) {
 
Index: /trunk/Ohana/src/gastro2/src/polyfit.c
===================================================================
--- /trunk/Ohana/src/gastro2/src/polyfit.c	(revision 8241)
+++ /trunk/Ohana/src/gastro2/src/polyfit.c	(revision 8241)
@@ -0,0 +1,409 @@
+# include "gastro2.h"
+
+static int NTERM, NPOWR, NPARS, NORDER, Npts;
+static double **sum, **xsum, **ysum;
+static double **matrix, **vector;
+
+void fit_init (int order) {
+
+  int i;
+
+  NORDER = order;
+  NPOWR = NORDER + 1;
+  NTERM = 2*NORDER + 1;
+  NPARS = (NORDER + 1)*(NORDER + 2) / 2;
+  Npts  = 0;
+
+  /* allocate arrays for fit solution */
+  ALLOCATE (sum, double *, NTERM);
+  ALLOCATE (xsum, double *, NTERM);
+  ALLOCATE (ysum, double *, NTERM);
+  for (i = 0; i < NTERM; i++) {
+    ALLOCATE (sum[i], double, NTERM);
+    bzero (sum[i], NTERM*sizeof(double));
+    ALLOCATE (xsum[i], double, NTERM);
+    bzero (xsum[i], NTERM*sizeof(double));
+    ALLOCATE (ysum[i], double, NTERM);
+    bzero (ysum[i], NTERM*sizeof(double));
+  }
+  ALLOCATE (matrix, double *, NPARS);
+  ALLOCATE (vector, double *, NPARS);
+  for (i = 0; i < NPARS; i++) {
+    ALLOCATE (matrix[i], double, NPARS);
+    ALLOCATE (vector[i], double, 2);
+    bzero (vector[i], 2*sizeof(double));
+    bzero (matrix[i], NPARS*sizeof(double));
+  }
+}
+
+/* */
+void fit_add (double x1, double y1, double x2, double y2, double wt) {
+
+  int n, m;
+  double xterm, yterm, term;
+
+  xterm = 1;
+  for (n = 0; n < NTERM; n++) {
+    yterm = 1;
+    for (m = 0; m < NTERM; m++) {
+      term = xterm*yterm;
+      if (n+m < NTERM) {
+	sum[n][m] += term;
+      }
+      if (n+m < NPOWR) {
+	xsum[n][m] += x2*term;
+	ysum[n][m] += y2*term;
+      }
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  Npts ++;
+
+}
+
+void fit_eval () {
+
+  int i, j, n, m, M, N;
+  double max;
+
+  i = 0;
+  for (m = 0; m < NPOWR; m++) {
+    for (n = 0; n < NPOWR - m; n++, i++) {
+      vector[i][0] = xsum[n][m];
+      vector[i][1] = ysum[n][m];
+    }	
+  }
+  j = 0;
+  for (M = 0; M < NPOWR; M++) {
+    for (N = 0; N < NPOWR - M; N++, j++) {
+      i = 0;
+      for (m = 0; m < NPOWR; m++) {
+	for (n = 0; n < NPOWR - m; n++, i++) {
+	  matrix[i][j] = sum[n+N][m+M];
+	}	
+      }
+    }
+  }       
+  max = 0.0;
+  for (i = 0; i < NPARS; i++) {
+    for (j = 0; j < NPARS; j++) {
+      max = MAX (max, fabs(matrix[i][j]));
+    }
+    max = MAX (max, fabs(vector[i][0]));
+    max = MAX (max, fabs(vector[i][1]));
+  }
+  for (i = 0; i < NPARS; i++) {
+    for (j = 0; j < NPARS; j++) {
+      matrix[i][j] /= max;
+    }
+    vector[i][0] /= max;
+    vector[i][1] /= max;
+  }
+  gaussj (matrix, NPARS, vector, 2); 
+  i = 0;
+  for (m = 0; m < NPOWR; m++) {
+    for (n = 0; n < NPOWR - m; n++, i++) {
+      xsum[n][m] = vector[i][0];
+      ysum[n][m] = vector[i][1];
+    }	
+  }
+  i = 0;
+  for (m = 0; m < NPOWR; m++) {
+    for (n = 0; n < NPOWR - m; n++, i++) {
+      if (VERBOSE) fprintf (stderr, "RA x^%dy^%d: %10.4g    DEC x^%dy^%d: %10.4g \n", n, m, vector[i][0], n, m, vector[i][1]);
+    }	
+  }
+}
+
+void fit_norm () { 
+
+  xsum[0][0] = 0;
+  xsum[1][0] = 1;
+  xsum[0][1] = 0;
+
+  ysum[0][0] = 0;
+  ysum[1][0] = 0;
+  ysum[0][1] = 1;
+}
+
+/* evaluate the fit at (X,Y) to yield (x,y) */
+void fit_apply (double *x, double *y, double X, double Y) {
+
+  int m, n;
+  double xterm, yterm;
+  double Xo, Yo;
+
+  Xo = Yo = 0;
+  yterm = 1;
+  for (m = 0; m < NPOWR; m++) { 
+    xterm = 1;
+    for (n = 0; n < NPOWR - m; n++) {
+      Xo += xterm*yterm*xsum[n][m];
+      Yo += xterm*yterm*ysum[n][m];
+      xterm *= X;
+    }	
+    yterm *= Y;
+  }
+  
+  *x = Xo;
+  *y = Yo;
+}
+
+/* evaluate the x-derivative of the fit at (X,Y) to yield (x,y) */
+void fit_apply_dx (double *x, double *y, double X, double Y) {
+
+  int m, n;
+  double xterm, yterm;
+  double Xo, Yo;
+
+  Xo = Yo = 0;
+  yterm = 1;
+  for (m = 0; m < NPOWR; m++) { 
+    xterm = 1;
+    for (n = 1; n < NPOWR - m; n++) {
+      Xo += n*xterm*yterm*xsum[n][m];
+      Yo += n*xterm*yterm*ysum[n][m];
+      xterm *= X;
+    }	
+    yterm *= Y;
+  }
+  
+  *x = Xo;
+  *y = Yo;
+}
+
+/* evaluate the y-derivative of the fit at (X,Y) to yield (x,y) */
+void fit_apply_dy (double *x, double *y, double X, double Y) {
+
+  int m, n;
+  double xterm, yterm;
+  double Xo, Yo;
+
+  Xo = Yo = 0;
+  yterm = 1;
+  for (m = 1; m < NPOWR; m++) { 
+    xterm = 1;
+    for (n = 0; n < NPOWR - m; n++) {
+      Xo += m*xterm*yterm*xsum[n][m];
+      Yo += m*xterm*yterm*ysum[n][m];
+      xterm *= X;
+    }	
+    yterm *= Y;
+  }
+  
+  *x = Xo;
+  *y = Yo;
+}
+
+/* measure the residual scatter in the fit */
+double fit_scat (StarData *st, StarData *sr, Coords *coords) {
+
+  int i;
+  int Npair, *idx1, *idx2;
+  double x, y, dx, dy, dX, dY, dX2, dY2, dR;
+  
+  Npair = pair_lists (&idx1, &idx2);
+
+  dX = dY = dX2 = dY2 = 0;
+  for (i = 0; i < Npair; i++) {
+
+    /* projection this direction includes the error introduced by
+       the interation on the nonlinear solution */
+    RD_to_XY (&x, &y, sr[idx2[i]].R, sr[idx2[i]].D, coords);
+    
+    dx = x - st[idx1[i]].X;
+    dy = y - st[idx1[i]].Y;
+    
+    dX += dx;
+    dY += dy;
+    dX2 += dx*dx;
+    dY2 += dy*dy;
+  }
+
+  /* scatter is measured on the tangent plane in degrees */
+  dX = dX / Npair;
+  dY = dY / Npair;
+  fprintf (stderr, "scatter: %f, %f\n", sqrt(dX2/Npair - dX*dX), sqrt(dY2/Npair - dY*dY));
+  fprintf (stderr, "Npts: %d\n", Npair);
+
+  dR = 0.5 * sqrt(fabs(dX2/Npair - dX*dX)) + 0.5 * sqrt (fabs(dY2/Npair - dY*dY));
+  return (dR);
+}
+
+/* convert fit terms to coords and polyterms */
+int fit_adjust (Coords *coords) {
+
+  int i, j, N;
+  double a00, a10, a01, a20, a11, a02, a30, a21, a12, a03, dA;
+  double b00, b10, b01, b20, b11, b02, b30, b21, b12, b03, dB;
+  double Xo, Yo, det;
+  int Np, Nv;
+  double **A, **B;
+    
+  /* start with the linear solution for Xo,Yo */
+  coords[0].cdelt1 = hypot (xsum[1][0], ysum[1][0]);
+  coords[0].cdelt2 = hypot (xsum[0][1], ysum[0][1]);
+
+  coords[0].cdelt1 = coords[0].cdelt2 = 1.0;
+  det = 1.0 / (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]);
+  Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]);
+  Yo = det*(xsum[0][0]*ysum[1][0] - ysum[0][0]*xsum[1][0]);
+
+  coords[0].Npolyterms = NORDER;
+
+  if (coords[0].Npolyterms > 1) {
+    /* use the linear solution as a starting guess */
+    /* solve for L(Xo,Yo) = 0, M(Xo,Yo) = 0 */
+    /* this is the Newton-Raphson method - it needs the high order terms to be small */
+    ALLOCATE (A, double *, 2);
+    ALLOCATE (B, double *, 2);
+    ALLOCATE (A[0], double, 2);
+    ALLOCATE (A[1], double, 2);
+    ALLOCATE (B[0], double, 1);
+    ALLOCATE (B[1], double, 1);
+
+    for (i = 0; i < 10; i++) {
+      fit_apply (&B[0][0], &B[1][0], Xo, Yo);
+      fit_apply_dx (&A[0][0], &A[0][1], Xo, Yo);
+      fit_apply_dy (&A[1][0], &A[1][1], Xo, Yo);
+      B[0][0] *= -1;
+      B[1][0] *= -1;
+      gaussj (A, 2, B, 1);
+      Xo += B[0][0]; 
+      Yo += B[1][0];
+    }
+    free (A[0]); free (B[0]);
+    free (A[1]); free (B[1]);
+    free (A); free (B);
+  }
+  coords[0].crpix1 = Xo;
+  coords[0].crpix2 = Yo;
+
+  switch (coords[0].Npolyterms) {
+    case 0:
+    case 1:
+      /* the linear solution can be analytically inverted */
+      coords[0].pc1_1 = xsum[1][0] / coords[0].cdelt1;
+      coords[0].pc1_2 = xsum[0][1] / coords[0].cdelt2;
+      coords[0].pc2_1 = ysum[1][0] / coords[0].cdelt1;
+      coords[0].pc2_2 = ysum[0][1] / coords[0].cdelt2;
+      for (i = 0; i < 7; i++) {
+	coords[0].polyterms[i][0] = coords[0].polyterms[i][1] = 0.0;
+      }
+      break;
+
+    case 2:
+      
+      a10 = xsum[1][0] + 2.0*xsum[2][0]*Xo + xsum[1][1]*Yo;
+      a01 = xsum[0][1] + 2.0*xsum[0][2]*Yo + xsum[1][1]*Xo;
+      a20 = xsum[2][0];
+      a11 = xsum[1][1];
+      a02 = xsum[0][2];
+
+      b10 = ysum[1][0] + 2.0*ysum[2][0]*Xo + ysum[1][1]*Yo;
+      b01 = ysum[0][1] + 2.0*ysum[0][2]*Yo + ysum[1][1]*Xo;
+      b20 = ysum[2][0];
+      b11 = ysum[1][1];
+      b02 = ysum[0][2];
+
+      coords[0].pc1_1 = a10 / coords[0].cdelt1;
+      coords[0].pc1_2 = a01 / coords[0].cdelt2;
+      coords[0].pc2_1 = b10 / coords[0].cdelt1;
+      coords[0].pc2_2 = b01 / coords[0].cdelt2;
+
+      coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1);
+      coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2);
+      coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2);
+
+      coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1);
+      coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2);
+      coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2);
+      for (i = 3; i < 7; i++) {
+	coords[0].polyterms[i][0] = coords[0].polyterms[i][1] = 0.0;
+      }
+      break;
+      
+    case 3:
+      /* the linear solution can be analytically inverted */
+      a10 = xsum[1][0] + 2*xsum[2][0]*Xo + xsum[1][1]*Yo + 3*xsum[3][0]*Xo*Xo + 2*xsum[2][1]*Xo*Yo + xsum[1][2]*Yo*Yo;
+      a01 = xsum[0][1] + 2*xsum[0][2]*Yo + xsum[1][1]*Xo + 3*xsum[0][3]*Yo*Yo + 2*xsum[1][2]*Xo*Yo + xsum[2][1]*Xo*Xo;
+      a20 = xsum[2][0] + 3*xsum[3][0]*Xo + xsum[2][1]*Yo;
+      a11 = xsum[1][1] +   xsum[2][1]*Xo + xsum[1][2]*Yo;
+      a02 = xsum[0][2] + 3*xsum[0][3]*Yo + xsum[1][2]*Xo;
+      a30 = xsum[3][0];
+      a21 = xsum[2][1];
+      a12 = xsum[1][2];
+      a03 = xsum[0][3];
+
+      b10 = ysum[1][0] + 2*ysum[2][0]*Xo + ysum[1][1]*Yo + 3*ysum[3][0]*Xo*Xo + 2*ysum[2][1]*Xo*Yo + ysum[1][2]*Yo*Yo;
+      b01 = ysum[0][1] + 2*ysum[0][2]*Yo + ysum[1][1]*Xo + 3*ysum[0][3]*Yo*Yo + 2*ysum[1][2]*Xo*Yo + ysum[2][1]*Xo*Xo;
+      b20 = ysum[2][0] + 3*ysum[3][0]*Xo + ysum[2][1]*Yo;
+      b11 = ysum[1][1] +   ysum[2][1]*Xo + ysum[1][2]*Yo;
+      b02 = ysum[0][2] + 3*ysum[0][3]*Yo + ysum[1][2]*Xo;
+      b30 = ysum[3][0];
+      b21 = ysum[2][1];
+      b12 = ysum[1][2];
+      b03 = ysum[0][3];
+
+      coords[0].pc1_1 = a10 / coords[0].cdelt1;
+      coords[0].pc1_2 = a01 / coords[0].cdelt2;
+      coords[0].pc2_1 = b10 / coords[0].cdelt1;
+      coords[0].pc2_2 = b01 / coords[0].cdelt2;
+
+      coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1);
+      coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2);
+      coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2);
+
+      coords[0].polyterms[3][0] = a30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1);
+      coords[0].polyterms[4][0] = a21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2);
+      coords[0].polyterms[5][0] = a12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1);
+      coords[0].polyterms[6][0] = a03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2);
+
+      coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1);
+      coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2);
+      coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2);
+
+      coords[0].polyterms[3][1] = b30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1);
+      coords[0].polyterms[4][1] = b21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2);
+      coords[0].polyterms[5][1] = b12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1);
+      coords[0].polyterms[6][1] = b03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2);
+      break;
+
+    default:
+      fprintf (stderr, "error: invalid order %d\n", coords[0].Npolyterms);
+      exit (2);
+  }
+
+  while (coords[0].crval1 < 0) coords[0].crval1 += 360.0;
+  while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
+
+  /* test for valid solution */
+  return (Npts);
+}
+
+int mk_polyterm (int n, int m, int norder) {
+  
+  int i, nt, N;
+  
+  N = 0;
+  nt = n + m;
+  for (i = 2; i < nt; i++) {
+    N += i + 1;
+  }
+  N += m;
+  return (N);
+}
+
+int mk_vector (int n, int m, int norder) {
+  
+  int i, N;
+  
+  N = 0;
+  for (i = 0; i < m; i++) {
+    N += (norder - i + 1);
+  }
+  N += n;
+  return (N);
+}
+