Index: /trunk/Ohana/src/uniphot/src/cam_zpt_correction.c
===================================================================
--- /trunk/Ohana/src/uniphot/src/cam_zpt_correction.c	(revision 39263)
+++ /trunk/Ohana/src/uniphot/src/cam_zpt_correction.c	(revision 39263)
@@ -0,0 +1,272 @@
+# include "setphot.h"
+
+// camera systematic (photoflat) correction functions:
+// load the correction set from a file
+
+// note that we are allocating pointers for XY00 - XY77, but only the octal elements are set (and not the 00,07,70,77 ones)
+
+static CamCorrection *cam = NULL;
+
+CamCorrection *get_cam_correction_ptr () {
+  return cam;
+}
+
+void free_cam_correction (CamCorrection *cam) {
+
+  if (!cam) return;
+
+  for (i = 0; i < cam->Nvalues; i++) {
+    FREE (cam->matrix[i]);
+  }
+
+  free (cam->matrix);
+  free (cam);
+}
+
+CamCorrection *alloc_cam_correction (int Nx, int Ny, int Nfilter, int Nseason) {
+
+  CamCorrection *cam;
+
+  ALLOCATE (cam, CamCorrection, 1);
+
+  ALLOCATE (cam->matrix, Matrix *, Nseason*Nfilter*Nx*Ny);
+
+  cam->Nx      = Nx;  	   // for gpc1, should be 8
+  cam->Ny      = Ny;  	   // for gpc1, should be 8
+  cam->Nfilter = Nfilter;  // for gpc1, should be 5
+  cam->Nseason = Nseason;  // for gpc1, should be 5 (1 for loading highres correction)
+
+  cam->Nchips  = Nx * Ny;
+  cam->Nflats  = Nx * Ny * Nfilter;
+  cam->Nvalues = Nx * Ny * Nfilter * Nseason;
+
+  ALLOCATE (cam->matrix, Matrix *, cam->Nvalues);
+  for (i = 0; i < cam->Nvalues; i++) {
+    cam->matrix[i] = NULL;
+  }
+
+  return cam;
+}
+
+int load_cam_correction (char *filename) {
+
+  int i, ix, iy, dir, filter;
+  Header header;
+
+  /* open file for input */
+  FILE *f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    gprint (GP_ERR, "can't open file for read : %s\n", filename);
+    return FALSE;
+  }
+
+  // read the PHU header
+  if (!gfits_load_header (f, &header)) return FALSE;
+  
+  int Nbytes = gfits_data_size (&header);
+  fseeko (f, Nbytes, SEEK_CUR);
+
+  int Ndir, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD;
+  if (!gfits_scan (&header, "NSEASON",  "%d", 1, &Nseason)) return FALSE;
+  if (!gfits_scan (&header, "NFILTER",  "%d", 1, &Nfilter)) return FALSE;
+  if (!gfits_scan (&header, "NX", 	"%d", 1, &Nx)) 	    return FALSE;
+  if (!gfits_scan (&header, "NY", 	"%d", 1, &Ny)) 	    return FALSE;
+  if (!gfits_scan (&header, "DX", 	"%d", 1, &dX)) 	    return FALSE;
+  if (!gfits_scan (&header, "DY", 	"%d", 1, &dY)) 	    return FALSE;
+
+  if (!gfits_scan (&header, "NX_CHIP", "%d", 1, &NxCCD))   return FALSE;
+  if (!gfits_scan (&header, "NY_CHIP", "%d", 1, &NyCCD))   return FALSE;
+
+  CamCorrection *cam = alloc_cam_correction (Nx, Ny, Nfilter, Nseason);
+
+  cam->dX      = dX;  	   // superpixel sampling, x-dir
+  cam->dY      = dY;  	   // superpixel sampling, y-dir
+
+  cam->NxCCD   = NxCCD;    // pixels per chip
+  cam->NyCCD   = NyCCD;    // pixels per chip
+
+  while (TRUE) {
+    Header theader;
+
+    // load data for this header : if not found, assume we hit the end of the file
+    if (!gfits_load_header (f, &theader)) break;
+    
+    if (!gfits_scan (&theader, "FILTER", "%d", 1, &filter)) return FALSE;
+    if (!gfits_scan (&theader, "SEASON", "%d", 1, &season)) return FALSE;
+
+    // XXX NOTE: astroflat.20150209.fits had ix and iy backwards in header
+    // double-check that the new flats are OK
+    
+    if (!gfits_scan (&theader, "X_CHIP",  "%d", 1, &ix))      return FALSE;
+    if (!gfits_scan (&theader, "Y_CHIP",  "%d", 1, &iy))      return FALSE;
+
+    Matrix *matrix = NULL;
+    ALLOCATE (matrix, Matrix, 1);
+    if (!gfits_load_matrix (f, matrix, &theader)) break;
+
+    int index = ix + iy*cam->Nx + filter*cam->Nchips + season*cam->Nflats;
+    myAssert (index >= 0, "index too small");
+    myAssert (index < cam->Nvalues, "index too big");
+
+    myAssert (!cam->matrix[index], "entry already assigned?");
+
+    // assert that cam->matrix[index] is NULL?
+    cam->matrix[index] = matrix;
+  }
+  return TRUE;
+}
+
+// input is:
+// chipID == octal chip ID (eg 40 for XY40)
+// filter == (01234 = grizy)
+int get_cam_correction (int chipID, int filter, float Xccd, float Yccd, double *dX, double *dY) {
+
+  *dX = 0.0;
+  *dY = 0.0;
+
+  int index, jx, jy, NxModel;
+  Matrix *matrix;
+  float *buffer, value;
+
+  // split out the chipID number (eg 43 for XY43) into X and Y elements
+  int ix = (int) (chipID / 10);
+  int iy = (int) (chipID % 10);
+
+  // dX (dir == 0)
+  index = ix + iy*cam->Nx + filter*cam->Nchips;
+
+  matrix = cam->matrix[index];
+  NxModel = cam->matrix[index]->Naxis[0];
+
+  buffer = (float *) matrix->buffer;
+
+  jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0);
+  jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0);
+
+  value = buffer[jx + jy*NxModel];
+
+  *dX = isnan(value) ? 0.0 : value;
+
+  // dY (dir == 1) uses the next block
+  index += cam->Ngroup;
+
+  matrix = cam->matrix[index];
+  NxModel = cam->matrix[index]->Naxis[0];
+
+  buffer = (float *) matrix->buffer;
+
+  jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0);
+  jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0);
+
+  value = buffer[jx + jy*NxModel];
+
+  *dY = isnan(value) ? 0.0 : value;
+
+  return TRUE;
+}
+
+float get_cam_correction_by_flatid (int flat_id, float Xccd, float Yccd) {
+
+  float dM = 0.0;
+
+  // validate the flat_id (not out of range?)
+  int seq = camcorr->IDtoSeq[flat_id];
+
+  Matrix *matrix = cam->matrix[seq];
+  NxModel = cam->matrix[index]->Naxis[0];
+
+  buffer = (float *) matrix->buffer;
+
+  jx = MAX(MIN(cam->NxCCD, Xccd / cam->dX), 0);
+  jy = MAX(MIN(cam->NyCCD, Yccd / cam->dY), 0);
+
+  value = buffer[jx + jy*NxModel];
+
+  dM = isnan(value) ? 0.0 : value;
+
+  return dM;
+}
+
+int merge_flatcorr_and_camcorr (FlatCorrectionTable *flatcorr, CamCorrection *camcorr) {
+
+  // this function is specifically tuned for merging the ubercal flats with the high-res flats
+
+  // we have an ubercal FlatCorrectionTable with Nseasons, Nfilter, Nx, Ny entries
+  // we have a highres flat with Nfilter, Nx, Ny images
+
+  myAssert (flatcorr->Nfilter <= camcorr->Nfilter, "unclear how to use more flatcorr filters");
+  myAssert (camcorr->Nseason == 1, "no code to merge seasons flats to seasonal flats");
+
+  // we need to create a new CamCorrection with Nx,Ny,Nfilter from camcorr and Nseason from flatcorr:
+  CamCorrection *newcorr = alloc_cam_correction (camcorr->Nx, camcorr->Ny, camcorr->Nfilter, flatcorr->Nseason);
+  newcorr->dX      = camcorr->dX;    // superpixel sampling, x-dir
+  newcorr->dY      = camcorr->dY;    // superpixel sampling, y-dir
+  newcorr->NxCCD   = camcorr->NxCCD; // pixels per chip
+  newcorr->NyCCD   = camcorr->NyCCD; // pixels per chip
+
+  for (Ns = 0; Ns < flatcorr->Nseason; Ns++) {
+    newcorr->tstart[Ns] = flatcorr->tstart[Ns];
+    newcorr->tstop[Ns]  = flatcorr->tstop[Ns];
+  }
+
+  // create the new flat images (1 per season)
+  for (Ns = 0; Ns < newcorr->Nseason; Ns++) {
+    for (Nf = 0; Nf < newcorr->Nfilter; Nf++) {
+      for (ix = 0; ix < newcorr->Nx; ix++) {
+	for (iy = 0; iy < newcorr->Ny; iy++) {
+	  int index = ix + iy*newcorr->Nx + Nf*newcorr->Nchips + Ns*newcorr->Nflats;
+	  ALLOCATE (newcorr->matrix[index], Matrix, 1);
+	  
+	  // camcorr only had one season
+	  int idxOld = ix + iy*camcorr->Nx + Nf*camcorr->Nchips;
+
+	  // allocates the buffer 
+	  gfits_copy_matrix (camcorr->matrix[idxOld], newcorr->matrix[index]);
+	}
+      }
+    }
+  }
+
+  // FlatCorrectionTable has:
+  // table->image[seq] with metadata about the specific image
+  // table->offset[seq][ix][iy] with dM values for image[seq] 
+
+  // loop over the flatcorr images and match from image[seq] to newcorr:
+
+  for (i = 0; i < flatcorr->Nimage; i++) {
+
+    int ID = flatcorr->image[i].ID;
+    myAssert (flatcorrTable->IDtoSeq[ID] == i, "oops");
+
+    int season = -1;
+    for (j = 0; j < newcorr->Nseason; j++) {
+      if (flatcorr->image[i].tstart != newcorr->tstart[j]) continue;
+      season = j;
+    }
+    myAssert (season >= 0, "oops");
+
+    int photcode = flatcorr->image[i].photcode;
+
+    iy = flatcorr->image[i].photcode % 10;
+    ix = (int)(flatcorr->image[i].photcode / 10) % 10;
+    filter = (int)(flatcorr->image[i].photcode / 100) % 10;
+	      
+    int seq = ix + iy * newcorr->Nx + filter * newcorr->Nchips + season * camcorr->Nflat;
+
+    // we now have newcorr->matrix[seq]
+    newcorr->matrix[seq];
+
+    // now loop over the newcorr pixels to get the flatcorr value:
+    for (jx = 0; jx < newcorr->NxCCD; jx++) {
+      for (jy = 0; jy < newcorr->NyCCD; jy++) {
+	
+	// convert jx,jy in newcorr pixels to Jx,Jy in flatcorr pixels
+	newcorr->matrix[seq]->buffer[jx + jy*NxCCD] += flatcorr->offset[i][Jx][Jy];
+
+      }
+    }
+  }    
+
+  return TRUE;
+}
+
Index: /trunk/Ohana/src/uniphot/src/match_camcorr_to_images.c
===================================================================
--- /trunk/Ohana/src/uniphot/src/match_camcorr_to_images.c	(revision 39263)
+++ /trunk/Ohana/src/uniphot/src/match_camcorr_to_images.c	(revision 39263)
@@ -0,0 +1,53 @@
+# include "setphot.h"
+
+// the date/time of the image is used to find the 'season'
+// the photcode is used to find the actual flat correction
+
+int match_camcorr_to_images (Image *image, off_t Nimage) {
+
+  int i, j;
+
+  CamCorrection *camcorr = get_cam_correction_ptr ();  
+
+  // we have an array of CamCorrection->matrix values, where each matrix is a flat-field image
+  // the sequence is seq = ix + Nx_chip*iy + filter*Nchips + season*(Nchips*Nfilter)
+
+  // we derive ix,iy,filter from photcode : 10143 -> ix = 4, iy = 3, filter = 1
+  // season comes from the date/time
+
+  int minCode = 10000;
+  int maxCode = 10576;
+
+  for (i = 0; i < Nimage; i++) {
+    if (!image[i].photcode) continue; // skip PHU images
+
+    if (image[i].photcode < minCode) continue; // skip non-gpc1-chip images
+    if (image[i].photcode > maxCode) continue; // skip non-gpc1-chip images
+
+    int found = FALSE;
+    for (j = 0; !found && (j < camcorr->Nseason); j++) {
+      if (image[i].tzero < camcorr->tstart[j]) continue;
+      if (image[i].tzero > camcorr->tstop[j]) continue;
+
+      int iy = image[i].photcode % 10;
+      int ix = (int)(image[i].photcode / 10) % 10;
+      int filter = (int)(image[i].photcode / 100) % 10;
+
+      int seq = ix + iy * camcorr->Nx + filter * camcorr->Nchips + j * camcorr->Ngroup;
+
+      image[i].photom_map_id = camcorr->imageID[seq].ID;
+      found = TRUE;
+    }
+    if (!found) {
+      Nmissed ++;
+      fprintf (stderr, "image does not match flatcorr Table ranges: %s\n", image[i].name);
+    }
+
+    if (Nmissed > 100)  {
+      fprintf (stderr, "something is wrong: too many images do not match\n");
+      exit (2);
+    }
+
+  }
+  return TRUE;
+}