Index: /trunk/Ohana/src/addstar/include/sedstar.h
===================================================================
--- /trunk/Ohana/src/addstar/include/sedstar.h	(revision 7692)
+++ /trunk/Ohana/src/addstar/include/sedstar.h	(revision 7693)
@@ -1,2 +1,13 @@
+
+typedef structu {
+  int Nfilter;
+  float *wavecode;
+  float *vegaToAB;
+  int *hashcode;
+  int codeP;
+  int codeM;
+  SEDtableRow **row;
+  int Nrow;
+} SEDtable;
 
 typedef struct {
Index: /trunk/Ohana/src/addstar/src/SEDfit.c
===================================================================
--- /trunk/Ohana/src/addstar/src/SEDfit.c	(revision 7693)
+++ /trunk/Ohana/src/addstar/src/SEDfit.c	(revision 7693)
@@ -0,0 +1,352 @@
+# include "addstar.h"
+# include "sedstar.h"
+
+int SEDfit (Catalog *outcat, Catalog *incat, SEDtable *table) {
+  
+  Graphdata graphdata;
+  KapaSection magSection, resSection;
+
+  SEDtableRow sourceValue, sourceError;
+  SEDfit minFit, testFit;
+
+  /* defaults */
+  outcat.average = NULL;
+  outcat.measure = NULL;
+  outcat.secfilt = NULL;
+
+  SEDtable = NULL;
+  SEDtableRaw = NULL;
+  sourceValue.mags = NULL;
+  sourceError.mags = NULL;
+  wavecode = NULL;
+  hashcode = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
+  magSection.name = NULL;
+  resSection.name = NULL;
+
+  oldsignal = signal (SIGINT, handle_interrupt);
+  interrupt = FALSE;
+
+  /* load photcode information */
+  if (!InitPhotcodes ()) goto escape;
+  Nsec = GetPhotcodeNsecfilt ();
+
+  /* interpret command-line options */
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 4)) goto usage;
+
+  PLOT = FALSE;
+  if ((N = get_argument (argc, argv, "-plot"))) {
+    remove_argument (N, &argc, argv);
+    PLOT = TRUE;
+  }
+
+  SAVEDIR = NULL;
+  if ((N = get_argument (argc, argv, "-save"))) {
+    remove_argument (N, &argc, argv);
+    SAVEDIR = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* interpret command-line options */
+  if (argc != 6) goto usage;
+
+  Nfit = 0;
+  colorP = GetPhotcodeCodebyName (argv[3]);
+  colorM = GetPhotcodeCodebyName (argv[5]);
+  if (!colorP || !colorM) goto color_undefined;
+
+  // artificially set USNOred and blu errors to 0.3
+  USNOred = GetPhotcodeCodebyName ("USNO_RED");
+  USNOblu = GetPhotcodeCodebyName ("USNO_BLUE");
+
+  // load SED table
+  f = fopen (argv[1], "r");
+  if (f == NULL) goto table_missing;
+
+  // XXX add error checks for header data
+  scan_line (f, line);
+  sscanf (line, "%*s %*s %d", &Nfilter);
+
+  // load SED table photcodes, generate the photcode hashtable
+  ALLOCATE (hashcode, int, 0x10000);
+  ALLOCATE (wavecode, float, Nfilter);
+  ALLOCATE (vegaToAB, float, Nfilter);
+
+  for (i = 0; i < 0x10000; i++) hashcode[i] = -1;
+  for (i = 0; i < Nfilter; i++) {
+    scan_line (f, line);
+    sscanf (line, "%*s %s %f %f", name, &wavecode[i], &vegaToAB[i]);
+    code = GetPhotcodeCodebyName (name);
+    if (code == 0) goto code_missing;
+    hashcode[code] = i;
+  }
+  codeP = hashcode[colorP];
+  codeM = hashcode[colorM];
+  if ((codeP == -1) || (codeM == -1)) goto color_missing;
+    
+  // skip remaining header lines
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  scan_line (f, line);
+  
+  // load the SED table data
+  Nrow = 0;
+  NROW = 100;
+  ALLOCATE (SEDtableRaw, SEDtableRow, NROW);
+  while (scan_line(f, line) != EOF) {
+    fparse (&SEDtableRaw[Nrow].Temp, 1, line);
+    fparse (&SEDtableRaw[Nrow].Av, 2, line);
+    ALLOCATE (SEDtableRaw[Nrow].mags, float, Nfilter);
+    for (i = 0; i < Nfilter; i++) {
+      fparse (&SEDtableRaw[Nrow].mags[i], i + 3, line);
+    }
+    SEDtableRaw[Nrow].color = SEDtableRaw[Nrow].mags[codeP] - SEDtableRaw[Nrow].mags[codeM];
+    Nrow ++;
+    CHECK_REALLOCATE (SEDtableRaw, SEDtableRow, NROW, Nrow, 100);
+  }      
+
+  // sort the SEDtable by the reference colors
+  SEDtable = sort_SEDtable (SEDtableRaw, Nrow);
+
+  // create holder for the source data
+  ALLOCATE (sourceValue.mags, float, Nfilter);
+  ALLOCATE (sourceError.mags, float, Nfilter);
+
+  if (PLOT) {
+    if (!GetGraph (&graphdata, &fd, NULL)) return (FALSE);
+    SetLimitsRaw (wavecode, NULL, Nfilter, &graphdata);
+    graphdata.style = 2;
+    graphdata.ptype = 2;
+    KapaClear (fd, TRUE);
+    magSection.name = strcreate ("mag");
+    magSection.x  = 0;
+    magSection.dx = 1;
+    magSection.y  = 0.5;
+    magSection.dy = 0.5;
+    resSection.name = strcreate ("res");
+    resSection.x  = 0.0;
+    resSection.dx = 1.0;
+    resSection.y  = 0.0;
+    resSection.dy = 0.5;
+    
+    KapaSetFont (fd, "helvetica", 14);
+    ALLOCATE (fitmags, float, Nfilter);
+    ALLOCATE (fiterrs, float, Nfilter);
+  }
+
+  /* load region corresponding to selection above */
+  if ((skylist = SelectRegions (RegionName, RegionList)) == NULL) goto escape;
+
+  /* loop over regions, extract data for each region */
+  // XXX add interrupt checks
+  fprintf (stderr, "using %d possible regions\n", skylist[0].Nregions);
+  for (k = 0; k < skylist[0].Nregions; k++) {
+    /* lock, load, unlock catalog */
+    catalog.filename = skylist[0].filename[k];
+    switch (lock_catalog (&catalog, LCK_SOFT)) {
+      case 2:
+	unlock_catalog (&catalog);
+      case 0:
+	catalog.Naverage = 0;
+	continue;
+    }
+    catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
+    if (!load_catalog (&catalog, TRUE)) {
+      catalog.Naverage = 0;
+      unlock_catalog (&catalog);
+      continue;
+    }
+
+    // perform the fit to all sources
+    for (i = 0; i < catalog.Naverage; i++) {
+
+      // blank out the source array
+      for (j = 0; j < Nfilter; j++) {
+	sourceValue.mags[j] = 100;
+      }	
+
+      // load the measurements for this source
+      m = catalog.average[i].offset;
+      for (j = 0; j < catalog.average[i].Nm; j++) {
+	idx = hashcode[catalog.measure[m+j].source];
+	if (idx == -1) continue;
+	// XXX do something more clever if more than one value exists per photcode
+	sourceValue.mags[idx] = catalog.measure[m+j].M_PS + vegaToAB[idx];
+	sourceError.mags[idx] = catalog.measure[m+j].dM_PS;
+	if ((catalog.measure[m+j].source == USNOred) || (catalog.measure[m+j].source == USNOblu)) {
+	  sourceError.mags[idx] = 0.3;
+	}
+      }
+
+      // XXX for the moment, skip sources without ref color
+      if (sourceValue.mags[codeP] > 50) continue;
+      if (sourceValue.mags[codeM] > 50) continue;
+      color = sourceValue.mags[codeP] - sourceValue.mags[codeM];
+
+      // XXX find tableRow within 0.1 mag of color
+      start = SEDcolorBracket (SEDtable, Nrow, color);
+      minFit = SEDchisq (SEDtable[start], &sourceValue, &sourceError, Nfilter);
+      minFit.row = start;
+
+      // search for min chisq backwards
+      done = FALSE;
+      row = start - 1;
+      while (!done && (row > 0)) {
+	testFit = SEDchisq (SEDtable[row], &sourceValue, &sourceError, Nfilter);
+	if (testFit.chisq < minFit.chisq) {
+	  minFit = testFit;
+	  minFit.row = row;
+	}
+	if (fabs(SEDtable[row][0].color - color) > 0.5) done = TRUE;
+	row --;
+      }
+
+      // search for min chisq forwards
+      done = FALSE;
+      row = start + 1;
+      while (!done && (row < Nrow)) {
+	testFit = SEDchisq (SEDtable[row], &sourceValue, &sourceError, Nfilter);
+	if (testFit.chisq < minFit.chisq) {
+	  minFit = testFit;
+	  minFit.row = row;
+	}
+	if (fabs(SEDtable[row][0].color - color) > 0.5) done = TRUE;
+	row ++;
+      }
+
+      Nfit ++;
+      // create the vectors for the example plots
+      if (PLOT) {
+	double tmp;
+	// find plot range
+	SetLimitsRaw (NULL, SEDtable[minFit.row][0].mags, Nfilter, &graphdata);
+	SWAP (graphdata.ymin, graphdata.ymax);
+
+	KapaClear (fd, TRUE);
+	KapaSetSection (fd, &magSection);
+    	KapaSetLimits (fd, &graphdata);
+	KapaBox (fd, &graphdata);
+	graphdata.color = KapaColorByName ("blue");
+	graphdata.etype = 0;
+	KapaPrepPlot (fd, Nfilter, &graphdata);
+	KapaPlotVector (fd, Nfilter, wavecode);
+	KapaPlotVector (fd, Nfilter, SEDtable[minFit.row][0].mags);
+	graphdata.color = KapaColorByName ("red");
+	graphdata.etype = 1;
+	for (j = 0; j < Nfilter; j++) {
+	  fitmags[j] = 100;
+	  fiterrs[j] = 0;
+	  if (sourceValue.mags[j] > 50) continue;
+	  fitmags[j] = sourceValue.mags[j] - minFit.Md;
+	  fiterrs[j] = sourceError.mags[j];
+	}
+	KapaPrepPlot (fd, Nfilter, &graphdata);
+	KapaPlotVector (fd, Nfilter, wavecode);
+	KapaPlotVector (fd, Nfilter, fitmags);
+	KapaPlotVector (fd, Nfilter, fiterrs);
+	KapaPlotVector (fd, Nfilter, fiterrs);
+	KapaSendLabel (fd, "model,fit (mags)", 1);
+
+	sprintf (line, "star: %10.6f %10.6f  T: %5.0fK  A_V|: %4.2f  M_D|: %5.2f  &sc&h^2|: %5.2f", 
+		 catalog.average[i].R, catalog.average[i].D, 
+		 SEDtable[minFit.row][0].Temp, SEDtable[minFit.row][0].Av, minFit.Md, minFit.chisq);
+	KapaSendLabel (fd, line, 2);
+	KapaSendLabel (fd, "model,fit (mags)", 1);
+
+	KapaSetSection (fd, &resSection);
+	graphdata.ymin = -1.0;
+	graphdata.ymax = +1.0;
+    	KapaSetLimits (fd, &graphdata);
+	KapaBox (fd, &graphdata);
+	graphdata.color = KapaColorByName ("red");
+	graphdata.etype = 1;
+
+	for (j = 0; j < Nfilter; j++) {
+	  fitmags[j] = 100;
+	  fiterrs[j] = 0;
+	  if (sourceValue.mags[j] > 50) continue;
+	  fitmags[j] = sourceValue.mags[j] - minFit.Md - SEDtable[minFit.row][0].mags[j];
+	  fiterrs[j] = sourceError.mags[j];
+	}
+	KapaPrepPlot (fd, Nfilter, &graphdata);
+	KapaPlotVector (fd, Nfilter, wavecode);
+	KapaPlotVector (fd, Nfilter, fitmags);
+	KapaPlotVector (fd, Nfilter, fiterrs);
+	KapaPlotVector (fd, Nfilter, fiterrs);
+	KapaSendLabel (fd, "wavelength (nm)", 0);
+	KapaSendLabel (fd, "resid (mags)", 1);
+
+	KiiCursorOn (fd);
+	while (KiiCursorRead (fd, &X, &Y, key)) {
+	  fprintf (stderr, "window: %f %f (%s)\n", X, Y, key);
+	  if (!strcasecmp (key, "Q")) {
+	    KiiCursorOff (fd);
+	    break;
+	  }
+	  if (!strcasecmp (key, "ESCAPE")) {
+	    KiiCursorOff (fd);
+	    unlock_catalog (&catalog);
+	    goto escape;
+	  }
+	}
+      }
+
+      // we now have the min chisq row. use this to supply the other filter values....
+    }
+    unlock_catalog (&catalog);
+  }
+  fprintf (stderr, "fitted %d stars\n", Nfit);
+  status = TRUE;
+  goto finish;
+  
+usage:
+  fprintf (stderr, "USAGE: fitset (sedtable) : (F) - (F)\n");
+  goto escape;
+
+table_missing:
+  fprintf (stderr, "ERROR: can't open the SED table\n");
+  goto escape;
+
+color_missing:
+  fprintf (stderr, "ERROR: reference color not in SED table\n");
+  goto escape;
+
+color_undefined:
+  fprintf (stderr, "ERROR: undefined photcode in reference color\n");
+  goto escape;
+
+code_missing:
+  fprintf (stderr, "ERROR: undefined photcode in SED table\n");
+  goto escape;
+
+escape:
+  status = FALSE;
+  goto finish;
+
+finish:
+  if (skylist != NULL) SkyListFree (skylist, ((RegionName != NULL) || (RegionList != NULL)));
+  if (catalog.average != NULL) free (catalog.average);
+  if (catalog.secfilt != NULL) free (catalog.secfilt);
+  if (catalog.measure != NULL) free (catalog.measure);
+
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
+  if (wavecode != NULL) free (wavecode);
+  if (hashcode != NULL) free (hashcode);
+  if (SEDtableRaw != NULL) {
+    for (i = 0; i < Nrow; i++) {
+      free (SEDtableRaw[i].mags);
+    }
+    free (SEDtableRaw);
+  }
+  if (SEDtable != NULL) free (SEDtable);
+  if (sourceValue.mags != NULL) free (sourceValue.mags);
+  if (sourceError.mags != NULL) free (sourceError.mags);
+
+  signal (SIGINT, oldsignal);
+  return (status);
+}
+
Index: /trunk/Ohana/src/addstar/src/SEDops.c
===================================================================
--- /trunk/Ohana/src/addstar/src/SEDops.c	(revision 7692)
+++ /trunk/Ohana/src/addstar/src/SEDops.c	(revision 7693)
@@ -1,2 +1,4 @@
+# include "addstar.h"
+# include "sedstar.h"
 
 // fit the data (with errors) to the given table row
Index: /trunk/Ohana/src/addstar/src/SEDtableLoad.c
===================================================================
--- /trunk/Ohana/src/addstar/src/SEDtableLoad.c	(revision 7692)
+++ /trunk/Ohana/src/addstar/src/SEDtableLoad.c	(revision 7693)
@@ -2,7 +2,14 @@
 # include "sedstar.h"
 
-SEDtable *SEDtableLoad (char *filename, int *nrow) {
+// XXX where do the colorP and colorM codes come from?
+SEDtable *SEDtableLoad (char *filename) {
 
   FILE *f;
+  char line[1024], name[64];
+  char colornameP[64], colornameM[64];
+  int colorP, colorM;
+  int i, Nrow, NROW;
+
+  ALLOCATE (table, SEDtable, 1);
 
   // load SED table
@@ -12,22 +19,28 @@
   // XXX add error checks for header data
   scan_line (f, line);
-  sscanf (line, "%*s %*s %d", &Nfilter);
+  sscanf (line, "%*s %*s %d", &table[0].Nfilter);
 
   // load SED table photcodes, generate the photcode hashtable
-  ALLOCATE (hashcode, int, 0x10000);
-  ALLOCATE (wavecode, float, Nfilter);
-  ALLOCATE (vegaToAB, float, Nfilter);
+  ALLOCATE (table[0].hashcode, int, 0x10000);
+  ALLOCATE (table[0].wavecode, float, table[0].Nfilter);
+  ALLOCATE (table[0].vegaToAB, float, table[0].Nfilter);
 
-  for (i = 0; i < 0x10000; i++) hashcode[i] = -1;
-  for (i = 0; i < Nfilter; i++) {
+  for (i = 0; i < 0x10000; i++) table[0].hashcode[i] = -1;
+  for (i = 0; i < table[0].Nfilter; i++) {
     scan_line (f, line);
-    sscanf (line, "%*s %s %f %f", name, &wavecode[i], &vegaToAB[i]);
+    sscanf (line, "%*s %s %f %f", name, &table[0].wavecode[i], &table[0].vegaToAB[i]);
     code = GetPhotcodeCodebyName (name);
     if (code == 0) goto code_missing;
-    hashcode[code] = i;
+    table[0].hashcode[code] = i;
   }
-  codeP = hashcode[colorP];
-  codeM = hashcode[colorM];
-  if ((codeP == -1) || (codeM == -1)) goto color_missing;
+
+  // define the selection color codes
+  sscanf (line, "# %s - %s", colornameP, colornameM);
+  colorP = GetPhotcodeCodebyName (colornameP);
+  colorM = GetPhotcodeCodebyName (colornameM);
+  table[0].codeP = table[0].hashcode[colorP];
+  table[0].codeM = table[0].hashcode[colorM];
+  if (table[0].codeP == -1) Shutdown ("missing positive color filter");
+  if (table[0].codeM == -1) Shutdown ("missing positive color filter");
     
   // skip remaining header lines
@@ -37,5 +50,5 @@
   scan_line (f, line);
   
-  // load the SED table data
+  // load the SED raw table rows
   Nrow = 0;
   NROW = 100;
@@ -44,6 +57,6 @@
     fparse (&SEDtableRaw[Nrow].Temp, 1, line);
     fparse (&SEDtableRaw[Nrow].Av, 2, line);
-    ALLOCATE (SEDtableRaw[Nrow].mags, float, Nfilter);
-    for (i = 0; i < Nfilter; i++) {
+    ALLOCATE (SEDtableRaw[Nrow].mags, float, table[0].Nfilter);
+    for (i = 0; i < table[0].Nfilter; i++) {
       fparse (&SEDtableRaw[Nrow].mags[i], i + 3, line);
     }
@@ -54,3 +67,6 @@
 
   // sort the SEDtable by the reference colors
-  SEDtable = sort_SEDtable (SEDtableRaw, Nrow);
+  table[0].row = sort_SEDtable (SEDtableRaw, Nrow);
+  table[0].Nrow = Nrow;
+  return (table);
+}
Index: /trunk/Ohana/src/addstar/src/args_sedstar.c
===================================================================
--- /trunk/Ohana/src/addstar/src/args_sedstar.c	(revision 7692)
+++ /trunk/Ohana/src/addstar/src/args_sedstar.c	(revision 7693)
@@ -34,38 +34,5 @@
     remove_argument (N, &argc, argv);
   }
-  /* override any header PHOTCODE values */
-  options.photcode = 0;
-  if ((N = get_argument (argc, argv, "-p"))) {
-    remove_argument (N, &argc, argv);
-    options.photcode = GetPhotcodeCodebyName (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
 
-  options.timeref = 0; 
-  options.mosaic = FALSE;
-  options.existing_regions = FALSE;
-  options.skip_missed = FALSE;
-  options.closest = FALSE;
-
-  /* only add to existing objects */
-  options.only_match = FALSE;
-  if ((N = get_argument (argc, argv, "-only-match"))) {
-    options.only_match = TRUE;
-    remove_argument (N, &argc, argv);
-  }
-  /* replace measurement, don't duplicate (ref/cat only) */
-  options.replace = FALSE;
-  if ((N = get_argument (argc, argv, "-replace"))) {
-    options.replace = TRUE;
-    remove_argument (N, &argc, argv);
-  }
-
-  /* select quality flags */
-  SELECT_2MASS_QUALITY = NULL;
-  if ((N = get_argument (argc, argv, "-2massquality"))) {
-    remove_argument (N, &argc, argv);
-    SELECT_2MASS_QUALITY = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
   /* extra error messages */
   VERBOSE = FALSE;
@@ -76,4 +43,11 @@
 
   /* other defaults */
+  options.timeref = 0; 
+  options.mosaic = FALSE;
+  options.existing_regions = FALSE;
+  options.skip_missed = FALSE;
+  options.closest = FALSE;
+  options.only_match = FALSE;
+  options.replace = FALSE;
   options.nosort = FALSE;
   options.update = FALSE;
@@ -87,6 +61,6 @@
   DUMP = NULL;
 
-  if (argc != 1) {
-    fprintf (stderr, "USAGE: load2mass\n");
+  if (argc != 3) {
+    fprintf (stderr, "USAGE: sedstar (sedtable) (outcatalog) [-region Rmin Rmax Dmin Dmax]\n");
     exit (2);
   }
@@ -97,11 +71,9 @@
 
   fprintf (stderr, "USAGE\n");
-  fprintf (stderr, "  load2mass");
-  fprintf (stderr, "     add data from 2MASS catalog to fullsky\n\n");
+  fprintf (stderr, "  sedstar (SEDtable)");
+  fprintf (stderr, "     fit objects to stellar SEDs\n\n");
 
   fprintf (stderr, "  optional flags:\n");
-  fprintf (stderr, "  -p (photcode)         	  : specify photcode (override header)\n");
-  fprintf (stderr, "  -only-match           	  : only add measurements to existing objects\n");
-  fprintf (stderr, "  -replace              	  : replace time/photcode measurements (no duplication)\n");
+  fprintf (stderr, "  -region Rmin Rmax Dmin Dmax : sky region for analysis\n");
   fprintf (stderr, "  -v                    	  : verbose mode\n");
   fprintf (stderr, "  -help                 	  : this list\n");
Index: /trunk/Ohana/src/addstar/src/sedstar.c
===================================================================
--- /trunk/Ohana/src/addstar/src/sedstar.c	(revision 7692)
+++ /trunk/Ohana/src/addstar/src/sedstar.c	(revision 7693)
@@ -7,4 +7,5 @@
   SkyTable *sky, *sky2mass;
   AddstarClientOptions options;
+  Catalog incatalog, outcatalog;
 
   // need to construct these options with args_load2mass...
@@ -16,8 +17,31 @@
   
   // select regions of interest
+  skylist = SkyListByPatch (sky, -1, &UserPatch);
 
   // load the SED data table
-  sedtable = SEDtableLoad (SEDfile, &Nsedtable);
+  sedtable = SEDtableLoad (argv[1]);
 
+  for (i = 0; i < skylist[0].Nregions; i++) {
+    incatalog.filename = skylist[0].filename[i];
+    load_pt_catalog (&incatalog, skylist[0].regions[i]);
+
+    // create output catalog filename
+    root = strstr (incatalog.filename, CATDIR);
+    if (root == NULL) Shutdown ("error with input catalog name");
+
+    extension = incatalog.filename + strlen(CATDIR);
+    ALLOCATE (outcatalog.filename, char, strlen(argv[2]) + strlen(extension) + 2);
+    sprintf (outcatalog.filename, "%s/%s", argv[2], extension);
+
+    status = lock_catalog (outcatalog, LCK_XCLD);
+    if (status != 2) Shutdown ("error: output catalog exists");
+    mkcatalog (region, catalog); /* fills in new header info */
+
+    SEDfit (&outcatalog, &incatalog, sedtable);
+    
+    wcatalog (outcatalog);
+    unlock_catalog (outcatalog);
+    free (outcatalog.filename);
+  }
 
   exit (0);