Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c	(revision 12068)
@@ -2,7 +2,6 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-FitPM (double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
 
-  PFfit *fit;
   double **A, **B;
 
@@ -57,6 +56,4 @@
   dgaussj (A, 4, B, 1);
 
-  ALLOCATE (fit, PMfit, 1);
-
   fit[0].Ro = B[0];
   fit[0].uR = B[1];
@@ -78,4 +75,4 @@
   free (B);
 
-  return (fit);
+  return (TRUE);
 }
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPMandPar.c	(revision 12068)
@@ -2,8 +2,6 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
 
-FitPM (double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
-
-  PFfit *fit;
   double **A, **B;
 
@@ -80,6 +78,4 @@
   dgaussj (A, 5, B, 1);
 
-  ALLOCATE (fit, PMfit, 1);
-
   fit[0].Ro = B[0];
   fit[0].uR = B[1];
@@ -101,4 +97,6 @@
   free (B);
 
-  return (fit);
+  /* get the chisq from the matrix values */
+
+  return (TRUE);
 }
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c	(revision 12068)
@@ -0,0 +1,29 @@
+# include "relastro.h"
+
+void FitSimple (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].R, ref[i].D);
+  }
+  fit_eval ();
+  fit_apply_coords (coords);
+}
+
+/* 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)
+*/
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c	(revision 12068)
@@ -12,4 +12,13 @@
 static Image        *image;
 static int          Nimage;
+
+Image *getimages (int *N) {
+  *N = Nimage;
+  return (image);
+}
+
+Image *getimage (int N) {
+  return (&image[N]);
+}
 
 void initImages (Image *input, int N) {
@@ -82,4 +91,5 @@
 }
 
+/* modify this function to use the measure->imageID field */
 void matchImage (Catalog *catalog, int meas, int cat) {
 
@@ -96,28 +106,4 @@
     if (measure[0].t > stop[i]) continue;
     
-# ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
-    if (USE_GRID) {
-
-      /* identify the ccd on the basis of the photcode name */
-      pname = GetPhotcodeNamebyCode (image[i].photcode);
-      filter = photcode[0].name;
-      sprintf (base, "%s.%s.", MOSAICNAME, filter);
-      if (strncmp (pname, base, strlen (base))) continue;
-      p = pname + strlen(base);
-      if (*p == 0) continue;
-      ccdnum = atoi (p);
-      /* some ambiguity here between seq number and id number */
-
-      /* add this measurement to the grid cell for this chip */
-      ave = measure[0].averef;
-      ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
-      dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
-       
-      /* X,Y always positive-definite in range 0,0 - dX, dY */
-      RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-      setGridMeasure (meas, cat, X, Y, ccdnum);
-    }
-# endif
-
     bin[cat][meas] = i;
 
@@ -133,18 +119,4 @@
     return;
   }
-  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
-}
-
-float getMcal (int meas, int cat) {
-
-  int i;
-  float value;
-
-  i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
-
-  if (image[i].code & IMAGE_BAD)  return (NO_MAG);  
-  value = image[i].Mcal;
-  return (value);
 }
 
@@ -156,128 +128,4 @@
   if (i == -1) return (NULL);
   return (&image[i].coords);
-}
-
-/* determine Mcal values for all images */
-void setMcal (Catalog *catalog, int PoorImages) {
-
-  int i, j, m, c, n, N, Nmax, mark, bad;
-  float Msys, Mrel, Mmos, Mgrid;
-  double *list, *dlist;
-  StatType stats;
-
-  if (FREEZE_IMAGES) return;
-
-  if (PoorImages) {
-    IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
-  }
-
-  Nmax = 0;
-  for (i = 0; i < Nimage; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  for (i = 0; i < Nimage; i++) {
-    
-    /* on PoorImages run, skip good images */
-    if (PoorImages) {
-      bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
-      if (!bad) continue;
-    }      
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measure[m].flags & MEAS_BAD) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-      
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-      list[N] = Msys - Mrel - Mmos - Mgrid;
-      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
-      N++;
-    }
-    /* Nlist[i] is all measurements, N is good measurements */
-
-    /* too few good measurements or too many bad measurements */
-    if (!PoorImages) {
-      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
-      if (mark) {
-	image[i].code |= ID_IMAGE_FEW;
-      } else {
-	image[i].code &= ~ID_IMAGE_FEW;
-      }      
-    }
-    
-    liststats (list, dlist, N, &stats);
-    image[i].Mcal  = stats.mean;
-    image[i].dMcal = stats.sigma;
-    image[i].Xm    = 100.0*log10(stats.chisq);
-  }
-  free (list);
-  free (dlist);
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
-    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
-  }
-  return;
-}
-
-/* mark image if: abs(Mcal) too large, dMcal too large */
-void clean_images () {
-
-  int i, N, mark, Nmark;
-  double *mlist, *slist, *dlist;
-  double MaxOffset, MaxScatter, MedOffset;
-  StatType stats;
-
-  if (FREEZE_IMAGES) return;
-
-  if (VERBOSE) fprintf (stderr, "marking poor images\n");
-
-  ALLOCATE (mlist, double, Nimage);
-  ALLOCATE (slist, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  for (i = N = 0; i < Nimage; i++) {
-    if (image[i].code & IMAGE_BAD) continue;
-    mlist[N] = fabs (image[i].Mcal);
-    slist[N] = image[i].dMcal;
-    dlist[N] = 1;
-    N++;
-  }
-  initstats ("MEAN");
-  liststats (mlist, dlist, N, &stats);
-  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
-  MedOffset = stats.median;
-  liststats (slist, dlist, N, &stats);
-  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
-  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
-  
-  Nmark = 0;
-  for (i = 0; i < Nimage; i++) {
-    mark = FALSE;
-    image[i].code &= ~ID_IMAGE_POOR;
-    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
-    if (mark) { 
-      Nmark ++;
-      image[i].code |= ID_IMAGE_POOR;
-    } else {
-      image[i].code &= ~ID_IMAGE_POOR;
-    }
-  }
-
-  fprintf (stderr, "%d images marked poor\n", Nmark);
-  initstats (STATMODE);
-  free (mlist);
-  free (slist);
-  free (dlist);
 }
 
@@ -330,134 +178,46 @@
 }
 
-StatType statsImageN (Catalog *catalog) {
-
-  int i, j, m, c, n, N;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-
-      m = mlist[i][j];
-      c = clist[i][j];
-
-      if (getMcal  (m, c) == NO_MAG) continue;
-      if (getMmos  (m, c) == NO_MAG) continue;
-      if (getMgrid (m, c) == NO_MAG) continue;
-      N++;
-    }
-    list[n] = N;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImageX (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = pow (10.0, 0.01*image[i].Xm);
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImageM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = image[i].Mcal;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImagedM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = image[i].dMcal;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-Image *getimages (int *N) {
-
-  *N = Nimage;
-  return (image);
-}
-
-Image *getimage (int N) {
-  return (&image[N]);
-}
-
+StarData *getImageRaw (int im, int *Nstars) {
+
+  StarData *raw;
+  
+  ALLOCATE (raw, StarData, Nlist[im]);
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    raw[i].X = catalog[c].measure[m].Xccd;
+    raw[i].Y = catalog[c].measure[m].Yccd;
+    
+    raw[i].Mag  = catalog[c].measure[m].M;
+    raw[i].dMag = catalog[c].measure[m].dM;
+
+    /* 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);
+  }
+  
+  *Nstars = Nlist[im];
+  return (raw);
+}
+
+StarData *getImageRef (int im, int *Nstars) {
+
+  StarData *raw;
+  
+  ALLOCATE (raw, StarData, Nlist[im]);
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+    n = catalog[c].measure[m].averef;
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    ref[i].R = catalog[c].average[n].R;
+    ref[i].D = catalog[c].average[n].D;
+  }
+  
+  *Nstars = Nlist[im];
+  return (raw);
+}
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/ParFactor.c	(revision 12068)
@@ -0,0 +1,60 @@
+# include "relastro.h"
+
+# if (0)
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+/* ra and dec are returned as decimal hours and decimal degrees. */
+void lpsun (double jd, double *ra, double *dec) {
+
+  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
+  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
+  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
+
+  // this is the conversion from ecliptic to celestial coords
+  x = cos(lambda);
+  y = cos(epsilon)*sin(lambda);
+  z = sin(epsilon)*sin(lambda);
+
+  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
+  *dec = (asin(z))*DEG_IN_RADIAN;
+}
+# endif
+
+/* code borrowed from Skycalc : fix this stuff XXX */
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
+
+  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)*RAD_DEG;
+  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
+  *beta = 0.0;					  // approx latitude
+  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
+  return TRUE;
+}
+
+/* given RA, DEC, Time, calculate the parallax factor */
+int ParFactor (double *pR, double *pD, double R, double D, time_t T) {
+
+  double jd;
+  double L, B;
+
+  /* given a time T in UNIX seconds, determine the solar longitude S */
+
+  jd = sec_to_jd (T);
+  sun_ecliptic (jd, &L, &B, &E);
+
+  e = E * RAD_DEG;
+  s = L * RAD_DEG;
+  r = R * RAD_DEG;
+  d = D * RAD_DEG;
+  
+  *pR =  +(cos(e)*sin(s)*cos(r) - cos(s)*sin(r));
+  *pD =  -(cos(e)*sin(s)*sin(r) + cos(s)*cos(r))*sin(d) + sin(e)*sin(s)*cos(d);
+  return TRUE;
+}
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c	(revision 12068)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+int UpdateChips () {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+
+  Image *image;
+  int Nimage;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (&image[i], &Nstars);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (&image[i], &Nstars);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitChip (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12068)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+int UpdateMosaic () {
+
+  /* we can measure new image parameters for each mosaic independently */
+
+  Image *image;
+  int Nimage;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip all except DIS images */
+    if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (&image[i], &Nstars);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (&image[i], &Nstars);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitMosaic (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c	(revision 12068)
@@ -2,9 +2,10 @@
 
 static int Nmax;
-static double *listX, *dlistX;
-static double *listY, *dlistY;
-static double *listR, *dlistR;
-static double *listD, *dlistD;
-static double *listT, *dlistT;
+static double *X, *dX;
+static double *Y, *dY;
+static double *R, *dR;
+static double *D, *dD;
+static time_t *T;
+static double *dT;
 
 void initObjectData (Catalog *catalog, int Ncatalog) {
@@ -19,6 +20,17 @@
   }
 
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
+  ALLOCATE (R, double, MAX (1, Nmax));
+  ALLOCATE (D, double, MAX (1, Nmax));
+  ALLOCATE (T, time_t, MAX (1, Nmax));
+  ALLOCATE (X, double, MAX (1, Nmax));
+  ALLOCATE (Y, double, MAX (1, Nmax));
+  ALLOCATE (pR, double, MAX (1, Nmax));
+  ALLOCATE (pD, double, MAX (1, Nmax));
+
+  ALLOCATE (dR, double, MAX (1, Nmax));
+  ALLOCATE (dD, double, MAX (1, Nmax));
+  ALLOCATE (dT, double, MAX (1, Nmax));
+  ALLOCATE (dX, double, MAX (1, Nmax));
+  ALLOCATE (dY, double, MAX (1, Nmax));
 }  
 
@@ -29,4 +41,5 @@
   StatType statsR, statsD;
   Coords coords;
+  PMFit fit;
 
   coords.crval1 = 0;
@@ -51,13 +64,13 @@
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
 	
-	listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	listT[N] = catalog[i].measure[m].t;
+	R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	T[N] = catalog[i].measure[m].t;
 
 	/* the astrometric errors are not being carried yet (but should be!) */
 	/* we use the photometric mag error as a weighting term */
-	dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	dlistT[N] = 1.0; /* XXX hard-wire this for now */
+	dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dT[N] = catalog[i].measure[m].dt;
 
 	N++;
@@ -65,4 +78,5 @@
       if (N < STAR_TOOFEW) { /* too few measurements */
 	catalog[i].average[j].code |= ID_STAR_FEW;
+	continue;
       } else {
 	catalog[i].average[j].code &= ~ID_STAR_FEW;
@@ -70,38 +84,58 @@
 
       /* we need to do the fit in a locally linear space; choose a ref coordinate */
-      coords.crval1 = listR[0];
-      coords.crval2 = listD[0];
+      coords.crval1 = R[0];
+      coords.crval2 = D[0];
       
       /* project all of the R,D coordinates to a plane centered on this coordinate */
       /* calculate pR[i], pD[i] for each point */
       for (k = 0; k < N; k++) {
-	  RD_to_XY (&listX[k], &listY[k], listR[k], listD[k], &coords);
-	  ParFactor (&pX[k], &pY[k], listR[k], listD[k]);
+	  RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
       }	  
 
-      /* in here, we should be fitting the parallax and proper-motion components */
-      /* XXX if I fit R,D,T as above, I will be getting the wrong answers near the pole */
-      if (FitProperMotion) {
-	  fit = FitPM (listX, dlistX, listY, dlistY, listT, pR, pD, N);
-      }
+      /* fit the model components as needed */
+      switch (FITTING_MODE) {
+	case FIT_AVERAGE:
+	  liststats (R, dR, N, &statsR);
+	  liststats (D, dD, N, &statsD);
 
-      if (FitProperMotion) {
-	  fit = FitPMandPar (listX, dlistX, listY, dlistY, listT, pR, pD, N);
-      }
+	  fit.Ro = statsR.mean;
+	  fit.dR = statsR.sigma;
 
-      if (AverageCoords) {
-	  liststats (listR, dlistR, N, &statsR);
-	  liststats (listD, dlistD, N, &statsD);
+	  fit.Do = statsD.mean;
+	  fit.dD = statsD.sigma;
 
-	  catalog[i].average[j].R = statsR.mean;
-	  catalog[i].average[j].dR = statsR.sigma;
+	  fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
+	  fit.Nfit = N;
+	  break;
 
-	  catalog[i].average[j].D = statsD.mean;
-	  catalog[i].average[j].dD = statsD.sigma;
+	case FIT_PM_ONLY:
+	  FitPM (&fit, X, dX, Y, dY, T, N);
+	  break;
 
-	  chisq = 0.5*(statsR.chisq + statsD.chisq);
-      }
+	case FIT_PM_AND_PAR:
+	  for (k = 0; k < N; k++) {
+	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
+	  }
+	  FitPMandPar (&fit, X, dX, Y, dY, T, pR, pD, N);
+	  break;
+	default:
+	  fprintf (stderr, "invalid fitting mode %d\n", FITTING_MODE);
+	  exit (2);
+      }	  
 
-      catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
+      catalog[i].average[j].R  	= fit.Ro;
+      catalog[i].average[j].D  	= fit.Do;
+      catalog[i].average[j].dR 	= fit.dRo;
+      catalog[i].average[j].dD 	= fit.dDo;
+
+      catalog[i].average[j].uR  = fit.uR;
+      catalog[i].average[j].uD  = fit.uD;
+      catalog[i].average[j].duR = fit.duR;
+      catalog[i].average[j].duD = fit.duD;
+
+      catalog[i].average[j].P   = fit.p;
+      catalog[i].average[j].dP  = fit.dp;
+
+      catalog[i].average[j].Xp  = (fit.Npts > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
     }
   }
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c	(revision 12068)
@@ -0,0 +1,37 @@
+# include "relastro.h"
+
+int UpdateSimple () {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+
+  Image *image;
+  int Nimage;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip WRP and DIS images */
+    if (!strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (&image[i], &Nstars);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (&image[i], &Nstars);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitSimple (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c	(revision 12068)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c	(revision 12068)
@@ -0,0 +1,284 @@
+# include "relastro.h"
+
+static int NTERM, NPOWER, NPARS, NORDER, Npts;
+static double **sum, **xsum, **ysum;
+static double **matrix, **vector;
+
+void 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;
+  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 < 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 () {
+
+  int i, j, n, m, M, N;
+
+  if (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
+
+  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
+}
+
+/* linear portion of fit : NORDER is 1 */
+void fit_apply_coords (Coords *coords) {
+
+  int i, j, Np, Nv, N;
+  double c11, c12;
+  double c21, c22;
+  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];
+      }
+    }
+  }
+
+  /* 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];
+
+  N = mkvector (1, 0, NORDER);
+  c11 = vector[N][0];  
+  c21 = vector[N][1];
+  N = mkvector (0, 1, NORDER);
+  c12 = vector[N][0];  
+  c22 = vector[N][1];
+  coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
+  R = 1 / coords[0].cdelt1;
+
+  coords[0].pc1_1  = c11*R;
+  coords[0].pc2_1  = c21*R;
+  coords[0].pc1_2  = c12*R;
+  coords[0].pc2_2  = c22*R;
+
+  coords[0].crpix1 = 0;
+  coords[0].crpix2 = 0;
+
+  coords[0].Npolyterms = NORDER;
+  strcpy (coords[0].ctype, "DEC--PLY");
+}
+
+/*
+  if we have just linear terms, the following holds for crpix1,2:
+  D = R / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+  coords[0].crpix1 = D * (coords[0].pc1_2*c20 - coords[0].pc2_2*c10);
+  coords[0].crpix2 = D * (coords[0].pc2_1*c10 - coords[0].pc1_1*c20);
+*/
+
+/* NORDER is order of gradient fit : Npolyterms is Norder + 1 */
+void fit_apply_grads (Coords *distort, Coords *project, int term) {
+
+  int i, j, Np, Nv1, Nv2;
+
+  distort[0].Npolyterms = NORDER + 1;
+  if (distort[0].Npolyterms > 1) strcpy (distort[0].ctype, "DEC--PLY");
+
+  for (i = 0; i < NPOWER + 1; i++) {
+    for (j = 0; j < (NPOWER + 1 - i); j++) {
+      
+      if (i + j < 2) continue;
+      Np  = mkpolyterm (i, j);
+      Nv1 = mkvector (i-1, j, NORDER);
+      Nv2 = mkvector (i, j-1, NORDER);
+
+      /** why do we have the negative sign? **/
+      if (j == 0) {
+	distort[0].polyterms[Np][term] = vector[Nv1][0] / i;
+      }
+      if (i == 0) {
+	distort[0].polyterms[Np][term] = vector[Nv2][1] / j;
+      }
+      if ((i > 0) && (j > 0)) {
+	distort[0].polyterms[Np][term] = 0.5*(vector[Nv1][0] / i + vector[Nv2][1] / j);
+      }
+    }
+  }
+
+  Nv1 = mkvector (0, 0, NORDER);
+  if (term == 0) {
+    project[0].pc1_1 = project[0].pc1_1 * (1 + vector[Nv1][0]);
+    project[0].pc1_2 = project[0].pc1_2 * (1 + vector[Nv1][1]);
+  } else {
+    project[0].pc2_1 = project[0].pc2_1 * (1 + vector[Nv1][0]);
+    project[0].pc2_2 = project[0].pc2_2 * (1 + vector[Nv1][1]);
+  }
+}
+
+void fit_correct_grads (Gradients *in, Gradients *out, int term) {
+
+  int i, k, m, n;
+  double x, y, dx, dy, dz1, dz2;
+
+  for (i = 0; i < in[0].Npts; i++) {
+    
+    dx = in[0].Lo[i];
+    dy = in[0].Mo[i];
+    dz1 = dz2 = 0.0;
+
+    k = 0;
+    x = y = 1;
+    for (m = 0; m < NPOWER; m++) {
+      x = y;
+      for (n = 0; n < NPOWER - m; n++, k++) {
+	dz1 += vector[k][0]*x;
+	dz2 += vector[k][1]*x;
+	x = x * dx / SCALE;
+      }
+      y = y * dy / SCALE;
+    }
+
+    out[0].Lo[i] = dx;
+    out[0].Mo[i] = dy;
+    if (term == 0) {
+      out[0].dPdL[i] = in[0].dPdL[i] - dz1;
+      out[0].dPdM[i] = in[0].dPdM[i] - dz2;
+    } else {
+      out[0].dQdL[i] = in[0].dQdL[i] - dz1;
+      out[0].dQdM[i] = in[0].dQdM[i] - dz2;
+    }
+  }
+}
+
+int mkvector (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 fit_free () {
+
+  int i;
+
+  for (i = 0; i < NTERM; i++) {
+    free (sum[i]);
+    free (xsum[i]);
+    free (ysum[i]);
+  }
+  free (sum);
+  free (xsum);
+  free (ysum);
+
+  for (i = 0; i < NPARS; i++) {
+    free (matrix[i]);
+    free (vector[i]);
+  }
+  free (matrix);
+  free (vector);
+}
+  
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/load_images.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/load_images.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/load_images.c	(revision 12068)
@@ -34,4 +34,7 @@
   initMosaics (subset, Nsubset);
   
+  /* unlock, if we can (else, unlocked below) */
+  if (!UPDATE) dvo_image_unlock (db); 
+
   return (skylist);
 }
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c	(revision 12068)
@@ -27,7 +27,4 @@
   skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
 
-  /* unlock, if we can (else, unlocked below) */
-  if (!UPDATE) dvo_image_unlock (&db); 
-
   /* load catalog data from region files */
   catalog = load_catalogs (skylist, &Ncatalog);
@@ -35,13 +32,8 @@
   /* match measurements with images, mosaics */
   initImageBins  (catalog, Ncatalog);
-  initMosaicBins (catalog, Ncatalog);
-  initGridBins   (catalog, Ncatalog);
-  initMrel (catalog, Ncatalog);
+  // initMosaicBins (catalog, Ncatalog);
 
   findImages (catalog, Ncatalog);
-  findMosaics (catalog, Ncatalog);  /* also sets Grid values */
-  SAVEPLOT = FALSE;
-
-  setExclusions (catalog, Ncatalog);
+  // findMosaics (catalog, Ncatalog);  /* also sets Grid values */
 
   if (PLOTSTUFF) {
@@ -54,6 +46,12 @@
     UpdateObjects (catalog, Ncatalog);
   }
-  if (UpdateImages) {
-    UpdateImages (catalog, Ncatalog);
+  if (UpdateSimple) {
+    UpdateSimple ();
+  }
+  if (UpdateChips) {
+    UpdateChips ();
+  }
+  if (UpdateMosaics) {
+    UpdateMosaics ();
   }
 
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/setExclusions.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/setExclusions.c	(revision 12067)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/setExclusions.c	(revision 12068)
@@ -1,4 +1,5 @@
 # include "relastro.h"
 
+/* XXX I think all of these, except for X,Y selection, are done / can be done in bcatalog */
 int setExclusions (Catalog *catalog, int Ncatalog) {