Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile	(revision 37445)
@@ -88,4 +88,5 @@
 $(SRC)/AstromOffsetMapIO.$(ARCH).o    \
 $(SRC)/AstromOffsetMapOps.$(ARCH).o \
+$(SRC)/AstromOffsetMapUtils.$(ARCH).o \
 $(SRC)/dvo_set_skyregion.$(ARCH).o \
 $(SRC)/dvo_set_catdir.$(ARCH).o \
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/libdvo_astro.h
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/libdvo_astro.h	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/libdvo_astro.h	(revision 37445)
@@ -57,5 +57,5 @@
   float **dXv;
   float **dYv;
-  unsigned int ID;
+  unsigned int tableID;
   unsigned int imageID;
 } AstromOffsetMap;
@@ -63,7 +63,8 @@
 typedef struct {
   int Nmap;
-  AstromOffsetMap *map;
+  AstromOffsetMap **map;
   int *IDtoSeq;
-  int MaxID;
+  int MaxImageID;
+  int MaxTableID;
 } AstromOffsetTable;
 
@@ -164,3 +165,8 @@
 int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir);
 
+/* in AstromOffsetMapUtils.c */
+int AstromOffsetTableNewMap (AstromOffsetTable *table, int order, Image *image);
+int AstromOffsetTableMatchChips (Image *images, off_t Nimages, AstromOffsetTable *table);
+AstromOffsetTable *AstromOffsetTableInit();
+
 # endif
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c	(revision 37445)
@@ -1,7 +1,3 @@
 # include "dvo.h"
-
-/* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric
- * corrections for each chip.  
- */
 
 AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) {
@@ -49,5 +45,6 @@
   }
 
-  // convert the disk table format to internal format:
+  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
+  // Convert the disk array of maps to then internal format in a rich structure:
   AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap);
   free (map_disk);
@@ -77,4 +74,6 @@
   gfits_create_matrix (&header, &matrix);
 
+  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
+  // Convert the internal format in a rich structure into a disk array of maps.
   AstromOffsetMap_Disk_6x6 *map_disk = AstromOffsetTableToMap (table, &Nmap); 
 
@@ -105,25 +104,27 @@
 
   table->Nmap = Nmap;
-  ALLOCATE (table->map, AstromOffsetMap, Nmap);
+  ALLOCATE (table->map, AstromOffsetMap *, Nmap);
 
   // find the max value of imageID
-  int MaxID = 0;
+  int MaxTableID = 0;
+  int MaxImageID = 0;
   for (i = 0; i < Nmap; i++) {
-    MaxID = MAX(map_disk[i].imageID, MaxID);
+    MaxTableID = MAX(map_disk[i].tableID, MaxTableID);
+    MaxImageID = MAX(map_disk[i].imageID, MaxImageID);
   }
-  MaxID ++; // we want the outer bound, not the last value
+  table->MaxTableID = MaxTableID;
+  table->MaxImageID = MaxImageID;
 
   // generate the index and init values to -1
-  ALLOCATE (table->IDtoSeq, int, MaxID);
-  for (i = 0; i < MaxID; i++) {
+  ALLOCATE (table->IDtoSeq, int, MaxImageID + 1);
+  for (i = 0; i <= MaxImageID; i++) {
     table->IDtoSeq[i] = -1;
   }
-  table->MaxID = MaxID - 1;
 
   // assign the ID values
   for (i = 0; i < Nmap; i++) {
-    int ID = map_disk[i].imageID;
-    myAssert (table->IDtoSeq[ID] == -1, "oops, duplicate image IDs");
-    table->IDtoSeq[ID] = i;
+    int ImageID = map_disk[i].imageID;
+    myAssert (table->IDtoSeq[ImageID] == -1, "oops, duplicate image IDs");
+    table->IDtoSeq[ImageID] = i;
   }
 
@@ -131,27 +132,28 @@
   // full 6x6, saving some memory while doing the analysis)
   for (i = 0; i < Nmap; i++) {
-    table->map[i].Nx 	  = map_disk[i].Nx;
-    table->map[i].Ny 	  = map_disk[i].Ny;
-    table->map[i].ID 	  = map_disk[i].ID;
-    table->map[i].imageID = map_disk[i].imageID;
+    ALLOCATE (table->map[i], AstromOffsetMap, 1);
+
+    table->map[i][0].Nx 	  = map_disk[i].Nx;
+    table->map[i][0].Ny 	  = map_disk[i].Ny;
+    table->map[i][0].tableID = map_disk[i].tableID;
+    table->map[i][0].imageID = map_disk[i].imageID;
     
     // STORE THESE VALUES?
-    table->map[i].dY 	  = map_disk[i].dX;
-    table->map[i].dX 	  = map_disk[i].dY;
+    table->map[i][0].dY 	  = map_disk[i].dX;
+    table->map[i][0].dX 	  = map_disk[i].dY;
 
-    ALLOCATE (table->map[i].dXv, float *, map_disk[i].Nx);
-    ALLOCATE (table->map[i].dYv, float *, map_disk[i].Nx);
+    ALLOCATE (table->map[i][0].dXv, float *, map_disk[i].Nx);
+    ALLOCATE (table->map[i][0].dYv, float *, map_disk[i].Nx);
 
     for (j = 0; j < map_disk[i].Nx; j++) {
-      ALLOCATE (table->map[i].dXv[j], float, map_disk[i].Ny);
-      ALLOCATE (table->map[i].dYv[j], float, map_disk[i].Ny);
+      ALLOCATE (table->map[i][0].dXv[j], float, map_disk[i].Ny);
+      ALLOCATE (table->map[i][0].dYv[j], float, map_disk[i].Ny);
 
       for (k = 0; k < map_disk[i].Ny; k++) {
-	table->map[i].dXv[j][k] = map_disk[i].dXv[j][k];
-	table->map[i].dYv[j][k] = map_disk[i].dYv[j][k];
+	table->map[i][0].dXv[j][k] = map_disk[i].dXv[j][k];
+	table->map[i][0].dYv[j][k] = map_disk[i].dYv[j][k];
       }
     }
   }
-
   return table;
 }
@@ -167,16 +169,16 @@
   // full 6x6, saving some memory while doing the analysis)
   for (i = 0; i < table->Nmap; i++) {
-    map_disk[i].Nx 	 = table->map[i].Nx;
-    map_disk[i].Ny 	 = table->map[i].Ny;
-    map_disk[i].ID 	 = table->map[i].ID;
-    map_disk[i].imageID  = table->map[i].imageID;
+    map_disk[i].Nx 	 = table->map[i][0].Nx;
+    map_disk[i].Ny 	 = table->map[i][0].Ny;
+    map_disk[i].tableID	 = table->map[i][0].tableID;
+    map_disk[i].imageID  = table->map[i][0].imageID;
     
-    map_disk[i].dX 	 = table->map[i].dX;
-    map_disk[i].dY 	 = table->map[i].dY;
+    map_disk[i].dX 	 = table->map[i][0].dX;
+    map_disk[i].dY 	 = table->map[i][0].dY;
 
     for (j = 0; j < map_disk[i].Nx; j++) {
       for (k = 0; k < map_disk[i].Ny; k++) {
-	map_disk[i].dXv[j][k] = table->map[i].dXv[j][k];
-	map_disk[i].dYv[j][k] = table->map[i].dYv[j][k];
+	map_disk[i].dXv[j][k] = table->map[i][0].dXv[j][k];
+	map_disk[i].dYv[j][k] = table->map[i][0].dYv[j][k];
       }
     }
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37445)
@@ -17,5 +17,5 @@
     float ymo = y * map->dY - 0.5;
     int   ymi = floor(ymo);
-    ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
+    ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 (Ny must be > 1)
 
     float ymf = ymo - ymi;
@@ -29,5 +29,5 @@
     float xmo = x * map->dX - 0.5;
     int   xmi = floor(xmo);
-    xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
+    xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-2  (Nx must be > 1)
 
     float xmf = xmo - xmi;
@@ -41,10 +41,10 @@
   float xmo = x * map->dX - 0.5;
   int   xmi = floor(xmo);
-  xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
+  xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-1
   float xmf = xmo - xmi;
 
   float ymo = y * map->dY - 0.5;
   int   ymi = floor(ymo);
-  ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
+  ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2
   float ymf = ymo - ymi;
 
@@ -113,5 +113,7 @@
   for (i = 0; i < Nx*Ny; i++) {
     ALLOCATE (A[i], float, Nx*Ny);
+    memset (A[i], 0, sizeof(float)*Nx*Ny);
     ALLOCATE (B[i], float, 1);
+    memset (B[i], 0, sizeof(float));
   }    
 
@@ -120,4 +122,12 @@
   // for (int n = 1; n < Nx - 1; n++) {
   // for (int m = 1; m < Ny - 1; m++) {
+
+  if (0) {
+    FILE *fd = fopen ("stats.dump.txt", "w");
+    for (i = 0; i < Npts; i++) {
+      fprintf (fd, "%d %f %f %f\n", i, x[i], y[i], f[i]);
+    }
+    fclose(fd);
+  }
 
   // float Total = 0.0;
@@ -267,5 +277,7 @@
 	  int J = (ix + jx) + Nx * (iy + jy);
 	  A[J][I] = sA[jx][jy];
-	  // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
+	  if (abs(A[J][I]) > 1000) {
+	    fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, ix, iy, ix + jx, iy + jy, sA[jx][jy]);
+	  }
 	  // Sum += sA[jn][jm];
 	}
@@ -283,5 +295,5 @@
   for (i = 0; i < Nx*Ny; i++) {
     Empty[i] = 0;
-    if (A[i][i] == 0.0) {
+    if (fabs(A[i][i]) < 1.e-2) {
       Empty[i] = 1;
       for (j = 0; j < Nx*Ny; j++) {
@@ -292,4 +304,15 @@
       B[i][0] = 0.0;
     }
+  }
+
+  if (0) {
+    FILE *fd = fopen ("matrix.dat", "w");
+    for (i = 0; i < Nx*Ny; i++) {
+      for (j = 0; j < Nx*Ny; j++) {
+	fprintf (fd, "%10.4f ", A[i][j]);
+      }
+      fprintf (fd, " : %10.4f\n", B[i][0]);
+    }
+    fclose (fd);
   }
 
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/astrom_maps.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/astrom_maps.c	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/astrom_maps.c	(revision 37445)
@@ -24,3 +24,47 @@
 }
 
-int 
+int AstromOffsetTableNewMap (AstromOffsetTable *table, int order, Image *image) {
+
+  int Nmap = table->Nmap;
+
+  table->Nmap++;
+  REALLOCATE (table->map, AstromOffsetMap, table->Nmap);
+
+  int Nx = order;
+  int Ny = order;
+
+  if (image->imageID > table->MaxID) {
+    int oldMaxID = table->MaxID;
+    table->MaxID = image->imageID;
+    REALLOCATE (table->IDtoSeq, int, table->MaxID + 1);
+    for (i = oldMaxID + 1; i < table->MaxID + 1; i++) {
+      table->IDtoSeq[i] = -1;
+    }
+  }
+  myAssert (table->IDtoSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
+  table->IDtoSeq[image->imageID] = Nmap;
+  
+  table->map[Nmap].Nx      = Nx;
+  table->map[Nmap].Ny      = Ny;
+  table->map[Nmap].ID      = table->maxID; table->maxID ++;
+  table->map[Nmap].imageID = image->imageID;
+  
+  table->map[Nmap].dX = Nx / image->Nx;
+  table->map[Nmap].dY = Ny / image->Ny;
+
+  ALLOCATE (table->map[Nmap].dXv, float *, Nx);
+  ALLOCATE (table->map[Nmap].dYv, float *, Nx);
+
+  for (j = 0; Nx; j++) {
+    ALLOCATE (table->map[i].dXv[j], float, Ny);
+    ALLOCATE (table->map[i].dYv[j], float, Ny);
+
+    for (k = 0; k < Ny; k++) {
+      table->map[i].dXv[j][k] = 0.0;
+      table->map[i].dYv[j][k] = 0.0;
+    }
+  }
+  image[0].coords.offsetMap = &table->map[Nmap];
+  return TRUE;    
+}
+
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/mosaic_astrom.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/mosaic_astrom.c	(revision 37444)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/mosaic_astrom.c	(revision 37445)
@@ -207,4 +207,6 @@
 # endif
 
+// XXX Note that this links (images[i].coords.mosaic, images[i].parent) will break
+// if we reallocate the Image array in the middle of program
 int BuildChipMatch (Image *images, off_t Nimages) {