Index: branches/ipp-magic-v0/Ohana/src/imregister/base/get_fwhm.c
===================================================================
--- trunk/Ohana/src/imregister/base/get_fwhm.c	(revision 21045)
+++ 	(revision )
@@ -1,130 +1,0 @@
-# include "imregister.h"
-
-/* see note at end of file */
-double get_fwhm (char *filename) {
-
-  FILE *f;
-  char line[1024];
-  double *fwhm, *mag, *Mag, *FWHM, *Nfwhm, flags;
-  double F, value, Mlim, Vmode, Fmin, Fmax;
-  int i, N, n, Nwant, bin, mode, Nmode;
-
-  /* get sextract file (has fwhm info) */
-  f = fopen (filename, "r");
-  if (f == (FILE *) NULL) {
-    fprintf (stderr, "ERROR: can't open file.sdat\n");
-    return (-1.0);
-  }
-
-  n = 0;
-  N = 1000;
-  ALLOCATE (fwhm, double, N);
-  ALLOCATE (mag, double, N);
-  
-  /* load in data from file */
-  while (scan_line (f, line) != EOF) {
-    dparse (&flags, 2, line);
-    if (flags > 0) continue;
-    dparse (&mag[n], 4, line);
-    if (mag[n] > 0) continue;
-    dparse (&fwhm[n], 1, line);
-    n++;
-    if (n == N) {
-      N += 1000;
-      REALLOCATE (fwhm, double, N);
-      REALLOCATE (mag, double, N);
-    }
-  }    
-  if (n == 0) { return (-1.0); }
-
-  /* find the mag threshold */
-  ALLOCATE (Mag, double, n);
-  memcpy (Mag, mag, n*sizeof(double));
-  dsort (Mag, n);
-  Nwant = MIN (MAX (0.25*n, 10), n-1);
-  Mlim = Mag[Nwant];
-
-  /* create subset of fwhm with mag < Mlim */
-  ALLOCATE (FWHM, double, n);
-  ALLOCATE (Nfwhm, double, 1024);
-  bzero (FWHM, n*sizeof(double));
-  bzero (Nfwhm, 1024*sizeof(double));
-  
-  N = 0;
-  for (i = 0; i < n; i++) {
-    if (mag[i] > Mlim) continue;
-    FWHM[N] = fwhm[i];
-    N++;
-
-    /* accumulate mag histogram */
-    bin = fwhm[i]*10.0;
-    if (bin < 0) continue;
-    if (bin > 999) continue;
-    Nfwhm[bin] ++;
-  }
-
-  /* find the mode (bin size = 0.1 pixels) */
-  mode = 0; 
-  Nmode = Nfwhm[mode];
-  for (i = 0; i < 1000; i++) {
-    if (Nfwhm[i] > Nmode) {
-      mode = i;
-      Nmode = Nfwhm[i];
-    }
-  }
-  Vmode = mode * 0.1;
-  Fmin = Vmode - 0.2;
-  Fmax = Vmode + 0.2;
-
-  /* average bins within 0.2 of mode */
-  F = 0;
-  n = 0;
-  for (i = 0; i < N; i++) {
-    if (FWHM[i] < Fmin) continue;
-    if (FWHM[i] > Fmax) continue;
-    n++;
-    F+=FWHM[i];
-  }
-  if (n > 3) {
-    value = F / n;
-    return (value);
-  }
-
-  /* if we don't get the answer from the above, resort to median */
-  dsort (FWHM, N);
-  F = 0;
-  n = 0;
-  for (i = 0.25*N; i < 0.75*N; i++) {
-    n++;
-    F+=FWHM[i];
-  }
-  if (n > 0) {
-    value = F / n;
-  } else {
-    value = -1;
-  }
-  return (value);
-
-}
-
-/* 
-
-   calculate the fwhm from the data in the given file.
-   The file contains output from sextractor (or other program, perhaps).
-   The columns that matter are: 
-     1) fwhm (pixels)
-     2) flags
-     4) mag
-
-   this function determines the magnitude range for the file,
-   selects stars in the upper 25% of that magnitude range (by number),
-   rejecting bad stars (flags > 0) along the way
-   
-   the resulting list of fwhms are examined for the mode, in 0.1 pixel bins.
-   stars with fwhm within 0.2 pix of that mode are averaged to find the
-   best fwhm value.
-
-   in the database, a value of 0 means the measurement has not been tried.
-   a value of -1 means the measurement failed
-
-*/
Index: branches/ipp-magic-v0/Ohana/src/imregister/detrend/usage.c
===================================================================
--- trunk/Ohana/src/imregister/detrend/usage.c	(revision 21045)
+++ 	(revision )
@@ -1,30 +1,0 @@
-# include "imregister.h"
-# include "detrend.h"
-
-int usage () {
-  
-  fprintf (stderr, "detsearch : select images from the detrend database\n");
-  fprintf (stderr, " USAGE: detsearch [-options...]\n");
-  fprintf (stderr, " -h : show this list\n");
-  fprintf (stderr, " --help : show this list\n");
-  fprintf (stderr, " -image filename (ccd) (mode) : find matching detrend data for this image \n");
-  fprintf (stderr, " -type (type) : limit selection by image type (bias, dark, flat, etc)\n");
-  fprintf (stderr, " -time yyyy/mm/dd,HH:MM:SS : specify time for selection\n");
-  fprintf (stderr, " -ccd (N) : limit selection to this ccd\n");
-  fprintf (stderr, " -filter (name) : limit selection to this filter\n");
-  fprintf (stderr, " -exptime (value) : limit selection to match this exposure time\n");
-  fprintf (stderr, " -tstop : display end of valid time range\n");
-  fprintf (stderr, " -treg  : display time of image registration\n");
-  fprintf (stderr, " -ve    : Elixir verbose (SUCCESS / ERROR) \n");
-  fprintf (stderr, " -quiet : Elixir quiet (no SUCCESS / ERROR) \n");
-  fprintf (stderr, " -close : select detrend data which is closest in time, overlap not forced\n");
-  fprintf (stderr, " -entry (value) : limit selection to match this entry\n");
-  fprintf (stderr, " -match (value) : only list the Nth matched entry\n");
-  fprintf (stderr, " -label (word)  : limit selection to match this label\n");
-  fprintf (stderr, " -select : force selection of the 'best' match\n");
-  fprintf (stderr, " -del    : delete the matched entries\n");
-  fprintf (stderr, " -delete : delete the matched entries\n");
-  fprintf (stderr, " -modify (entry) (value) : change entry in the selection to value\n");
-  fprintf (stderr, "   possible -modify entries: label, order, tstart, tstop\n");
-  exit (2);
-}
Index: branches/ipp-magic-v0/Ohana/src/imregister/imreg/load_probes.c
===================================================================
--- trunk/Ohana/src/imregister/imreg/load_probes.c	(revision 21045)
+++ 	(revision )
@@ -1,108 +1,0 @@
-# include "imregister.h"
-# define NBYTES 1024
-# define WINDOW 1800
-/* probes are updated every 600 seconds.  give us some leeway here */
-
-int load_probes (char *filename, unsigned long tzero, int *wantprobe, double *values, int Nprobe) {
-
-  FILE *f;
-  char line[256], closeline[256], *buffer, *c, *p;
-  int nprobe, probe[10];
-  int i, Nstart, Nread, Nbytes, Nshift, Nleft, done, close_enough;
-  double time, jdstart, tmp;
-  int Nsec, sec;
-
-  f = fopen (filename, "r");
-  if (f == (FILE *) NULL) {
-    fprintf (stderr, "can't open DataLogger file: %s\n", filename);
-    return (FALSE);
-  }
-
-  /* assume fixed format */
-  scan_line (f, line);
-  scan_line (f, line);
-  scan_line (f, line);
-  sscanf (line, "%*s %*s %*s %*s %lf", &jdstart);
-  Nsec = (jdstart - 2440587.5)*86400;
-  if (tzero < Nsec - WINDOW) {
-    fprintf (stderr, "missing time in DataLogger file\n");
-    return (FALSE);
-  }
-
-  scan_line (f, line);
-  scan_line (f, line);
-  scan_line (f, line);
-  sscanf (line, "%*s %*s %d", &nprobe);
-  for (i = 0; (i < nprobe) && (i < 10); i++) {
-    dparse (&tmp, i+4, line);
-    probe[i] = tmp;
-  }
-
-  /* check that probes match desired probes -- demand same order in file */
-  if (nprobe != Nprobe) {
-    fprintf (stderr, "wrong number of probes\n");
-    return (FALSE);
-  }
-  for (i = 0; (i < nprobe) && (i < 10); i++) {
-    if (probe[i] != wantprobe[i]) {
-      fprintf (stderr, "probe mismatch\n");
-      return (FALSE);
-    }
-  }
-  
-  ALLOCATE (buffer, char, NBYTES+1);
-    
-  Nstart = 0;
-  done = FALSE;
-  close_enough = FALSE;
-  while (!done && ((Nread = fread (&buffer[Nstart], 1, NBYTES - Nstart, f)) > 0)) {
-    Nbytes = Nread + Nstart;
-    buffer[Nbytes] = 0;
-    
-    /* we are using strchr lib functions - buffer can't have NULLs */
-    c = strrchr (buffer, '\n');
-    if (c == (char *) NULL) c = &buffer[Nbytes-1];
-    Nshift = c - buffer + 1;
-    
-    /* limit searches to the range buffer[0] - buffer[Nshift] */
-    p = &buffer[0];
-    while (1) {
-      c = strchr (p, '\n');
-      if (c == (char *) NULL) break;
-      if (c > &buffer[Nshift]) break;
-      *c = 0;
-      dparse (&time, 0, p);
-      sec = Nsec + time*86400;
-      if (sec + WINDOW > tzero) {
-	strcpy (closeline, p);
-	close_enough = TRUE;
-      }
-      if (sec > tzero) {
-	for (i = 0; i < Nprobe; i++) {
-	  dparse (&values[i], i+5, p);
-	}
-	free (buffer);
-	return (TRUE);
-      }
-      *c = '\n';
-      p = c + 1;
-    }
-    
-    Nleft  = Nbytes - Nshift;
-    if (Nleft > 0) memmove (buffer, &buffer[Nshift], Nleft);
-    
-    Nstart = Nleft;
-  }
-
-  if (close_enough) {
-    for (i = 0; i < Nprobe; i++) {
-      dparse (&values[i], i+5, closeline);
-    }
-    free (buffer);
-    return (TRUE);
-  }
-
-  fprintf (stderr, "time not in DataLogger file %ld\n", tzero);
-  free (buffer);
-  return (FALSE);
-}
Index: branches/ipp-magic-v0/Ohana/src/opihi/cmd.astro/fixcols.c
===================================================================
--- branches/ipp-magic-v0/Ohana/src/opihi/cmd.astro/fixcols.c	(revision 21062)
+++ branches/ipp-magic-v0/Ohana/src/opihi/cmd.astro/fixcols.c	(revision 21062)
@@ -0,0 +1,246 @@
+# include "astro.h"
+
+// this function is a hack to mask bad columns / rows in an image.  we assume that there
+// are columns or rows which are, as a whole, significant deviant from the others.
+// in fact, we measure this by segments
+
+// SIGMA is a robust measure of the sigma, STDEV is the formal stdev of the values in range
+int vector_stats (float *vect, int Nvect, float *MEDIAN, float *SIGMA, float *STDEV);
+
+int fixrows (int argc, char **argv) {
+  
+  int ix, iy, Nx, Ny, Nvect, start, stop;
+  float *Vin, *vect, *meds, median, stdev, sigma, Nsigma, value;
+  Buffer *in;
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: fixrows (buffer) (col_start) (col_stop) (Nsigma) (value)\n");
+    return (FALSE);
+  }
+
+  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  start = atoi (argv[2]);
+  stop = atoi (argv[3]);
+  Nsigma = atof (argv[4]);
+  value = atof (argv[5]);
+
+  Nx = in[0].matrix.Naxis[0];
+  Ny = in[0].matrix.Naxis[1];
+
+  ALLOCATE (meds, float, Ny);
+  ALLOCATE (vect, float, stop - start);
+
+  Vin  = (float *)in[0].matrix.buffer;
+
+  // measure the median for each column segment, accumulate the stats
+  for (iy = 0; iy < Ny; iy++) {
+
+    Nvect = 0;
+    for (ix = start; ix < stop; ix++) {
+      vect[Nvect] = Vin[ix + iy*Nx];
+      Nvect ++;
+    }
+
+    // get the median, robust sigma, and formal stdev of the column segment
+    vector_stats (vect, Nvect, &median, &sigma, &stdev);
+    meds[iy] = median;
+    // fprintf (stderr, "median : %f +/- %f (%f)\n", median, sigma, stdev);
+  }
+
+  // get the median, robust sigma, and formal stdev of the medians
+  vector_stats (meds, Nx, &median, &sigma, &stdev);
+  // fprintf (stderr, "median of medians : %f +/- %f (%f)\n", median, sigma, stdev);
+
+  // mask values which are more than Nsigma out of range
+  for (iy = 0; iy < Ny; iy++) {
+    if (fabs(meds[iy] - median) < Nsigma * sigma) continue;
+    for (ix = start; ix < stop; ix++) {
+      Vin[ix + iy*Nx] = value;
+    }
+  }
+  
+  free (meds);
+  free (vect);
+
+  return (TRUE);
+}
+
+int fixcols (int argc, char **argv) {
+  
+  int ix, iy, Nx, Ny, Nvect, start, stop;
+  float *Vin, *vect, *meds, median, stdev, sigma, Nsigma, value;
+  Buffer *in;
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: fixcols (buffer) (row_start) (row_stop) (Nsigma) (value)\n");
+    return (FALSE);
+  }
+
+  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  start = atoi (argv[2]);
+  stop = atoi (argv[3]);
+  Nsigma = atof (argv[4]);
+  value = atof (argv[5]);
+
+  Nx = in[0].matrix.Naxis[0];
+  Ny = in[0].matrix.Naxis[1];
+
+  ALLOCATE (meds, float, Nx);
+  ALLOCATE (vect, float, stop - start);
+
+  Vin  = (float *)in[0].matrix.buffer;
+
+  // measure the median for each column segment, accumulate the stats
+  for (ix = 0; ix < Nx; ix++) {
+
+    Nvect = 0;
+    for (iy = start; iy < stop; iy++) {
+      vect[Nvect] = Vin[ix + iy*Nx];
+      Nvect ++;
+    }
+
+    // get the median, robust sigma, and formal stdev of the column segment
+    vector_stats (vect, Nvect, &median, &sigma, &stdev);
+    meds[ix] = median;
+    // fprintf (stderr, "median : %f +/- %f (%f)\n", median, sigma, stdev);
+  }
+
+  // get the median, robust sigma, and formal stdev of the medians
+  vector_stats (meds, Nx, &median, &sigma, &stdev);
+  // fprintf (stderr, "median of medians : %f +/- %f (%f)\n", median, sigma, stdev);
+
+  // mask values which are more than Nsigma out of range
+  for (ix = 0; ix < Nx; ix++) {
+    if (fabs(meds[ix] - median) < Nsigma * sigma) continue;
+    for (iy = start; iy < stop; iy++) {
+      Vin[ix + iy*Nx] = value;
+    }
+  }
+  
+  free (meds);
+  free (vect);
+
+  return (TRUE);
+}
+
+// this macro is accumulating the histogram until the sum passes the desired threshold
+// point the fractional position of that threshold is found by interpolation
+# define SUMPOINT(POINT, NPOINT, FRACTION)				\
+  if ((NPOINT) == -1) {							\
+    (POINT) += Nval[i];							\
+    if ((POINT) >= FRACTION*N) {					\
+      (NPOINT) = i;							\
+      value = i * dx + min;						\
+      if (Nval[i] != 0) {						\
+	/* POINT for the bin center is not exactly FRACTION*N */	\
+	/* interpolate to the true value, but limit the scale of the */	\
+	/* interpolation to a single bin */				\
+	adjust = dx * (FRACTION*N - (POINT)) / Nval[i];			\
+	value += MAX (MIN (adjust, dx), -dx);				\
+      }									\
+      (POINT) = value;							\
+    }									\
+  }
+
+// for testing, this is a useful line:
+// fprintf (stderr, "for point %f (%f): bin: %d = %f, dx: %f, adjust: %f, value: %f, POINT: %f, Nval[i]: %d, Nval[i-1]: %d\n",
+// FRACTION, FRACTION*N, i, i*dx + min, dx, adjust, value, POINT, Nval[i], Nval[i-1]);
+
+int vector_stats (float *vect, int Nvect, float *MEDIAN, float *SIGMA, float *STDEV) {
+  
+  int i, iter, N, done;
+  double max, min, sum, var, dvar, mean, stdev, sigma, adjust, value;
+  int Nval[1001], bin, Nmode, N_MED, N_LQ, N_UQ, N_SM, N_SP;
+  double dx, mode, MED, UQ, LQ, SM, SP;
+  float *X;
+
+  /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */
+
+  /* calculate max, min, mean, sum, npix */
+  max = -HUGE_VAL;
+  min = HUGE_VAL;
+  sum = N = 0;
+  X = vect;
+  for (i = 0; i < Nvect; i++, X++) {
+    if (!finite (*X)) continue;
+    max = MAX (*X, max);
+    min = MIN (*X, min);
+    sum += *X;
+    N++;
+  }      
+  mean = sum / N;
+
+  // calculate median and mode with initial resolution dx = (max - min) / 1000 
+  // we adjust the resolution and retry if the resulting robust sigma < dx
+  dx = (max - min) / 1000;
+  if (dx == 0) {
+    MED = mode = min;
+    stdev = 0.0;
+    sigma = 0.0;
+    goto skip;
+  }
+
+  done = FALSE;
+  for (iter = 0; (iter < 3) && !done; iter ++) {
+
+    bzero (Nval, 1001*sizeof(int));
+    var = 0;
+    X = vect;
+    for (i = 0; i < Nvect; i++, X++) {
+      if (!finite (*X)) continue;
+      bin = MAX (0, MIN (1000, (*X - min) / dx));
+      Nval[bin] ++;
+      if (bin == 0) continue;
+      if (bin == 1000) continue;
+      // only count the values in range for the stdev
+      dvar = (*X - mean);
+      var += dvar*dvar;
+    }      
+    stdev = sqrt (var / N);
+
+    Nmode = 0;
+    mode = Nval[Nmode];
+    MED = SM = SP = LQ = UQ = 0;
+    N_MED = N_SM = N_SP = N_LQ = N_UQ = -1;
+    for (i = 0; i < 1001; i++) {
+      SUMPOINT (MED, N_MED, 0.5000);
+      SUMPOINT (LQ,  N_LQ,  0.2500);
+      SUMPOINT (UQ,  N_UQ,  0.7500);
+      SUMPOINT (SM,  N_SM,  0.1587);
+      SUMPOINT (SP,  N_SP,  0.8413);
+
+      if (mode < Nval[i]) {
+	Nmode = i;
+	mode = Nval[Nmode];
+      }
+    }
+
+    sigma = 0.5*(SP - SM);
+    mode = Nmode * dx + min;
+
+    // if (sigma < dx), we adjust the range to be centered on the median
+    if (sigma < dx) {
+      dx *= 0.1;
+      min = MED - 500*dx;
+    } else {
+      done = TRUE;
+    }
+  }
+
+skip:
+
+  *MEDIAN = MED;
+  *SIGMA = sigma;
+  *STDEV = stdev;
+
+  // set_variable ("MIN",      min);
+  // set_variable ("MAX",      max);
+  // set_variable ("MEAN",     mean);
+  // set_variable ("MODE",     mode);
+  // set_variable ("TOTAL",    sum);
+  // set_int_variable ("NPIX", N);
+  // set_variable ("SIGMA",    stdev);
+
+  return (TRUE);
+}
+