Index: trunk/Ohana/src/relphot/src/GridOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/GridOps.c	(revision 20535)
+++ trunk/Ohana/src/relphot/src/GridOps.c	(revision 20536)
@@ -19,7 +19,7 @@
 static int    **Ymeas;
 
-static int    **clist;
-static int    **mlist;
-static int     *Nlist;
+static int    **clist; // link from measurement on a cell to catalog containing measurement
+static int    **mlist; // link from measurement on a cell to measurement in a catalog
+static int     *Nlist; // list of measurements for each grid cell
 static int     *NLIST;
 
@@ -282,4 +282,239 @@
 }
 
+/* direct (non-iterative) solution for Mgrid values for all grid bins */
+void setMgridDirect (Catalog *catalog, int Ncatalog) {
+
+  int **gotstar, **gridmeas;
+  int i, j, k, N, Nmax, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
+  double **A, **B, *Mjx, *Wjx;
+  float Msys, Mrel, Mcal, Mmos, Merr, Wsys;
+  double Mj, Wj;
+  
+  if (!USE_GRID) return;
+
+  ALLOCATE (A, double *, Ngrid);
+  ALLOCATE (B, double *, Ngrid);
+  for (i = 0; i < Ngrid; i++) {
+    ALLOCATE (A[i], double, Ngrid);
+    ALLOCATE (B[i], double, 1);
+    memset (A[i], 0, Ngrid*sizeof(double));
+    memset (B[i], 0, sizeof(double));
+  }
+
+  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = 0;
+
+  ALLOCATE (gotstar, int *, Ncatalog);
+  for (i = 0; i < Ncatalog; i++) {
+    ALLOCATE (gotstar[i], int, catalog[i].Naverage);
+  }  
+
+  // set up gridmeas table : grid index for each measurement
+  ALLOCATE (gridmeas, int *, Ncatalog);
+  for (i = 0; i < Ncatalog; i++) {
+    ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure);
+    for (j = 0; j < catalog[i].Nmeasure; j++) {
+      gridmeas[i][j] = -1;
+    }
+  }  
+  for (i = 0; i < Ngrid; i++) {
+    for (j = 0; j < Nlist[i]; j++) {
+      int m, c;
+      m = mlist[i][j];
+      c = clist[i][j];
+      gridmeas[c][m] = i;
+    }
+  }
+
+  // as we loop over the grid cells, we need to accumulate the values of Wjx for each star
+  ALLOCATE (Wjx, double, Ngrid);
+  ALLOCATE (Mjx, double, Ngrid);
+
+  // accumulate the elements of the matrix equation.  We have an equation of the form: Ax = B
+  // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix,
+  // and B is an Ngrid vector.  For a cell A(x,y), we need the following elements from each
+  // star which touches grid cell (x):
+  // 
+  // Mj  : sum over all measurements of Msys / dMsys^2 
+  // Wj  : sum over all measurements of 1.0 / dMsys^2 
+  // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2
+  // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2
+
+  // Mjx and Wjx can be calculated by summing over all measurements which touch the cell 
+  // Mj requires looping over stars which touch (x)
+
+  // this is tricky because we need to know both the measurements which touch a cell
+  // and the stars for which any measurement touches a cell.  we need to accumulate
+  // sums for each star which touches as cell on both bases.
+
+  for (i = 0; i < Ngrid; i++) {
+    
+    for (j = 0; j < Ncatalog; j++) {
+      memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int));
+    }  
+
+    // we are looping over the stars, but doing so by looping over the set of measurements:
+    // every star which touches this grid cell has a measurement in Nlist[i]
+    for (j = 0; j < Nlist[i]; j++) {
+      
+      int mx, c, n, m0, Npts;
+
+      mx = mlist[i][j];
+      c  = clist[i][j];
+      n  = catalog[c].measure[mx].averef;
+      
+      // if we have already visited this star, skip the stuff below
+      if (gotstar[c][n]) continue;
+      gotstar[c][n] = TRUE;
+
+      // skip stars marked as BAD
+      if (catalog[c].average[n].code & STAR_BAD) {
+	Nrel ++;
+	continue;
+      }
+
+      m0 = catalog[c].average[n].measureOffset;
+
+      // we accumuate an entry for each cell
+      memset (Wjx, 0, Ngrid*sizeof(double));
+      memset (Mjx, 0, Ngrid*sizeof(double));
+      Npts = Mj = Wj = 0.0;
+
+      // if we have not yet visited this star, accumulate the Mj, Wj entries
+      for (k = 0; k < catalog[c].average[n].Nmeasure; k++) {
+
+	int m, Ng;
+
+	m = m0 + k;
+
+	// skip measurements marked as BAD
+	if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+	  Nbad ++;
+	  continue;
+	}
+
+	// skip images marked as BAD
+	Mcal = getMcal  (m, c);
+	if (isnan(Mcal)) {
+	  Ncal ++;
+	  continue;
+	}
+
+	// skip mosaics marked as BAD
+	Mmos = getMmos  (m, c);
+	if (isnan(Mmos)) {
+	  Nmos ++;
+	  continue;
+	}
+
+	// select the color- and airmass-corrected observed magnitude for this star
+	// XXX need to be able to turn off the color-correction until initial average mags are found
+	// Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+	Msys = PhotCat (&catalog[c].measure[m]);
+	if (isnan(Msys)) {
+	  Nsys++;
+	  continue;
+	}
+
+	// mag-error for this measurement
+	Merr =  MAX (catalog[c].measure[m].dM, MIN_ERROR);
+
+	// Wsys = 1.0 / SQ(Merr);
+	Wsys = 1.0;
+
+	Ng = gridmeas[c][m];
+	if (Ng == -1) continue;  // skip measurements which do not touch any cell
+
+	Mj += Msys * Wsys;  // we are only including measurements touching this cell
+	Wj += Wsys;  // we are only including measurements touching this cell
+	Npts ++;
+	Ngood ++;
+
+	Mjx[Ng] += Msys * Wsys;  // we are only including measurements touching cell (x)
+	Wjx[Ng] += Wsys;  // we are only including measurements touching cell (x)
+      }
+
+      // some stars will not have any valid measurements, skip these
+      if (Npts == 0) continue;
+
+      B[i][0] += Mj*Wjx[i]/Wj - Mjx[i];
+      A[i][i] -= Wjx[i];
+      for (k = 0; k < Ngrid; k++) {
+	A[i][k] += Wjx[i]*Wjx[k]/Wj;
+	// fprintf (stderr, "%3.0f ", Wjx[k]);
+      }
+    }
+  }
+
+  if (1) {
+
+    FILE *f;
+    Header header, theader;
+    Matrix matrix;
+
+    /* we are writing to this file */
+    f = fopen ("matrix.fits", "w");
+    if (f == (FILE *) NULL) { 
+      fprintf (stderr, "cannot open matrix.fits for output\n");
+      return;
+    }
+
+    /* save grid mag values */
+    gfits_init_header (&theader);
+    theader.Naxes = 2;
+    theader.Naxis[0] = Ngrid;
+    theader.Naxis[1] = Ngrid;
+    theader.bitpix   = -32;
+    gfits_create_Theader (&theader, "IMAGE");
+    gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX");
+    gfits_create_matrix  (&theader, &matrix);
+    for (i = 0; i < Ngrid; i++) {
+      for (j = 0; j < Ngrid; j++) {
+	gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]);
+      }
+    }
+    gfits_fwrite_header (f, &theader);
+    gfits_fwrite_matrix (f, &matrix);
+    gfits_free_matrix (&matrix);
+
+    gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE");
+    gfits_create_matrix  (&theader, &matrix);
+    for (i = 0; i < Ngrid; i++) {
+      for (j = 0; j < Ngrid; j++) {
+	gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]);
+      }
+    }
+    gfits_fwrite_header (f, &theader);
+    gfits_fwrite_matrix (f, &matrix);
+    gfits_free_matrix (&matrix);
+    fclose (f);
+  }
+
+  dgaussjordan (A, B, Ngrid, 1);
+
+  fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys);
+
+  for (i = 0; i < Ngrid; i++) {
+    gridM[i] = B[i][0];
+    gridS[i] = sqrt(A[i][i]);
+    gridN[i] = N;
+  }
+
+  free (Wjx);
+  free (Mjx);
+
+  for (i = 0; i < Ngrid; i++) {
+    free (A[i]);
+  }
+  free (A);
+  free (B);
+
+  for (i = 0; i < Ncatalog; i++) {
+    free (gotstar[i]);
+    free (gridmeas[i]);
+  }  
+  free (gotstar);
+  free (gridmeas);
+}
+
 /* determine Mgrid values for all grid bins */
 void setMgrid (Catalog *catalog) {
@@ -437,5 +672,5 @@
 void dump_grid () { 
 
-    int i, j, Nimage, Nbytes, Nformat;
+  int i, j, Nimage, Nbytes, Nformat;
   int *imlist;
   FILE *f;