Index: /trunk/Ohana/src/opihi/dvo/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/dvo/Makefile	(revision 3687)
+++ /trunk/Ohana/src/opihi/dvo/Makefile	(revision 3688)
@@ -57,4 +57,5 @@
 $(SDIR)/dmt.$(ARCH).o		  	\
 $(SDIR)/elixir.$(ARCH).o                \
+$(SDIR)/fitcolors.$(ARCH).o             \
 $(SDIR)/gcat.$(ARCH).o		  	\
 $(SDIR)/gimages.$(ARCH).o	  	\
Index: /trunk/Ohana/src/opihi/dvo/fitcolors.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/fitcolors.c	(revision 3688)
+++ /trunk/Ohana/src/opihi/dvo/fitcolors.c	(revision 3688)
@@ -0,0 +1,228 @@
+# include "dvo1.h"
+# define NMIN_PTS 100
+
+void init_catalog (Catalog *catalog);
+void free_catalog (Catalog *catalog, int Ncatalog);
+
+/* this function takes a photcode (and camera name?) and measures the  *
+ * chip-to-chip slopes for all DEP photcodes equiv to the PRI/SEC code */
+
+int fitcolors (int argc, char **argv) {
+  
+  int *list, Nlist, Ncode;
+  int i, k, m, NP1, NP2, NP, Np, Npts, NPTS;
+  int N1, N2, i1, i2, mode[4];
+  int Nsec, Ngraph, Nregions, Ncatalog, status;
+  char catdir[256], filename[256];
+  char *cmd, *outcmd, *camera;
+  double Radius;
+  double *M1, *M2;
+  float *out;
+  Catalog *catalog;
+  Graphdata graphsky;
+  RegionFile *regions;
+  PhotCode **codelist, *tcode, *code[4];
+  Vector *xvec, *yvec;
+  Buffer *buf;
+
+  /* defaults */
+  regions  = NULL;
+  catalog  = NULL;
+  codelist = NULL;
+  xvec = yvec = NULL;
+
+  /* find CATDIR in config system */
+  VarConfig ("CATDIR", "%s", catdir);
+
+  /* load photcode information */
+  if (!InitPhotcodes ()) goto escape;
+  Nsec = GetPhotcodeNsecfilt ();
+
+  /*** determine appropriate catalog files (find_regions) ***/
+  /* load data about plot windows */
+  Ngraph = 0;
+  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
+  /* find catalog files which overlap this region */
+  Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
+  regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
+
+  /* interpret command-line options */
+  /*** determine any restrictions (time, ?) ***/
+  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (argc != 4) goto usage;
+
+  /* vectors to save data */
+  Npts = 0;
+  NPTS = 1;
+  if ((xvec = SelectVector ("tmp_x", ANYVECTOR, TRUE)) == NULL) goto escape;
+  if ((yvec = SelectVector ("tmp_y", ANYVECTOR, TRUE)) == NULL) goto escape;
+
+  /* determine relevant photcodes, colors */
+  if (!(Np = GetPhotcodeCodebyName (argv[2]))) {
+    fprintf (stderr, "ERROR: photcode not found in photcode table\n");
+    return (FALSE);
+  }
+  camera = argv[3];
+
+  /* reduce the list of codes */
+  list = GetPhotcodeEquivList (Np, &Nlist);
+  ALLOCATE (codelist, PhotCode *, Nlist);
+  for (i = NP = 0; i < Nlist; i++) {
+    tcode = GetPhotcodebyCode (list[i]);
+    if (strncmp (tcode[0].name, camera, strlen(camera))) continue;
+    codelist[NP] = tcode;
+    NP++;
+  }
+  mode[0] = mode[1] = MAG_REL;  /* we should be applying any relative photometry corrections here */
+  mode[2] = mode[3] = MAG_AVE;
+  /* set the reference colors */
+  code[2] = GetPhotcodebyCode (codelist[0][0].c1);
+  code[3] = GetPhotcodebyCode (codelist[0][0].c2);
+  /* all codes must have the same colors (validate) */
+  for (i = 0; i < NP; i++) {
+    if (codelist[i][0].c1 != codelist[0][0].c1) goto color_mismatch;
+    if (codelist[i][0].c2 != codelist[0][0].c2) goto color_mismatch;
+  }
+  fprintf (stderr, "using %d photcodes\n", NP);
+
+  /* output is a named buffer */
+  if ((buf = SelectBuffer (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  buf[0].matrix.Naxis[0] = NP;
+  buf[0].matrix.Naxis[1] = NP;
+
+  /* I should encapsulate this in a create_default_buffer */
+  fits_free_matrix (&buf[0].matrix);
+  fits_free_header (&buf[0].header);
+  buf[0].header.bitpix = buf[0].bitpix = -32;
+  buf[0].header.unsign = buf[0].unsign = FALSE;
+  buf[0].header.bscale = buf[0].bscale = 1.0;
+  buf[0].header.bzero  = buf[0].bzero  = 0.0;
+  buf[0].header.Naxes = 2;
+  buf[0].header.Naxis[0] = NP;
+  buf[0].header.Naxis[1] = NP;
+  fits_create_header (&buf[0].header);
+  fits_create_matrix (&buf[0].header, &buf[0].matrix);
+  out = (float *) buf[0].matrix.buffer;
+  /* we set a default flag value of -1 */
+  for (i = 0; i < NP*NP; i++) {
+    out[i] = -1;
+  }
+
+  /* loop over regions, extract data for each region */
+  ALLOCATE (catalog, Catalog, Nregions);
+  for (k = 0; k < Nregions; k++) {
+    /* lock, load, unlock catalog */
+    sprintf (filename, "%s/%s", catdir, regions[k].name);
+    catalog[k].filename = filename;
+    switch (lock_catalog (&catalog[k], LCK_SOFT)) {
+      case 2:
+	unlock_catalog (&catalog[k]);
+      case 0:
+	catalog[k].Naverage = 0;
+	continue;
+    }
+    if (!load_catalog (&catalog[k], LOAD_AVES | LOAD_MEAS | LOAD_SECF, TRUE)) {
+      catalog[k].Naverage = 0;
+      unlock_catalog (&catalog[k]);
+      continue;
+    }
+    unlock_catalog (&catalog[k]);
+  }
+  fprintf (stderr, "using %d regions\n", Nregions);
+
+  /*** generate the color-color vectors for the pairs ***/
+  /* loop over chip photcode pairs */
+  for (NP1 = 0; NP1 < NP; NP1++) {
+    for (NP2 = NP1 + 1; NP2 < NP; NP2++) {
+      code[0] = codelist[NP1];
+      code[1] = codelist[NP2];
+      
+      /* extract all magnitude pairs from catalog tables */
+      Npts = 0;
+      for (k = 0; k < Nregions; k++) {
+	if (catalog[k].Naverage == 0) continue;
+
+	/* get correct mags, convert to X,Y */
+	for (i = 0; i < catalog[k].Naverage; i++) {
+	  M1 = M2 = NULL;
+	  m = catalog[k].average[i].offset;
+
+	  M1 = ExtractDMag (&code[0], &mode[0], &catalog[k].average[i], &catalog[k].secfilt[i*Nsec], &catalog[k].measure[m], &N1);
+	  if (N1 == 0) goto skip_star;
+
+	  M2 = ExtractDMag (&code[2], &mode[2], &catalog[k].average[i], &catalog[k].secfilt[i*Nsec], &catalog[k].measure[m], &N2);
+	  if (N2 == 0) goto skip_star;
+
+	  for (i1 = 0; i1 < N1; i1++) {
+	    for (i2 = 0; i2 < N2; i2++) {
+	      yvec[0].elements[Npts] = M1[i1];
+	      xvec[0].elements[Npts] = M2[i2];
+	      Npts++;
+	      if (Npts == NPTS) {
+		NPTS += 2000;
+		REALLOCATE (xvec[0].elements, float, NPTS);
+		REALLOCATE (yvec[0].elements, float, NPTS);
+	      }
+	    }
+	  }
+	skip_star:
+	  if (M1 != NULL) free (M1);
+	  if (M2 != NULL) free (M2);
+	}
+      }
+
+      if (Npts < NMIN_PTS) continue;
+      xvec[0].Nelements = Npts;
+      yvec[0].Nelements = Npts;
+
+      /* perform robust fit on dmag vs color */
+      cmd = strcreate ("fit tmp_x tmp_y 1 -clip 3 3 -quiet");
+      status = command (cmd, &outcmd);
+      if (outcmd != NULL) free (outcmd);
+      
+      /* do something useful with the results (stored in Cn, C0, C1, etc) */
+      fprintf (GetOutfile(), "%s - %s : ", code[0][0].name, code[1][0].name);
+      fprintf (GetOutfile(), "%7.4f %7.4f   %7.4f   ", 
+	       get_double_variable ("C0"), get_double_variable ("C1"), get_double_variable ("dC"));
+      fprintf (GetOutfile(), "%5s of %5d\n", get_variable ("Cnv"), Npts);
+      out[NP1 + NP2*NP] = get_double_variable ("C1");
+    }
+  }
+  return (TRUE);
+
+usage:
+  fprintf (stderr, "USAGE: chipcolors (output) (photcode) (camera)\n");
+  goto escape;
+
+color_mismatch:
+  fprintf (stderr, "error: all chips must have the same colors\n");
+
+escape:
+  if (regions != NULL) free (regions);
+  if (codelist != NULL) free (codelist);
+  DeleteVector (xvec);
+  DeleteVector (yvec);
+  free_catalog (catalog, Ncatalog);
+}
+
+void free_catalog (Catalog *catalog, int Ncatalog) {
+
+  int i;
+
+  if (catalog == NULL) return;
+  for (i = 0; i < Ncatalog; i++) {
+    if (catalog[i].average != NULL) free (catalog[i].average);
+    if (catalog[i].secfilt != NULL) free (catalog[i].secfilt);
+    if (catalog[i].measure != NULL) free (catalog[i].measure);
+  }
+  free (catalog);
+}
+
+void init_catalog (Catalog *catalog) {
+  catalog[0].average = NULL; 
+  catalog[0].secfilt = NULL;
+  catalog[0].measure = NULL;
+  catalog[0].Naverage = 0; 
+  catalog[0].Nsecfilt = 0;
+  catalog[0].Nmeasure = 0;
+}
Index: /trunk/Ohana/src/opihi/dvo/init.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/init.c	(revision 3687)
+++ /trunk/Ohana/src/opihi/dvo/init.c	(revision 3688)
@@ -18,4 +18,5 @@
 int dmt		    PROTO((int, char **));
 int elixir          PROTO((int, char **));
+int fitcolors       PROTO((int, char **));
 int gcat	    PROTO((int, char **));
 int getxtra	    PROTO((int, char **));
@@ -69,4 +70,5 @@
   {"dmt",         dmt,          "plot mag scatter"},
   {"elixir",      elixir,       "talk to elixir"},
+  {"fitcolors",   fitcolors,    "fit chip-to-chip color terms"},
   {"gcat",    	  gcat,         "get catalog at location"},
   {"gimages", 	  gimages,      "get images at location"},
