Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c	(revision 12205)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitChip.c	(revision 12205)
@@ -0,0 +1,34 @@
+# include "relastro.h"
+
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+
+  fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
+  }
+  fit_eval ();
+  fit_apply_coords (coords);
+  /* FitChip and FitSimple update the coords in different ways? maybe not... */
+}
+
+/* XXX the supporting fitting functions are not quite right
+   I probably need to apply the projection and then do the fit 
+   to L,M.  I also need to incorporate the transformations 
+   determined for psastro */
+
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c	(revision 12205)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitMosaic.c	(revision 12205)
@@ -0,0 +1,34 @@
+# include "relastro.h"
+
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+
+  fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    fit_add (raw[i].L, raw[i].M, ref[i].P, ref[i].Q);
+  }
+  fit_eval ();
+  fit_apply_coords (coords);
+  /* FitChip and FitSimple update the coords in different ways? maybe not... */
+}
+
+/* XXX the supporting fitting functions are not quite right
+   I probably need to apply the projection and then do the fit 
+   to L,M.  I also need to incorporate the transformations 
+   determined for psastro */
+
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c	(revision 12205)
@@ -8,5 +8,5 @@
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
-    fit_add (raw[i].X, raw[i].Y, ref[i].R, ref[i].D);
+    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
   }
   fit_eval ();
@@ -28,2 +28,6 @@
    P,Q -> X,Y (polynomial transformation)
 */
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c	(revision 12205)
@@ -178,9 +178,13 @@
 }
 
-StarData *getImageRaw (int im, int *Nstars) {
+StarData *getImageRaw (int im, int *Nstars, int isMosaic) {
 
   StarData *raw;
   
   ALLOCATE (raw, StarData, Nlist[im]);
+
+  if (isMosaic) {
+    mosaic = getMosaicForImage (im);
+  }
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -196,5 +200,13 @@
 
     /* note that for a Simple image, L,M = P,Q */
-    XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].X, raw[i].Y, &image[im].coords);
+    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
+    if (isMosaic) {
+      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, &mosaic);
+      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &mosaic);
+    } else {
+      raw[i].P = raw[i].L;
+      raw[i].Q = raw[i].M;
+      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
+    }
   }
   
@@ -203,9 +215,13 @@
 }
 
-StarData *getImageRef (int im, int *Nstars) {
-
-  StarData *raw;
-  
-  ALLOCATE (raw, StarData, Nlist[im]);
+StarData *getImageRef (int im, int *Nstars, int isMosaic) {
+
+  StarData *ref;
+  
+  ALLOCATE (ref, StarData, Nlist[im]);
+
+  if (isMosaic) {
+    mosaic = getMosaicForImage (im);
+  }
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -217,7 +233,18 @@
     ref[i].R = catalog[c].average[n].R;
     ref[i].D = catalog[c].average[n].D;
+
+    /* note that for a Simple image, L,M = P,Q */
+    RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
+    if (isMosaic) {
+      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, &mosaic);
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &mosaic);
+    } else {
+      ref[i].L = ref[i].P;
+      ref[i].M = ref[i].Q;
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+    }
   }
   
   *Nstars = Nlist[im];
-  return (raw);
-}
+  return (ref);
+}
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c	(revision 12205)
@@ -18,8 +18,8 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars);
+    raw = getImageRaw (&image[i], &Nstars, TRUE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars);
+    ref = getImageRef (&image[i], &Nstars, TRUE);
 
     /* XXX grab this code from mosastro */
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12205)
@@ -18,8 +18,8 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars);
+    raw = getImageRaw (&image[i], &Nstars, TRUE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars);
+    ref = getImageRef (&image[i], &Nstars, TRUE);
 
     /* XXX grab this code from mosastro */
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c	(revision 12205)
@@ -19,8 +19,8 @@
 
     /* convert measure coordinates to raw entries */
-    raw = getImageRaw (&image[i], &Nstars);
+    raw = getImageRaw (&image[i], &Nstars, FALSE);
 
     /* convert average coordinates to ref entries */
-    ref = getImageRef (&image[i], &Nstars);
+    ref = getImageRef (&image[i], &Nstars, FALSE);
 
     /* XXX grab this code from mosastro */
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c	(revision 12204)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c	(revision 12205)
@@ -1,127 +1,119 @@
 # include "relastro.h"
 
-static int NTERM, NPOWER, NPARS, NORDER, Npts;
-static double **sum, **xsum, **ysum;
-static double **matrix, **vector;
-
-void fit_init (int order) {
+/* these functions support simultaneous 2D fits to
+   x2 = \sum a_i,j x1^i y1^j 
+   y2 = \sum b_i,j x1^i y1^j 
+
+   the order of the fit (largest coefficient) is fixed to a single
+   value for both x1,y1 terms and for x2,y2 fits
+
+   the code is currently confusing because we limit to i+j <= order.
+   this could be cleaner if we used masks and allowed i <= order, j <= order
+*/
+
+// XXX define a fit structure and drop the file static variables?
+CoordFit *fit_init (int order) {
 
   int i;
 
-  Npts  = 0;
-  NORDER = order;
-  NPOWER = NORDER + 1;
-  NTERM = 2*NORDER + 1;
-  NPARS = (NORDER + 1)*(NORDER + 2) / 2;
-
-  /* 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));
-  }
-
-}
-
-# define SCALE 1.0
-void fit_add (double x1, double y1, double x2, double y2) {
-
-  int n, m;
+  ALLOCATE (fit, CoordFit, 1);
+
+  fit[0].Npts   = 0;
+  fit[0].Norder = order;
+  fit[0].Nterms = order + 1;
+  fit[0].Nsums  = 2*order + 1;
+  fit[0].Nelems = SQ(order + 1);
+
+  /* summing arrays for fit solution */
+
+  // xsum[i][j] holds \sum (x2 wt x1^i y1^j)
+  // ysum[i][j] holds \sum (y2 wt x1^i y1^j)
+  ALLOCATE (fit[0].xsum, double *, fit[0].Nterms);
+  ALLOCATE (fit[0].ysum, double *, fit[0].Nterms);
+  for (i = 0; i < fit[0].Nterms; i++) {
+    bzero (fit[0].xsum[i], fit[0].Nterms*sizeof(double));
+    bzero (fit[0].ysum[i], fit[0].Nterms*sizeof(double));
+  }      
+
+  // sum[i][j] holds \sum (wt x1^i y1^j)
+  ALLOCATE (fit[0].sum, double *, fit[0].Nsums);
+  for (i = 0; i < fit[0].Nsums; i++) {
+    ALLOCATE (fit[0].sum[i], double, fit[0].Nsums);
+    bzero (fit[0].sum[i], fit[0].Nsums*sizeof(double));
+  }
+
+  // matrix, vector hold the final linear system
+  ALLOCATE (fit[0].matrix, double *, fit[0].Nelems);
+  ALLOCATE (fit[0].vector, double *, fit[0].Nelems);
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ALLOCATE (fit[0].matrix[i], double, fit[0].Nelems);
+    ALLOCATE (fit[0].vector[i], double, fit[0].Nelems);
+    bzero (fit[0].matrix[i], fit[0].Nelems*sizeof(double));
+    bzero (fit[0].vector[i], 2*sizeof(double));
+  }
+  return (fit);
+}
+
+// XXX use implicit masks as below or explicit masks (with function to set?)
+// XXX eg, add a global mask to this file and 
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) {
+
+  int ix, iy;
   double xterm, yterm, term;
 
   xterm = 1;
-  for (n = 0; n < NTERM; n++) {
+  for (ix = 0; ix < fit[0].Nsums; ix++) {
     yterm = 1;
-    for (m = 0; m < NTERM; m++) {
-      term = xterm*yterm;
-      if (n+m < NTERM) {
-	sum[n][m] += term;
-      }
-      if (n+m < NPOWER) {
-	xsum[n][m] += x2*term;
-	ysum[n][m] += y2*term;
-      }
-      yterm *= y1/SCALE;
-    }
-    xterm *= x1/SCALE;
-  }
-  Npts ++;
-}
-
-/** I am renormalizing here by the max pivots to keep gaussj sane **
- ** would not be needed if the fit used scaled ind. variables **/
-void fit_eval () {
+    for (iy = 0; iy < fit[0].Nsums; iy++) {
+      term = xterm*yterm*wt;
+      fit[0].sum[ix][iy] += term;
+      if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) {
+	fit[0].xsum[ix][iy] += x2*term;
+	fit[0].ysum[ix][iy] += y2*term;
+      }
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  fit[0].Npts ++;
+}
+
+/* convert the xsum,ysum,sum terms into vector,matrix and solve */
+void fit_eval (CoordFit *fit) {
 
   int i, j, n, m, M, N;
 
-  if (Npts == 0) {
+  if (fit[0].Npts == 0) {
     fprintf (stderr, "warning: no valid pts\n");
   }
 
-  i = 0;
-  for (m = 0; m < NPOWER; m++) {
-    for (n = 0; n < NPOWER - m; n++, i++) {
-      vector[i][0] = xsum[n][m];
-      vector[i][1] = ysum[n][m];
-    }	
-  }
-  j = 0;
-  for (M = 0; M < NPOWER; M++) {
-    for (N = 0; N < NPOWER - M; N++, j++) {
-      i = 0;
-      for (m = 0; m < NPOWER; m++) {
-	for (n = 0; n < NPOWER - m; n++, i++) {
-	  matrix[i][j] = sum[n+N][m+M];
-	}	
-      }
-    }
-  }       
-# if (0)
-  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;
-  }
-# endif
+  /* remap the xsum,ysum terms into the vector */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    vector[i][0] = xsum[ix][iy];
+    vector[i][1] = ysum[ix][iy];
+
+    for (j = 0; j < fit[0].Nelems; j++) {
+      jx = j % fit[0].Nterms;
+      jy = j / fit[0].Nterms;
+      matrix[i][j] = sum[ix+jx][iy+jy];
+    }
+  }
 
   dgaussj (matrix, NPARS, vector, 2); 
 
-# if (0)
-  i = 0;
-  for (m = 0; m < NPOWER; m++) {
-    for (n = 0; n < NPOWER - 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]);
-    }	
-  }
-# endif
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+	     ix, iy, vector[i][0], ix, iy, vector[i][1]);
+  }	
 }
 
 /* linear portion of fit : NORDER is 1 */
-void fit_apply_coords (Coords *coords) {
+/* this should only apply to the polynomial, not the projection terms */
+/* compare with psastro supporting code */
+void fit_apply_coords (CoordFit *fit, Coords *coords) {
 
   int i, j, Np, Nv, N;
@@ -130,19 +122,26 @@
   double R;
 
-  /* update the higher order terms */
-  if (NORDER > 1) {
-    for (i = 0; i < NPOWER; i++) {
-      for (j = 0; j < (NPOWER - i); j++) {
-	if (i + j < 2) continue;
-	Np = mkpolyterm (i, j);
-	Nv = mkvector (i, j, NORDER);
-	coords[0].polyterms[Np][0] = vector[Nv][0];
-	coords[0].polyterms[Np][1] = vector[Nv][1];
-      }
-    }
-  }
+  /* I have L,M = fit(X,Y). set corresponding terms for coords */
+
+  // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ...
+  // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
+
+  coords[0].crpix1,crpix2 = CoordsGetCenter (fit);
+  modfit = CoordsSetCenter (fit, xo, yo);
+
+  // set pc1_1, pc1_2, pc2_1, pc2_2
+  fit_to_coordterms (coords, modfit, 0, 1);
+  fit_to_coordterms (coords, modfit, 1, 0);
+
+  fit_to_coordterms (coords, modfit, 0, 2);
+  fit_to_coordterms (coords, modfit, 1, 1);
+  fit_to_coordterms (coords, modfit, 2, 0);
+
+  fit_to_coordterms (coords, modfit, 0, 3);
+  fit_to_coordterms (coords, modfit, 1, 2);
+  fit_to_coordterms (coords, modfit, 2, 1);
+  fit_to_coordterms (coords, modfit, 3, 0);
 
   /* get the correct vector entries for the linear terms */
-  N = mkvector (0, 0, NORDER);
   coords[0].crval1 = vector[N][0];  
   coords[0].crval2 = vector[N][1];
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c	(revision 12205)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/mkpolyterm.c	(revision 12205)
@@ -0,0 +1,68 @@
+# include "mosastro.h"
+
+fit_to_coordterms (Coords *coords, CoordFit *fit, int nx, int ny) {
+
+  int nsum;
+
+  nsum = nx + ny;
+
+  switch (nsum) {
+    case 0:
+      xterm = &coords[0].crpix1;
+      yterm = &coords[0].crpix2;
+      return;
+    case 1:
+      if (ny == 0) {
+	xterm = &coords[0].pc1_1;
+	yterm = &coords[0].pc2_1;
+      } else {
+	xterm = &coords[0].pc1_2;
+	yterm = &coords[0].pc2_2;
+      }
+      return;
+    case 2:
+      xterm = &coords[0].polyterm[ny][0];
+      yterm = &coords[0].polyterm[ny][1];
+      return;
+    case 3:
+      xterm = &coords[0].polyterm[ny+3][0];
+      yterm = &coords[0].polyterm[ny+3][1];
+      return;
+    default:
+      fprintf (stderr, "programming error %s:%d\n", __func__, __line__);
+      exit (2);
+  }
+  fprintf (stderr, "programming error %s:%d\n", __func__, __line__);
+  exit (2);
+}
+
+int mkpolyterm (int n, int m) {
+  
+  int i, nt, N;
+  
+  N = 0;
+  nt = n + m;
+  for (i = 2; i < nt; i++) {
+    N += i + 1;
+  }
+  N += m;
+  return (N);
+}
+
+/*
+  coords uses an odd, rigid sequence for the coefficients:
+
+  x^0 y^0 : crpix1,2 (but note this is applied before higher order terms)
+
+  x^1 y^0 : pc1_1, pc2_1
+  x^0 y^1 : pc1_2, pc2_2
+
+  x^2 y^0 : polyterm[0][0,1]
+  x^1 y^1 : polyterm[1][0,1]
+  x^0 y^2 : polyterm[2][0,1]
+
+  x^3 y^0 : polyterm[3][0,1]
+  x^2 y^1 : polyterm[4][0,1]
+  x^1 y^2 : polyterm[5][0,1]
+  x^0 y^3 : polyterm[6][0,1]
+*/