Changeset 41341 for trunk/Ohana/src/opihi

trunk/Ohana/src/opihi/cmd.astro/cdensify.c

-              r41340
+              r41341
+  }
+  int binning = 1;
+  if ((N = get_argument (argc, argv, "-binning"))) {
+    remove_argument (N, &argc, argv);
+    binning = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
   int PSFTYPE = IS_DOT;
   if ((N = get_argument (argc, argv, "-psf"))) {
 …
   Ymax = graphmode.ymax;
   Ymin = graphmode.ymin;
+  dX = (Xmax - Xmin) / (Xpix - 1);
+  dY = (Ymax - Ymin) / (Ypix - 1);
+  // (dX, dY) are the pixel scale of the output image (output pixels / input pixels)
+  dX = binning * (Xmax - Xmin) / (Xpix - 1);
+  dY = binning * (Ymax - Ymin) / (Ypix - 1);
   CHECKVAL(Xmin);

trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c

-              r41340
+              r41341
   int kapa;
   Graphdata graphmode;
+  int binning = 1;
+  if ((N = get_argument (argc, argv, "-binning"))) {
+    remove_argument (N, &argc, argv);
+    binning = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
   range = NULL;
 …
   Ymax = graphmode.ymax;
   Ymin = graphmode.ymin;
+  dX = (Xmax - Xmin) / (Xpix - 1);
+  dY = (Ymax - Ymin) / (Ypix - 1);
+  // (dX, dY) are the pixel scale of the output image (output pixels / input pixels)
+  dX = binning * (Xmax - Xmin) / (Xpix - 1);
+  dY = binning * (Ymax - Ymin) / (Ypix - 1);
   CHECKVAL(Xmin);
 …
   CHECKVAL(dY);
   int Nx = (Xmax - Xmin) / dX + 1;
   int Ny = (Ymax - Ymin) / dY + 1;
+  int Nx = abs((Xmax - Xmin) / dX) + 1;
+  int Ny = abs((Ymax - Ymin) / dY) + 1;
   Coords newcoords = graphmode.coords;

trunk/Ohana/src/opihi/cmd.astro/cgrid.c

-              r41340
+              r41341
   Graphdata graphmode;
+  if ((N = get_argument (argc, argv, "-h"))) goto usage;
+  if ((N = get_argument (argc, argv, "--help"))) goto usage;
   RAbyHour = FALSE;
   if ((N = get_argument (argc, argv, "-ra-by-hour"))) {
 …
+  }
   int JustifyRA = 5;
   double JustifyDEC = 5;
+  int JustifyRA = 8;
+  int JustifyDEC = 2;
   if ((N = get_argument (argc, argv, "-justify-ra"))) {
     remove_argument (N, &argc, argv);
 …
   if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;
+  if (argc != 1) {
+    gprint (GP_ERR, "USAGE: cgrid [style] [-ra-by-hour] [-labels]\n");
+    return (FALSE);
+  }
+  if (argc != 1) goto usage;
   /* are we plotting one of the poles? */
 …
+    }
+  }
   /* set spacings for DEC */
 …
   return (TRUE);
+ usage:
+  gprint (GP_ERR, "USAGE: cgrid [style] [options]\n");
+  gprint (GP_ERR, "  options:\n");
+  gprint (GP_ERR, "  -h, --help: show this list\n");
+  gprint (GP_ERR, "  -ra-by-hour : RA grid lines will be space on rounded hour lines (default is degrees)\n");
+  gprint (GP_ERR, "  -labels : add RA & dec coordinates\n");
+  gprint (GP_ERR, "  -major-spacing : grid lines drawn at major tickmarks, not minor tickmarks\n");
+  gprint (GP_ERR, "  -ra-spacing : specify size of RA grid steps in degrees (ignores -ra-by-hour)\n");
+  gprint (GP_ERR, "  -dec-spacing : specify size of dec grid steps in degrees\n");
+  gprint (GP_ERR, "  -justify-ra : choose how RA labels are justified (see below)\n");
+  gprint (GP_ERR, "  -justify-dec : choose how dec labels are justified (see below)\n");
+  gprint (GP_ERR, "  -label-ra : RA coordinate of the Dec-line labels\n");
+  gprint (GP_ERR, "  -label-dec : Dec coordinate of the RA-line labels\n");
+  gprint (GP_ERR, "  -label-color : color for the labels (independent of grid lines)\n");
+  gprint (GP_ERR, "  text justification: text is justified horizontally and vertically based on the following numbers:\n");
+  gprint (GP_ERR, "   6 7 8\n");
+  gprint (GP_ERR, "   3 4 5\n");
+  gprint (GP_ERR, "   0 1 2\n");
+  return (FALSE);
+}

trunk/Ohana/src/opihi/cmd.astro/csystem.c

-              r41340
+              r41341
   Vector *uxvec = NULL;
   Vector *uyvec = NULL;
+  int Quiet = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
   // transform proper motions at the same time
 …
     if (uXname) {
+      gprint (GP_LOG, "%10.6f %10.6f @ %10.6f %10.6f\n", x, y, ux, uy);
+      if (!Quiet) {
+        gprint (GP_LOG, "%10.6f %10.6f @ %10.6f %10.6f\n", x, y, ux, uy);
+        switch (output) {
+          case COORD_CELESTIAL:
+            set_variable ("uR", ux);
+            set_variable ("uD", uy);
+            break;
+          case COORD_ECLIPTIC:
+            set_variable ("uL", ux);
+            set_variable ("uB", uy);
+            break;
+          case COORD_GALACTIC:
+          case COORD_GALACTIC_REID_2004:
+            set_variable ("uLon", ux);
+            set_variable ("uLat", uy);
+            break;
+          default:
+            break;
+        }
+      }
     } else {
+      gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
+      if (!Quiet) {
+        gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
+      }
+    }
+    switch (output) {
+      case COORD_CELESTIAL:
+        set_variable ("RA", x);
+        set_variable ("DEC", y);
+        break;
+      case COORD_ECLIPTIC:
+        set_variable ("Lambda", x);
+        set_variable ("Beta", y);
+        break;
+      case COORD_GALACTIC:
+      case COORD_GALACTIC_REID_2004:
+        set_variable ("gLon", x);
+        set_variable ("gLat", y);
+        break;
+      default:
+        break;
+    }
     free (transform);
 …
   gprint (GP_ERR, " e.g.   : csystem C G R D -pm uR uD -> R D will contain glon, glat; uR, uD will contain uL, uB\n");
   gprint (GP_ERR, " e.g.   : csystem C G R D -pmback uL uB -> R D will contain glon, glat; uL, uB will contain uR, uD\n");
+  gprint (GP_ERR, "        :  if input values are vectors, their values are replaced with the transformed values\n");
+  gprint (GP_ERR, "        :  if input values are scalars, the transformed coordinates are printed (disable with -q)\n");
+  gprint (GP_ERR, "        :    and the following output variables will be set (names depend on output system)\n");
+  gprint (GP_ERR, "        :  output positions (RA, DEC), (Lambda, Beta), (gLon, gLat)\n");
+  gprint (GP_ERR, "        :  output motions   (uR, uD),  (uL, uB),       (uLon, uLat)\n");
   return FALSE;
+}

trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

-              r41340
+              r41341
     remove_argument (N, &argc, argv);
     if ((mvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  Vector *WxOutV = NULL; // mask vector
+  if ((N = get_argument (argc, argv, "-Wxvec"))) {
+    remove_argument (N, &argc, argv);
+    if ((WxOutV = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  Vector *WyOutV = NULL; // mask vector
+  if ((N = get_argument (argc, argv, "-Wyvec"))) {
+    remove_argument (N, &argc, argv);
+    if ((WyOutV = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
+  }
 …
       mask = mvec->elements.Int;
+    }
+  }
+  if (WxOutV) {
+    ResetVector (WxOutV, OPIHI_FLT, Ntotal);
+    for (i = 0; i < Ntotal; i++) { WxOutV->elements.Flt[i] = NAN; }
+  }
+  if (WyOutV) {
+    ResetVector (WyOutV, OPIHI_FLT, Ntotal);
+    for (i = 0; i < Ntotal; i++) { WyOutV->elements.Flt[i] = NAN; }
+  }
 …
+        }
+      }
+    }
+  }
+  if (WxOutV) {
+    // save the calculated weights
+    for (i = 0; i < fitdata.Npts; i++) {
+      // fitdata only includes the previously unmasked points
+      int n = fitdata.index[i];
+      WxOutV->elements.Flt[n] = fitdata.Wx[i];
+    }
+  }
+  if (WyOutV) {
+    // save the calculated weights
+    for (i = 0; i < fitdata.Npts; i++) {
+      // fitdata only includes the previously unmasked points
+      int n = fitdata.index[i];
+      WyOutV->elements.Flt[n] = fitdata.Wy[i];
+    }
+  }

trunk/Ohana/src/opihi/cmd.astro/mkgauss.c

-              r41340
+              r41341
   double x, y, r, f, Xo, Yo;
   Buffer *buf;
+  int Normalize = FALSE;
+  if ((N = get_argument (argc, argv, "-norm"))) {
+    Normalize = TRUE;
+    remove_argument (N, &argc, argv);
+  }
   // this should be Nx/2, Ny/2 if not set
 …
   /* f = exp (-r), r = (x^2 / 2Sx) + (y^2 / 2Sy) + Sxy*x*y */
+  double Io = Normalize ? 1.0 / (2.0 * M_PI * Sig_x * Sig_y) : 1.0;
   in = (float *) buf[0].matrix.buffer;
   for (j = 0; j < Ny; j++) {
 …
       y = j - Yo;
       r = 0.5*x*x/Sx + 0.5*y*y/Sy + x*y*Sxy;
       f = exp (-r);
+      f = Io * exp (-r);
       *in += f;
+    }

trunk/Ohana/src/opihi/cmd.astro/region.c

-              r41340
+              r41341
 int region (int argc, char **argv) {
   double Ra, Dec, Radius;
+  double Ra, Dec;
   float dx, dy;
   int N, kapa, NoClear, dXpix, dYpix;
 …
+  }
+  float XSIZE = NAN;
+  if ((N = get_argument (argc, argv, "-xsize"))) {
+    remove_argument (N, &argc, argv);
+    XSIZE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    if (XSIZE <= 0) {
+      gprint (GP_ERR, "ERROR: xsize <= 0\n");
+      return FALSE;
+    }
+  }
+  float YSIZE = NAN;
+  if ((N = get_argument (argc, argv, "-ysize"))) {
+    remove_argument (N, &argc, argv);
+    YSIZE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    if (YSIZE <= 0) {
+      gprint (GP_ERR, "ERROR: xsize <= 0\n");
+      return FALSE;
+    }
+  }
   if ((N = get_argument (argc, argv, "-ew"))) {
     remove_argument (N, &argc, argv);
 …
+  }
+  if ((argc != 4) && (argc != 5)) {
+    gprint (GP_ERR, "USAGE: region Ra Dec Radius [projection] [orientation]\n");
+  int RadiusArg, CtypeArg;
+  if (!isnan(XSIZE) || !isnan(YSIZE)) {
+    RadiusArg = -1;
+    CtypeArg = 3;
+    if ((argc != 3) && (argc != 4)) {
+    gprint (GP_ERR, "USAGE: region Ra Dec [projection] [-xsize deg] [-ysize deg]\n");
+    gprint (GP_ERR, "   OR: region Ra Dec Radius [projection]\n");
     gprint (GP_ERR, "  [-image] [-ew] [+ew] [-ns] [+ns] [-no-clear] [-angle theta]\n");
     gprint (GP_ERR, " current: %f %f (%f x %f) (%s)\n",
 …
              &graphmode.coords.ctype[5]);
     return (FALSE);
+  }
+    }
+  } else {
+    RadiusArg = 3;
+    CtypeArg = 4;
+    if ((argc != 4) && (argc != 5)) {
+      gprint (GP_ERR, "USAGE: region Ra Dec Radius [projection]\n");
+      gprint (GP_ERR, "   OR: region Ra Dec [projection] [-xsize deg] [-ysize deg]\n");
+      gprint (GP_ERR, "  [-image] [-ew] [+ew] [-ns] [+ns] [-no-clear] [-angle theta]\n");
+      gprint (GP_ERR, " current: %f %f (%f x %f) (%s)\n",
+              graphmode.coords.crval1, graphmode.coords.crval2,
+              fabs(graphmode.xmax - graphmode.xmin),
+              fabs(graphmode.ymax - graphmode.ymin),
+              &graphmode.coords.ctype[5]);
+      return (FALSE);
+    }
+  }
   if (!ohana_str_to_radec (&Ra, &Dec, argv[1], argv[2])) return (FALSE);
+  Radius = atof (argv[3]);
+  // region 0 0 sin -- should raise an error (radius = 0 or non-numeric)
+  // I want to be able to support the old style call in which both of these were valid:
+  // region 0 0 90 ait
+  // region 0 0 90      <- uses existing, sticky projection type
+  // but I also want to be able to use
+  // region 0 0 ait -xsize 5
+  // region 0 0 ait -ysize 5
+  double Radius = NAN;
+  if (RadiusArg >= 0) {
+    Radius = atof (argv[RadiusArg]);
+  }
   InitCoords (&graphmode.coords, "DEC--TAN");
+  if (argc == 5) {
+    if (!strcasecmp (argv[4], "TAN"))
+      strcpy (graphmode.coords.ctype, "DEC--TAN");
+    if (!strcasecmp (argv[4], "SIN"))
+      strcpy (graphmode.coords.ctype, "DEC--SIN");
+    if (!strcasecmp (argv[4], "ARC"))
+      strcpy (graphmode.coords.ctype, "DEC--ARC");
+    if (!strcasecmp (argv[4], "STG"))
+      strcpy (graphmode.coords.ctype, "DEC--STG");
+    if (!strcasecmp (argv[4], "ZEA"))
+      strcpy (graphmode.coords.ctype, "DEC--ZEA");
+    if (!strcasecmp (argv[4], "AIT"))
+      strcpy (graphmode.coords.ctype, "DEC--AIT");
+    if (!strcasecmp (argv[4], "GLS"))
+      strcpy (graphmode.coords.ctype, "DEC--GLS");
+    if (!strcasecmp (argv[4], "PAR"))
+      strcpy (graphmode.coords.ctype, "DEC--PAR");
+  }
+  if (argc == CtypeArg + 1) {
+    if (!strcasecmp (argv[CtypeArg], "TAN")) { strcpy (graphmode.coords.ctype, "DEC--TAN"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "SIN")) { strcpy (graphmode.coords.ctype, "DEC--SIN"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "ARC")) { strcpy (graphmode.coords.ctype, "DEC--ARC"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "STG")) { strcpy (graphmode.coords.ctype, "DEC--STG"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "ZEA")) { strcpy (graphmode.coords.ctype, "DEC--ZEA"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "AIT")) { strcpy (graphmode.coords.ctype, "DEC--AIT"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "GLS")) { strcpy (graphmode.coords.ctype, "DEC--GLS"); goto got_ctype; }
+    if (!strcasecmp (argv[CtypeArg], "PAR")) { strcpy (graphmode.coords.ctype, "DEC--PAR"); goto got_ctype; }
+    gprint (GP_ERR, "ERROR: invalid projection type %s\n", argv[CtypeArg]);
+    gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR\n");
+    return FALSE;
+  }
+got_ctype:
   graphmode.coords.crval1 = Ra;
 …
   graphmode.coords.pc2_2 =  cos(Angle*RAD_DEG)*pc2_2;
+  /* ask kapa for coordinate limits, to get the right aspect ratio */
+  // ask kapa for coordinate limits, to get the right aspect ratio
+  // dx, dy are the size of the graph region in pixels
   KapaGetLimits (kapa, &dx, &dy);
   dx = fabs (dx);
 …
   /* define limits for Ra, Dec at center, grid in degrees */
+  if (dy < dx) {
+    graphmode.xmin = -(dx/dy)*Radius;
+    graphmode.ymin = -Radius;
+    graphmode.xmax = (dx/dy)*Radius;
+    graphmode.ymax = Radius;
+  if (RadiusArg >= 0) {
+    // force non-anamorphic projection with Radius set to smaller axis
+    if (dy < dx) {
+      graphmode.xmin = -(dx/dy)*Radius;
+      graphmode.ymin = -Radius;
+      graphmode.xmax = (dx/dy)*Radius;
+      graphmode.ymax = Radius;
+    } else {
+      graphmode.xmin = -Radius;
+      graphmode.ymin = -(dy/dx)*Radius;
+      graphmode.xmax = Radius;
+      graphmode.ymax = (dy/dx)*Radius;
+    }
   } else {
+    graphmode.xmin = -Radius;
+    graphmode.ymin = -(dy/dx)*Radius;
+    graphmode.xmax = Radius;
+    graphmode.ymax = (dy/dx)*Radius;
+  }
+    if (isnan(XSIZE)) {
+      graphmode.xmin = -(dx/dy)*YSIZE/2.0;
+      graphmode.ymin = -YSIZE/2.0;
+      graphmode.xmax = (dx/dy)*YSIZE/2.0;
+      graphmode.ymax = YSIZE/2.0;
+    }
+    if (isnan(YSIZE)) {
+      graphmode.xmin = -XSIZE/2.0;
+      graphmode.ymin = -(dy/dx)*XSIZE/2.0;
+      graphmode.xmax = XSIZE/2.0;
+      graphmode.ymax = (dy/dx)*XSIZE/2.0;
+    }
+    // anamorphic projection:
+    if (!isnan(XSIZE) && !isnan(YSIZE)) {
+      graphmode.xmin = -XSIZE/2.0;
+      graphmode.ymin = -YSIZE/2.0;
+      graphmode.xmax =  XSIZE/2.0;
+      graphmode.ymax =  YSIZE/2.0;
+    }
+  }
   set_variable ("XMIN", graphmode.xmin);

trunk/Ohana/src/opihi/cmd.basic/init.c

r41340	r41341
73	73	{1, "??", list_vars, "list variables *"},
74	74	{1, "#", nop, "a NOP function"},
	75	{1, "##", nop, "a NOP function"},
	76	{1, "###", nop, "a NOP function"},
75	77	{1, "local", local, "define local variables"},
76	78	{1, "macro", macro, "deal with the macros *"},

trunk/Ohana/src/opihi/cmd.basic/list.c

-              r41340
+              r41341
     gprint (GP_ERR, "USAGE: list (root) -split (words)                         : create list from words\n");
     gprint (GP_ERR, "USAGE: list (root) -splitbychar (char) (word) [(words)..] : create list from words\n");
+    gprint (GP_ERR, "USAGE: list (root) -join (varname)                        : convert a list to a single variable (space separated)\n");
     gprint (GP_ERR, "USAGE: list (root) -add (words)                           : extend a list\n");
     gprint (GP_ERR, "USAGE: list (root) -copy (list)                           : copy a list to a new name\n");
     gprint (GP_ERR, "USAGE: list (root) -del (word)                            : delete the entry by value\n");
+    gprint (GP_ERR, "USAGE: list (root) -glob (word)                           : create a list from a file glob\n");
+    gprint (GP_ERR, "USAGE: list (root) -file (word)                           : create a list from lines in a file\n");
+    gprint (GP_ERR, "OPTIONS: -excel                                           : generate names of the for foo:A\n");
+    gprint (GP_ERR, "OPTIONS: -excel-style                                     : generate names of the for foo:A\n");
     return (FALSE);
+  }

trunk/Ohana/src/opihi/cmd.basic/memory.c

-              r41340
+              r41341
   if (argc < 2) goto usage;
+  if (!strcasecmp (argv[1], "max-lines")) {
+    if (argc != 3) goto usage;
+    int MaxLines = atoi(argv[2]);
+    ohana_memdump_set_maxlines (MaxLines);
+    return TRUE;
+  }
   if (!strcasecmp (argv[1], "all")) {
 …
+  }
+  if (!strncasecmp ("strings", argv[1], strlen(argv[1]))) {
+    ohana_memdump_strings_file (stderr, FALSE);
+    return (TRUE);
+  }
 usage:
+  gprint (GP_ERR, "USAGE: memory (all/leaks)\n");
+  gprint (GP_ERR, "USAGE: memory (mode) [-max-lines]\n");
+  gprint (GP_ERR, "USAGE: memory (max-lines) (N)\n");
+  gprint (GP_ERR, " mode options: all, leaks, check, checkfree, variables, vectors, buffers, macros, commands, strings\n");
   return (FALSE);
+}

trunk/Ohana/src/opihi/cmd.basic/test/if.sh

-              r41340
+              r41341
  end
 end
+# check memleaks
+macro memtest1
+ local a b i N
+ $i = 0
+ $a = 1
+ memory check
+ output /dev/null
+ for i 0 1000
+   if ($a == 1)
+    echo "run"
+   end
+ end
+ output stdout
+ memory check
+end
+# check memleaks
+macro memtest2
+ local a b i N
+ $i = 0
+ $a = 1
+ $b = 2
+ memory check
+ output /dev/null
+ for i 0 1000
+   if (($a == 1) && ($b == 2))
+    echo "run"
+   end
+ end
+ output stdout
+ memory check
+end
+# memory test using continue
+macro memtest3
+ local a b i N
+ $i = 0
+ $a = 1
+ $b = 2
+ memory check
+ output /dev/null
+ for i 0 1000
+   if (($a == 1) && ($b == 2))
+    continue
+    echo "run"
+   end
+ end
+ output stdout
+ memory check
+end
+# memory test using continue
+macro memtest4
+ local a b i N bool
+ $i = 0
+ $a = 1
+ $b = 2
+ $bool = 0
+ memory check
+ output /dev/null
+ for i 0 1000
+   $bool = (($a == 1) && ($b == 2))
+ end
+ output stdout
+ memory check
+end

trunk/Ohana/src/opihi/cmd.basic/test/macro.sh

-              r41340
+              r41341
 # Is the macro working?
 macro test_prep
+ $test_var1 = check1
+ $test_var2 = check2
+ $var_count = $0
+end
+# Is the macro working?
+macro test_local
+ local test_var test_var2 var_count
  $test_var1 = check1
  $test_var2 = check2
 …
  end
 end
+# Memory Test for macro
+macro memtest2
+ $i = 0
+ test_prep
+ memory check
+ for i 0 10000
+  test_prep
+ end
+ memory check
+end
+# Memory Test for macro
+macro memtest3
+ $i = 0
+ test_local
+ memory check
+ for i 0 10000
+  test_local
+ end
+ memory check
+end

trunk/Ohana/src/opihi/cmd.data/Makefile

-              r41340
+              r41341
 $(SRC)/init.$(ARCH).o          \
 $(SRC)/accum.$(ARCH).o          \
+$(SRC)/antialias.$(ARCH).o              \
 $(SRC)/applyfit1d.$(ARCH).o     \
 $(SRC)/applyfit2d.$(ARCH).o     \
 …
 $(SRC)/join.$(ARCH).o           \
 $(SRC)/jpeg.$(ARCH).o           \
+$(SRC)/kapamemory.$(ARCH).o             \
 $(SRC)/kern.$(ARCH).o           \
 $(SRC)/keyword.$(ARCH).o        \
 …
 $(SRC)/lookup.$(ARCH).o \
 $(SRC)/matrix.$(ARCH).o \
+$(SRC)/match1d.$(ARCH).o        \
 $(SRC)/match2d.$(ARCH).o        \
 $(SRC)/mkrgb.$(ARCH).o  \
 $(SRC)/mcreate.$(ARCH).o        \
 $(SRC)/medacc.$(ARCH).o \
+$(SRC)/mgaussdev.$(ARCH).o      \
 $(SRC)/mget.$(ARCH).o           \
 $(SRC)/mget3d.$(ARCH).o         \
 …
 $(SRC)/medimage.$(ARCH).o       \
 $(SRC)/medimage_commands.$(ARCH).o \
+$(SRC)/medimage_calc.$(ARCH).o \
 $(SRC)/mset.$(ARCH).o           \
 $(SRC)/needles.$(ARCH).o        \
 …
 $(SRC)/point.$(ARCH).o          \
 $(SRC)/ps.$(ARCH).o             \
+$(SRC)/pdf.$(ARCH).o            \
 $(SRC)/print_vectors.$(ARCH).o  \
 $(SRC)/mprint.$(ARCH).o         \
 …
 $(SRC)/queuesubstr.$(ARCH).o    \
 $(SRC)/queueinit.$(ARCH).o      \
+$(SRC)/queue2book.$(ARCH).o     \
 $(SRC)/radial.$(ARCH).o \
 $(SRC)/rd.$(ARCH).o             \
 …
 $(SRC)/vgauss.$(ARCH).o            \
 $(SRC)/vlorentz.$(ARCH).o          \
+$(SRC)/vsigmoid.$(ARCH).o          \
 $(SRC)/vellipse.$(ARCH).o          \
 $(SRC)/vmaxwell.$(ARCH).o          \
 …
 $(SRC)/vzload.$(ARCH).o            \
 $(SRC)/vpeaks.$(ARCH).o            \
+$(SRC)/vtransitions.$(ARCH).o     \
 $(SRC)/vpop.$(ARCH).o              \
 $(SRC)/vroll.$(ARCH).o             \
 …
 $(SRC)/vsmooth.$(ARCH).o           \
 $(SRC)/vstats.$(ARCH).o            \
+$(SRC)/virls.$(ARCH).o             \
+$(SRC)/vwtmean.$(ARCH).o \
 $(SRC)/xsection.$(ARCH).o          \
 $(SRC)/vsh.$(ARCH).o               \

trunk/Ohana/src/opihi/cmd.data/create.c

-              r41340
+              r41341
   Vector *vec;
+  // create a vector of empty strings
+  if ((N = get_argument (argc, argv, "-str"))) {
+    remove_argument (N, &argc, argv);
+    char *stringValue = NULL;
+    if ((N = get_argument (argc, argv, "-value"))) {
+      remove_argument (N, &argc, argv);
+      stringValue = strcreate (argv[N]);
+      remove_argument (N, &argc, argv);
+    }
+    if (argc != 3) {
+      gprint (GP_ERR, "USAGE: create vector Nelements -value word\n");
+      return (FALSE);
+    }
+    if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    int Nelements = atoi (argv[2]);
+    // a newly reset vector has NULL-valued pointers
+    ResetVector (vec, OPIHI_STR, Nelements);
+    for (int i = 0; i < Nelements; i++) {
+      vec[0].elements.Str[i] = stringValue ? strcreate (stringValue) : strcreate ("");
+    }
+    return TRUE;
+  }
   INT = FALSE;
   if ((N = get_argument (argc, argv, "-int"))) {
 …
     gprint (GP_ERR, "USAGE: create vector start end [delta] [-int]\n");
     gprint (GP_ERR, " -int : resulting vector is integer type (delta must be integer)\n");
+    gprint (GP_ERR, " -str : resulting vector is string type (only give number of elements)\n");
     return (FALSE);
+  }

trunk/Ohana/src/opihi/cmd.data/cut.c

-              r41340
+              r41341
 # include "data.h"
+enum {SUM, MEAN, MEDIAN};
+enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
+double cutstat (float *value, int Nvalue, int mode) {
+  double output = 0.0;
+  if (Nvalue == 0) return NAN;
+  if (mode == MEDIAN) {
+    fsort (value, Nvalue);
+    if (Nvalue % 2) {
+      int Ncenter = Nvalue / 2;
+      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue), "cutstat");
+      output = value[Ncenter];
+    } else {
+      int Ncenter = Nvalue / 2 - 1;
+      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue - 1), "cutstat");
+      output = 0.5*(value[Ncenter] + value[Ncenter + 1]);
+    }
+    return output;
+  }
+  if (mode == INNER) {
+    fsort (value, Nvalue);
+    int Ns = 0, Ne = 0;
+    if (Nvalue % 2) {
+      // for an odd number of points take the same number below
+      // and above the center value
+      int Ncenter = Nvalue / 2;
+      int Nquarter = 0.25*Nvalue;
+      Ns = Ncenter - Nquarter;
+      Ne = Ncenter + Nquarter;
+    } else {
+      // for an even number, the middle lies between two points
+      // take the same number below as above
+      int Ncenter = (int)(Nvalue / 2) - 1;
+      int Nquarter = 0.25*Nvalue;
+      Ns = Ncenter     - Nquarter;
+      Ne = Ncenter + 1 + Nquarter;
+    }
+    int Nv = 0;
+    for (int i = Ns; i <= Ne; i++) {
+      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
+      output += value[i];
+      Nv ++;
+    }
+    output /= Nv;
+    return output;
+  }
+  if (mode == NGOOD) {
+    return Nvalue;
+  }
+  if ((mode == MEAN) || (mode == SUM)) {
+    for (int i = 0; i < Nvalue; i++) {
+      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
+      output += value[i];
+    }
+    if (mode == MEAN) { output /= Nvalue; }
+    return output;
+  }
+  return NAN;
+}
 int cut (int argc, char **argv) {
   int i, j, N, Nx, Ny, Mode;
   float *Vin, *Vbuf, value;
+  float *Vin, *Vbuf;
   int sx, sy, nx, ny;
   Vector *xvec, *yvec;
 …
   Mode = SUM;
+  if ((N = get_argument (argc, argv, "-sum"))) {
+    remove_argument (N, &argc, argv);
+    Mode = SUM;
+  }
   if ((N = get_argument (argc, argv, "-median"))) {
     remove_argument (N, &argc, argv);
 …
     remove_argument (N, &argc, argv);
     Mode = MEAN;
+  }
+  if ((N = get_argument (argc, argv, "-inner"))) {
+    remove_argument (N, &argc, argv);
+    Mode = INNER;
+  }
+  if ((N = get_argument (argc, argv, "-ngood"))) {
+    remove_argument (N, &argc, argv);
+    Mode = NGOOD;
+  }
 …
   Ny = buf[0].matrix.Naxis[1];
+  if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
+    gprint (GP_ERR, "region out of range\n");
+  // ny & nx do not need to be constrained by the buffer, but they do need to be sensible
+  if ((nx < 0) || (ny < 0) || (nx > 1e8) || (ny > 1e8)) {
+    gprint (GP_ERR, "warning : extraction size is crazy: %d,%d\n", nx, ny);
     return (FALSE);
+  }
 …
     for (i = 0; i < nx; i++) {
+      // for out-of-range areas, set the output pixels to NAN
+      if (i + sx < 0) {
+        yvec[0].elements.Flt[i] = NAN;
+        continue;
+      }
+      if (i + sx >= Nx) {
+        yvec[0].elements.Flt[i] = NAN;
+        continue;
+      }
       /* accumulate values */
+      // skip out-of-range areas in y
+      if (sy < 0) {
+        ny += sy;
+        sy = 0;
+      }
+      if (sy > Ny) sy = Ny;
+      if (sy + ny >= Ny) ny = Ny - sy;
       Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i;
+      int Npix = 0;
       for (j = 0; j < ny; j++, Vin += Nx) {
+        Vbuf[j] = *Vin;
+      }
+      /* apply stat of choice */
+      if (Mode == MEDIAN) {
+        fsort (Vbuf, ny);
+        value = Vbuf[(int)(0.5*ny)];
+      } else {
+        value = 0;
+        for (j = 0; j < ny; j++) {
+          value += Vbuf[j];
+        }
+        if (Mode == MEAN) { value /= ny; }
+      }
+      yvec[0].elements.Flt[i] = value;
+        if (!isfinite(*Vin)) continue;
+        Vbuf[Npix] = *Vin;
+        Npix ++;
+      }
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
+    }
     free (Vbuf);
 …
     for (i = 0; i < ny; i++) {
+      // for out-of-range areas, set the output pixels to NAN
+      if (i + sy < 0) {
+        yvec[0].elements.Flt[i] = NAN;
+        continue;
+      }
+      if (i + sy >= Ny) {
+        yvec[0].elements.Flt[i] = NAN;
+        continue;
+      }
       /* accumulate values */
+      // skip out-of-range areas in x
+      if (sx < 0) {
+        nx += sx;
+        sx = 0;
+      }
+      if (sx > Nx) sx = Nx;
+      if (sx + nx >= Nx) nx = Nx - sx;
       Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx;
+      int Npix = 0;
       for (j = 0; j < nx; j++, Vin ++) {
+        Vbuf[j] = *Vin;
+      }
+      /* apply stat of choice */
+      if (Mode == MEDIAN) {
+        fsort (Vbuf, nx);
+        value = Vbuf[(int)(0.5*nx)];
+      } else {
+        value = 0;
+        for (j = 0; j < nx; j++) {
+          value += Vbuf[j];
+        }
+        if (Mode == MEAN) { value /= nx; }
+      }
+      yvec[0].elements.Flt[i] = value;
+        if (!isfinite(*Vin)) continue;
+        Vbuf[Npix] = *Vin;
+        Npix ++;
+      }
+      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
+    }
     free (Vbuf);

trunk/Ohana/src/opihi/cmd.data/dbselect.c

-              r41340
+              r41341
       case FIELD_TYPE_DOUBLE:
       case FIELD_TYPE_TIME:
+      case FIELD_TYPE_TIMESTAMP:
       case FIELD_TYPE_DATE:
       case FIELD_TYPE_DATETIME:
 …
             break;
           case FIELD_TYPE_TIME:
+          case FIELD_TYPE_TIMESTAMP:
           case FIELD_TYPE_DATE:
           case FIELD_TYPE_DATETIME:

trunk/Ohana/src/opihi/cmd.data/extract.c

-              r41340
+              r41341
 # include "data.h"
+/* <from> : source image, must exist
+   <to>   : target image -- if it does not exist, it is created of size (Nx,Ny)
+   sx, sy : source starting coordinate -- need not be on a valid image pixel
+   nx, ny : number of source pixels -- currently must not go out of bounds -> allow out-of-bounds
+   Sx, Sy : target starting coordinate -- need not be on a valid image pixel
+   Nx, Ny : redundant information UNLESS <to> does exist (in which case it is required information) -> make Nx,Ny optional if <to> exists?
+ */
 int extract (int argc, char **argv) {
+  int i, j;
+  float *Vin, *Vout;
+  int sx, sy, nx, ny, NX, NY;
+  int Sx, Sy, Nx, Ny;
+  int N;
   Buffer *in, *out;
+  float initValue = 0.0;
+  int initOutput = FALSE;
+  if ((N = get_argument (argc, argv, "-init"))) {
+    remove_argument (N, &argc, argv);
+    initValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    initOutput = TRUE;
+  }
   if (argc != 11) {
 …
   if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
   NX = in[0].matrix.Naxis[0];
   NY = in[0].matrix.Naxis[1];
+  int NX = in[0].matrix.Naxis[0];
+  int NY = in[0].matrix.Naxis[1];
   sx = atof (argv[3]);
   sy = atof (argv[4]);
   nx = atof (argv[5]);
   ny = atof (argv[6]);
+  int sx = atof (argv[3]);
+  int sy = atof (argv[4]);
+  int nx = atof (argv[5]);
+  int ny = atof (argv[6]);
   Sx = atof (argv[7]);
   Sy = atof (argv[8]);
   Nx = atof (argv[9]);
   Ny = atof (argv[10]);
+  int Sx = atof (argv[7]);
+  int Sy = atof (argv[8]);
+  int Nx = atof (argv[9]);
+  int Ny = atof (argv[10]);
   if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
     gprint (GP_ERR, "mismatch between source and dest regions\n");
+    gprint (GP_ERR, "source pixels extend beyond target pixels\n");
     gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
     return (FALSE);
+    // return (FALSE);
+  }
+  // XXX : allow source region to fall outside source image
   /* region is not on first image */
 …
       (sy > in[0].matrix.Naxis[1])) {
     gprint (GP_ERR, "region outside of source image\n");
     return (FALSE);
+    // return (FALSE);
+  }
-  if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
-    gprint (GP_ERR, "source region larger than dest region\n");
-    return (FALSE);
+  }
   if ((Sx < 0) || (Sy < 0)) {
     gprint (GP_ERR, "dest region out of range\n");
     return (FALSE);
+    // return (FALSE);
+  }
 …
+  }
+  for (j = 0; j < ny; j++) {
+    if (j + sy < 0) continue;
+    if (j + sy >= NY) continue;
+    Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx;
+    Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;
+    for (i = 0; i < nx; i++, Vin++, Vout++) {
+  // (NX,NY) : source image dimensions
+  // (Nx,Ny) : target image dimensions
+  // allow (sx, sy), (Sx, Sy) to be off image
+  // allow (sx+nx, sy+ny) to be off image
+  // if (sx < 0) {
+  //   nx += sx;
+  //   sx = 0;
+  // }
+  // if (sx > NX) sx = NX;
+  //
+  // if (Sx < 0) {
+  //   Nx += sx;
+  //   sx = 0;
+  // }
+  // if (sx > NX) sx = NX;
+  float *Vin = (float *)(in[0].matrix.buffer);
+  float *Vout = (float *)(out[0].matrix.buffer);
+  if (initOutput) {
+    for (int j = 0; j < Ny; j++) {
+      for (int i = 0; i < Nx; i++) {
+        Vout[i + Nx*j] = initValue;
+      }
+    }
+  }
+  int Nps = NX*NY;
+  int Npt = Nx*Ny;
+  for (int j = 0; j < ny; j++) {
+    if (j + sy < 0) continue; // not yet on source image
+    if (j + Sy < 0) continue; // not yet on target image
+    if (j + sy >= NY) continue; // past edge of source image
+    if (j + Sy >= Ny) continue; // past edge of target image
+    // Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx;
+    // Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;
+    for (int i = 0; i < nx; i++) {
       if (i + sx < 0) continue;
       if (i + sx >= NX) continue;
+      *Vout = *Vin;
+      if (i + Sx < 0) continue;
+      if (i + Sx >= Nx) continue;
+      int ps = i + sx + (j + sy)*NX;
+      int pt = i + Sx + (j + Sy)*Nx;
+      myAssert (pt >= 0, "oops 1");
+      myAssert (pt < Npt, "oops 2");
+      myAssert (ps >= 0, "oops 3");
+      myAssert (ps < Nps, "oops 4");
+      Vout[pt] = Vin[ps];
+    }
+  }
 …
+}

trunk/Ohana/src/opihi/cmd.data/init.c

-              r41340
+              r41341
 int accum            PROTO((int, char **));
+int antialias        PROTO((int, char **));
 int applyfit         PROTO((int, char **));
 int applyfit1d       PROTO((int, char **));
 …
 int join             PROTO((int, char **));
 int jpeg             PROTO((int, char **));
+int kapamemory       PROTO((int, char **));
 int kern             PROTO((int, char **));
 int keyword          PROTO((int, char **));
 …
 int lookup           PROTO((int, char **));
 int matrix           PROTO((int, char **));
+int match1d          PROTO((int, char **));
 int match2d          PROTO((int, char **));
 int mkrgb            PROTO((int, char **));
 int mcreate          PROTO((int, char **));
 int medacc           PROTO((int, char **));
+int mgaussdev        PROTO((int, char **));
 int mget             PROTO((int, char **));
 int mget3d           PROTO((int, char **));
 …
 int parity           PROTO((int, char **));
 int point            PROTO((int, char **));
+int pdf              PROTO((int, char **));
 int ps               PROTO((int, char **));
 int vprint           PROTO((int, char **));
 …
 int queuesubstr      PROTO((int, char **));
 int queuesize        PROTO((int, char **));
+int queue2book       PROTO((int, char **));
 int rd               PROTO((int, char **));
 int rdseg            PROTO((int, char **));
 …
 int vgauss           PROTO((int, char **));
 int vlorentz         PROTO((int, char **));
+int vsigmoid         PROTO((int, char **));
 int vellipse         PROTO((int, char **));
 int vmaxwell         PROTO((int, char **));
 …
 int vzload           PROTO((int, char **));
 int vstats           PROTO((int, char **));
+int vstats           PROTO((int, char **));
+int virls            PROTO((int, char **));
+int vwtmean          PROTO((int, char **));
 int vroll            PROTO((int, char **));
 int vshift           PROTO((int, char **));
 int vpeaks           PROTO((int, char **));
+int vtransitions     PROTO((int, char **));
 int vpop             PROTO((int, char **));
 int vsmooth          PROTO((int, char **));
 …
 static Command cmds[] = {
   {1, "accum",        accum,            "accumulate vector values in another vector"},
+  {1, "antialias",    antialias,        "set anti-alias sigma value for display"},
   {1, "applyfit",     applyfit1d,       "apply 1-d fit to new vector"},
   {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
 …
   {1, "join",         join,             "find the join of two ID vectors"},
   {1, "jpeg",         jpeg,             "convert display image to JPEG"},
+  {1, "kapamemory",   kapamemory,       "manage kapa memory dump options"},
   {1, "kern",         kern,             "convolve with 3x3 kernel"},
   {1, "keyword",      keyword,          "extract a FITS keyword from image header"},
 …
   {1, "imcreate",     mcreate,          "create an image"},
   {1, "medacc",       medacc,           "accumulate vector values in another vector"},
+  {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
   {1, "mget",         mget,             "extract a vector from an image"},
   {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
 …
   {1, "medimage",     medimage_command, "median image manipulation"},
   {1, "matrix",       matrix,           "matrix math operations"},
+  {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
   {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
   {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
 …
   {1, "ppm",          jpeg,             "convert display graphic to PPM"},
   {1, "ps",           ps,               "convert display to PostScript"},
+  {1, "pdf",          pdf,              "convert display to PDF"},
   {1, "print_vectors", vprint,          "print a set of vectors"},
   {1, "vprint",       vprint,           "print a set of vectors"},
 …
   {1, "queuesize",    queuesize,        "show queue size"},
   {1, "queuesubstr",  queuesubstr,      "bulk replace strings in queue"},
+  {1, "queue2book",   queue2book,       "convert queue with ipptool output to book"},
+  {1, "ipptool2book", queue2book,       "convert queue with ipptool output to book"},
   {1, "rd",           rd,               "load fits image"},
   {1, "rdseg",        rdseg,            "read a segment of an image from a file"},
 …
   {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
   {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
+  {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
   {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
   {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
 …
   {1, "vload",        vload,            "load vectors as overlay on image display"},
   {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a vector"},
+  {1, "vpeaks",       vpeaks,           "fine coord and flux of peaks in vector"},
+  {1, "vpeaks",       vpeaks,           "find coord and flux of peaks in vector"},
+  {1, "vtransitions", vtransitions,     "find points in vector that cross the transition value"},
   {1, "vpop",         vpop,             "remove first element of a vector"},
   {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
 …
   {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
   {1, "vstats",       vstats,           "statistics on a vector"},
+  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
+  {1, "virls",        virls,            "IRLS mean of a vector"},
   {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
   {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

trunk/Ohana/src/opihi/cmd.data/interpolate_presort.c

-              r41340
+              r41341
   double x0, dx, dy, y0;
   Vector *xout, *yout, *xin, *yin;
+  int N;
+  int FillEnds = FALSE;
+  if ((N = get_argument (argc, argv, "-fill-ends"))) {
+    FillEnds = TRUE;
+    remove_argument (N, &argc, argv);
+  }
   /** check basic syntax **/
 …
   y0 = yin[0].elements.Flt[0];
+  /* in and out vectors are sorted */
+  j = 0;
+  // fill in the start with NANs
+  for (i = 0; (i < xout[0].Nelements) && (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]); i++) {
+    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN;
+  }
   // every output pixel should get a value unless
   // we are below the lowest in or above the highest in
   for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
+  for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
     yout[0].elements.Flt[i] = NAN;
 …
   // fill in the rest with NANs
+  int NinLast = xin[0].Nelements - 1;
   while (i < yout[0].Nelements) {
     yout[0].elements.Flt[i] = NAN;
+    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN;
     i++;
+  }
   return (TRUE);
+}

trunk/Ohana/src/opihi/cmd.data/list_vectors.c

-              r41340
+              r41341
+  }
+  if (vec->type == OPIHI_FLT) {
+    if (Variable) {
+      set_str_variable (Variable, "FLT");
+    } else {
+      gprint (GP_LOG, "%s : FLT\n", argv[1]);
+    }
+  } else {
+    if (Variable) {
+      set_str_variable (Variable, "INT");
+    } else {
+      gprint (GP_LOG, "%s : INT\n", argv[1]);
+    }
+  switch (vec->type) {
+    case OPIHI_FLT:
+      if (Variable) {
+        set_str_variable (Variable, "FLT");
+      } else {
+        gprint (GP_LOG, "%s : FLT\n", argv[1]);
+      }
+      break;
+    case OPIHI_INT:
+      if (Variable) {
+        set_str_variable (Variable, "INT");
+      } else {
+        gprint (GP_LOG, "%s : INT\n", argv[1]);
+      }
+      break;
+    case OPIHI_STR:
+      if (Variable) {
+        set_str_variable (Variable, "STR");
+      } else {
+        gprint (GP_LOG, "%s : STR\n", argv[1]);
+      }
+      break;
+  }
   if (Variable) free (Variable);

trunk/Ohana/src/opihi/cmd.data/medimage_commands.c

-              r41340
+              r41341
 int medimage_add (int argc, char **argv) {
+  int N;
   Buffer *image;
+  Buffer *var = NULL;
+  if ((N = get_argument (argc, argv, "-variance"))) {
+    remove_argument (N, &argc, argv);
+    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
   if (argc != 3) {
     gprint (GP_ERR, "USAGE: medimage add (name) (image)\n");
+    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n");
     gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
+    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
     return FALSE;
+  }
   if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if (var) {
+    if ((var->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
+        (var->matrix.Naxis[1] != image->matrix.Naxis[1])) {
+      gprint (GP_ERR, "variance buffer does not match image buffer dimensions\n");
+      return FALSE;
+    }
+  }
   MedImageType *median = FindMedImage (argv[1]);
 …
   // new image should match existing medimage dimensions
+  // AddMedImage (median, image, var);
   int Ninput = median->Ninput;
   median->Ninput ++;
+  REALLOCATE (median->buffers, float *, median->Ninput);
+  REALLOCATE (median->flx, float *, median->Ninput);
+  REALLOCATE (median->var, float *, median->Ninput);
   ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny);
   memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
+  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  return TRUE;
+}
+int medimage_calc (int argc, char **argv) {
+  int ix, iy, n, N;
+  Buffer *output;
+  int CALC_MEAN = FALSE;
+  if ((N = get_argument (argc, argv, "-mean"))) {
+    CALC_MEAN = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
+    gprint (GP_ERR, "       calculate the median image for the median image set\n");
+    return FALSE;
+  }
+  MedImageType *median = FindMedImage (argv[1]);
+  if (!median) {
+    gprint (GP_ERR, "median image %s not found\n", argv[1]);
+    return FALSE;
+  }
+  if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  int Ninput = median->Ninput;
+  int Nx = median->Nx;
+  int Ny = median->Ny;
+  float *value = NULL;
+  ALLOCATE (value, float, Ninput);
+  gfits_free_matrix (&output->matrix);
+  gfits_free_header (&output->header);
+  if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
+  float *outvalue = (float *) output->matrix.buffer;
+  for (iy = 0; iy < Ny; iy++) {
+    for (ix = 0; ix < Nx; ix++) {
+      int N = 0;
+      int Npix = ix + Nx*iy;
+      for (n = 0; n < Ninput; n++) {
+        float v = median->buffers[n][Npix];
+        if (!isfinite(v)) continue;
+        value[N] = v;
+        N++;
+      }
+      if (N == 0) continue;
+      if (CALC_MEAN) {
+        float sum = 0.0;
+        for (n = 0; n < N; n++) {
+          sum += value[n];
+        }
+        outvalue[Npix] = sum / (float) N;
+      } else {
+        fsort (value, N);
+        outvalue[Npix] = value[(int)(0.5*N)];
+      }
+    }
+  median->var[Ninput] = NULL;
+  if (var) {
+    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
+    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  }
 …
 /*
 int medimage_save (int argc, char **argv) {
+   int medimage_save (int argc, char **argv) {
   int N;
+   int N;
   int APPEND = FALSE;
   if ((N = get_argument (argc, argv, "-append"))) {
     APPEND = TRUE;
     remove_argument (N, &argc, argv);
+  }
+   int APPEND = FALSE;
+   if ((N = get_argument (argc, argv, "-append"))) {
+   APPEND = TRUE;
+   remove_argument (N, &argc, argv);
+   }
   if (argc != 3) {
     gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
     return FALSE;
+  }
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
+   return FALSE;
+   }
   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
     gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
     return (FALSE);
+  }
   return TRUE;
+}
+   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
+   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
 int medimage_load (int argc, char **argv) {
+   int medimage_load (int argc, char **argv) {
   if (argc != 3) {
     gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
     return FALSE;
+  }
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
+   return FALSE;
+   }
   if (!LoadMedImage(argv[2], argv[1])) {
     gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
     return (FALSE);
+  }
   return TRUE;
+}
+   if (!LoadMedImage(argv[2], argv[1])) {
+   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
 */

trunk/Ohana/src/opihi/cmd.data/peak.c

-              r41340
+              r41341
       if (*X < start) continue;
       if (*X > end) continue;
+      if (!isfinite(*Y)) continue;
       if (*Y < ymax) continue;
       xmax = *X;
 …
       if (*X < start) continue;
       if (*X > end) continue;
+      if (!isfinite(*Y)) continue;
       if (*Y < ymax) continue;
       xmax = *X;

trunk/Ohana/src/opihi/cmd.data/print_vectors.c

-              r41340
+              r41341
         if (vec[i][0].type == OPIHI_FLT) {
           gprint (GP_LOG, "%f ", vec[i][0].elements.Flt[j]);
+        } else {
+        }
+        if (vec[i][0].type == OPIHI_INT) {
           gprint (GP_LOG, OPIHI_INT_FMT" ", vec[i][0].elements.Int[j]);
+        }
+        }
+        if (vec[i][0].type == OPIHI_STR) {
+          gprint (GP_LOG, "%s ", vec[i][0].elements.Str[j]);
+        }
+      }
+    }

trunk/Ohana/src/opihi/cmd.data/read_vectors.c

-              r41340
+              r41341
 // vector types
 enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_CHAR};
+enum {COLTYPE_NONE, COLTYPE_FLT, COLTYPE_INT, COLTYPE_TIME, COLTYPE_DATE, COLTYPE_CHAR, COLTYPE_STR, COLTYPE_HMS};
 static int      Nvec     = 0;
 …
     gprint (GP_ERR, "     -skip N : skip N lines before reading\n");
     gprint (GP_ERR, "     column names may include a type: name:type\n");
     gprint (GP_ERR, "       type is int, float, char, time\n");
+    gprint (GP_ERR, "       type is int, float, char, str, time, date, hms\n");
     gprint (GP_ERR, "       for char, values are placed into a list $name:0 - $name:n\n");
+    gprint (GP_ERR, "       for str, values are placed into a string-typed vector\n");
+    gprint (GP_ERR, "       for hms, values are sexigesimal values HH:MM:SS (good for degrees as well)\n");
+    gprint (GP_ERR, "       for date, values are human-readable date strings YYYY/MM/DD\n");
     gprint (GP_ERR, "       for time, values are human-readable date/time strings YYYY/MM/DD,hh:mm:ss\n");
     gprint (GP_ERR, "       the resulting vector is a float based on the TIMEFORMAT, TIMEREF values\n");
+    gprint (GP_ERR, "         date & time values are stored as a float using the TIMEFORMAT, TIMEREF values for the conversion\n");
     gprint (GP_ERR, "     -fits options:\n");
 …
       if (!strcasecmp(ptr, "int"))   { coltype[i] = COLTYPE_INT; }
       if (!strcasecmp(ptr, "char"))  { coltype[i] = COLTYPE_CHAR; }
+      if (!strcasecmp(ptr, "str"))   { coltype[i] = COLTYPE_STR; }
+      if (!strcasecmp(ptr, "date"))  { coltype[i] = COLTYPE_DATE; }
       if (!strcasecmp(ptr, "time"))  { coltype[i] = COLTYPE_TIME; }
+      if (!strcasecmp(ptr, "hms"))   { coltype[i] = COLTYPE_HMS; }
       if (!coltype[i]) goto bad_colname;
+    }
 …
         break;
       case COLTYPE_FLT:
+      case COLTYPE_HMS:
       case COLTYPE_TIME:
+        // note that COLTYPE_TIME is a type of float
+      case COLTYPE_DATE:
         ResetVector (vec[i], OPIHI_FLT, NELEM);
+        break;
+      case COLTYPE_STR:
+        ResetVector (vec[i], OPIHI_STR, NELEM);
         break;
       case COLTYPE_CHAR:
 …
         time_t tvalue;
         int readStatus = FALSE;
+        int dataStatus = FALSE;
         // need to make the if cases for coltype[i]
         switch (coltype[i]) {
 …
               break;
+            }
+          case COLTYPE_STR:
+            {
+              // I need to get an isolated word in 'value' with the string value of this field
+              char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0);
+              char *value = NULL;
+              if (IsCSV) {
+                char *end = parse_nextword_csv (ptr);
+                if (end) {
+                  value = end ? strncreate (ptr, end - ptr) : strcreate (ptr);
+                }
+              } else {
+                value = thisword(ptr);
+              }
+              vec[i][0].elements.Str[Nelem] = value; // needs to be freed later.
+              break;
+            }
           case COLTYPE_FLT:
             readStatus = IsCSV ? dparse_csv (&dvalue, col[i], c0) : dparse (&dvalue, col[i], c0);
 …
             vec[i][0].elements.Flt[Nelem] = readStatus ? dvalue : NAN;
             break;
+          case COLTYPE_DATE: {
+            char *string = NULL;
+            if (IsCSV) {
+              char *ptr = ptrparse_csv (col[i], c0);
+              string = strcreate (ptr); // create separate string (NULL-safe)
+              ptr = (string == NULL) ? NULL : strchr (string, ',');
+              if (ptr != NULL) *ptr = 0; // place an EOL here so parsing does not go past the comma
+            } else {
+              char *ptr = ptrparse (col[i], c0); // NULL-safe
+              string = getword (ptr); // NULL-safe
+            }
+            tvalue = ohana_date_to_sec (string); // NULL-safe (returns 0)
+            dvalue = TimeValue (tvalue, TimeReference, TimeFormat);
+            vec[i][0].elements.Flt[Nelem] = string ? dvalue : NAN;
+            FREE (string);
+            break;
+          }
+          case COLTYPE_HMS: {
+            char *ptr = IsCSV ? ptrparse_csv (col[i], c0) : ptrparse (col[i], c0);
+            char *string = strcreate (ptr); // make a copy so we can munge it in ohana_dms_to_ddd
+            dataStatus = ohana_dms_to_ddd (&dvalue, string);
+            vec[i][0].elements.Flt[Nelem] = ((string != NULL) && dataStatus) ? dvalue : NAN;
+            FREE (string);
+            break;
+          }
+        }
         if (!readStatus && VERBOSE) {
 …
               break;
             case COLTYPE_FLT:
+            case COLTYPE_HMS:
             case COLTYPE_TIME:
+            case COLTYPE_DATE:
               ResetVector (vec[i], OPIHI_FLT, NELEM);
+              break;
+            case COLTYPE_STR:
+              ResetVector (vec[i], OPIHI_STR, NELEM);
               break;
             case COLTYPE_CHAR:
 …
         break;
       case COLTYPE_FLT:
+      case COLTYPE_HMS:
       case COLTYPE_TIME:
+      case COLTYPE_DATE:
         ResetVector (vec[i], OPIHI_FLT, Nelem);
+        break;
+      case COLTYPE_STR:
+        ResetVector (vec[i], OPIHI_STR, Nelem);
         break;
       case COLTYPE_CHAR:
 …
+  }
+  // if CharAsVectors, char fields will be saved as vectors NAME:0 -- NAME:n for n characters
+  // else char fields will be saved as $NAME:0 - $NAME:m for m rows
+  // by default, we now store a string-type field as a string-type vector.
+  // if CharAsList is selected (and My < 10000), char fields will be saved as $NAME:0 - $NAME:m for m rows
+  // if CharAsVectors is selected, char fields will be saved as vectors NAME:0 -- NAME:n for n characters
   // if (Ny > 10000), force CharAsVectors
   int CharAsVectors = FALSE;
 …
     remove_argument (N, &argc, argv);
     CharAsVectors = TRUE;
+  }
+  int CharAsList = FALSE;
+  if ((N = get_argument (argc, argv, "-char-list"))) {
+    remove_argument (N, &argc, argv);
+    CharAsList = TRUE;
+  }
 …
     if (!FITS_TRANSPOSE) {
       // read string column into a list rather than a vector
+      if (!strcmp (type, "char") && !CharAsVectors && (Ny < 3000)) {
+        char *fieldName = argv[i];
+        char *Ptr = data;
+        char varname[1024];  // used as a buffer for the names of string fields
+        for (j = 0; j < Ny; j++) {
+          set_list_varname (varname, fieldName, j, FALSE);
+          char *value = strncreate (&Ptr[j*Nval], Nval);
+          // replace instances of $ with _
+          char *p = strchr (value, '$');
+          while (p) {
+            *p = '_';
+            p = strchr (p, '$');
+      if (!strcmp (type, "char")) {
+        // save char-type field as a List:
+        if (CharAsList && (Ny < 3000)) {
+          char *fieldName = argv[i];
+          char *Ptr = data;
+          char varname[1024];  // used as a buffer for the names of string fields
+          for (j = 0; j < Ny; j++) {
+            set_list_varname (varname, fieldName, j, FALSE);
+            char *value = strncreate (&Ptr[j*Nval], Nval);
+            // replace instances of $ with _
+            char *p = strchr (value, '$');
+            while (p) {
+              *p = '_';
+              p = strchr (p, '$');
+            }
+            set_str_variable (varname, value);
+            free (value);
+          }
+          set_str_variable (varname, value);
+          free (value);
+          sprintf (varname, "%s:n", fieldName);
+          set_int_variable (varname, Ny);
+          continue;
+        }
+        sprintf (varname, "%s:n", fieldName);
+        set_int_variable (varname, Ny);
+        continue;
+        // save char-type field as a string-type vector:
+        if (!CharAsVectors) {
+          Vector *myVector = NULL;
+          sprintf (name, "%s", argv[i]);
+          if ((myVector = SelectVector (name, ANYVECTOR, TRUE)) == NULL) ESCAPE ("bad vector name\n");
+          ResetVector (myVector, OPIHI_STR, Ny);
+          char *Ptr = data;
+          for (j = 0; j < Ny; j++) {
+            myVector[0].elements.Str[j] = strncreate (&Ptr[j*Nval], Nval);
+            // replace instances of $ with _
+            char *p = strchr (myVector[0].elements.Str[j], '$');
+            while (p) { *p = '_'; p = strchr (p, '$'); }
+          }
+          continue;
+        }
+      }
       // define the multifield vector names (Nval vectors x Ny elements)
+      // CharAsVectors is handled below automatically
       ALLOCATE (vec, Vector *, Nval);
       for (j = 0; j < Nval; j++) {

trunk/Ohana/src/opihi/cmd.data/roll.c

-              r41340
+              r41341
 int roll (int argc, char **argv) {
-  int Nbytes, Nextra;
-  int dX, dx, dy, nx, ny;
   Buffer *buf;
+  // this function is probably not needed
+  gprint (GP_ERR, "ERROR: use 'shift' instead of 'roll'\n");
+  return (FALSE);
   if (argc != 4) {
 …
   if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  dx = atof (argv[2]);
+  dy = atof (argv[3]);
+  if (dy != 0) {
+    gprint (GP_ERR, "only x rolls implemented for now\n");
+  int dx = atoi (argv[2]);
+  int dy = atoi (argv[3]);
+  /* if (dx,dy < 0), we are moving the start position backwards,
+     if (dx,dy > 0), we are moving the start position forward */
+  int Nx = buf[0].matrix.Naxis[0];
+  int Ny = buf[0].matrix.Naxis[1];
+  if (dy == 0) {
+    // moves in just dx can be done row-by-row
+    int dX = abs(dx);
+    int Nbytes = (Nx * Ny - dX) * sizeof (float);
+    int Nextra = dX * sizeof (float);
+    if (dx < 0) {
+      memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
+      memset (&buf[0].matrix.buffer[Nbytes], 0, Nextra);
+    } else {
+      memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
+      memset (buf[0].matrix.buffer, 0, Nextra);
+    }
+    return TRUE;
+  }
+  /* if (dx < 0), we are moving the start position back by dx pixels,
+     if (dx > 0), we are moving the start position forward by dx pixels */
+  nx = buf[0].matrix.Naxis[0];
+  ny = buf[0].matrix.Naxis[1];
+  if (dx == 0) {
+    // moves in just dy can be done row-by-row
+  dX = abs(dx);
+  Nbytes = nx * ny * sizeof (float);
+  Nextra = (nx * ny + dX) * sizeof (float);
+    ALLOCATE_PTR (output, float, Nx*Ny);
   REALLOCATE (buf[0].matrix.buffer, char, Nextra);
+    float *input = (float *) buf[0].matrix.buffer;
+  if (dx < 0) {
+    memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
+  } else {
+    memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
+    for (int iy = 0; iy < Ny; iy++) {
+      if (iy + dy < 0) continue;
+      if (iy + dy >= Ny) continue;
+      memcpy (&output[iy*Nx], &input[(iy + dy)*Nx], Nx);
+    }
+    if (dy > 0) {
+      memset (output, 0, dy*Nx*sizeof(float));
+    } else {
+      int Nbytes = Nx * (Ny - abs(dy)) * sizeof (float);
+      int Nextra = abs(dy) * Nx * sizeof (float);
+      memset (&output[Nbytes], 0, Nextra);
+    }
+    free (buf[0].matrix.buffer);
+    buf[0].matrix.buffer = (char *) output;
+    return TRUE;
+  }
+  REALLOCATE (buf[0].matrix.buffer, char, Nbytes);
+  ALLOCATE_PTR (output, float, Nx*Ny);
+  float *input = (float *) buf[0].matrix.buffer;
+  return (TRUE);
+  for (int iy = 0; iy < Ny; iy++) {
+    if (iy + dy < 0) continue;
+    if (iy + dy >= Ny) continue;
+    for (int ix = 0; ix < Nx; ix++) {
+      if (ix + dx < 0) continue;
+      if (ix + dx >= Nx) continue;
+      output[ix + iy*Nx] = input[(ix + dx) + (iy + dy)*Nx];
+    }
+  }
+  free (buf[0].matrix.buffer);
+  buf[0].matrix.buffer = (char *) output;
+  return TRUE;
+}

trunk/Ohana/src/opihi/cmd.data/rotate.c

-              r41340
+              r41341
 int rotate (int argc, char **argv) {
   int i, j, NX, NY, X, Y, Lx, Ly, N;
+  int i, j, NX, NY;
   float *in_buff, *out_buff, *c;
   double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1;
+  double Xo, Yo, x, y, X1, Y1;
   double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22;
   Buffer *buf;
   Xo = 0;
   Yo = 0;
   if ((N = get_argument (argc, argv, "-center"))) {
     remove_argument (N, &argc, argv);
     Xo  = atof(argv[N]);
     remove_argument (N, &argc, argv);
     Yo  = atof(argv[N]);
     remove_argument (N, &argc, argv);
+  }
+//  Xo = 0;
+//  Yo = 0;
+//  if ((N = get_argument (argc, argv, "-center"))) {
+//    remove_argument (N, &argc, argv);
+//    Xo  = atof(argv[N]);
+//    remove_argument (N, &argc, argv);
+//    Yo  = atof(argv[N]);
+//    remove_argument (N, &argc, argv);
+//  }
   if (argc != 3) {
 …
+  }
   angle = atof (argv[2]);
   CosAngle = cos (angle*RAD_DEG);
   SinAngle = sin (angle*RAD_DEG);
+  double angle = atof (argv[2]);
+  double CosAngle = cos (angle*RAD_DEG);
+  double SinAngle = sin (angle*RAD_DEG);
   gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);
   Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
   Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
   dX = MAX(0,NY*SinAngle);
+  dY = MAX(0,-NX*SinAngle);
   /*
   gprint (GP_ERR, "%f %f  -->  ", Xo, Yo);
   X1 = Xo*CosAngle - Yo*SinAngle + dX;
+  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
   gprint (GP_ERR, "%f %f\n", X1, Y1);
   */
+  // we are rotating about the center pixel, (NX/2, NY/2),
+  // but we are then putting the result in a new image
+  // of size (Lx,Ly).
+  // (x,y) = (i - Nx/2),(j - Ny/2)
+  // (x',y') = R(theta) (x,y)
+  // (I,J) = (x' + Lx/2),(y' + Ly/2)
+  int Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
+  int Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
   /* fix reference pixel */
   gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
   gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
+  /*
+  X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle;
+  Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle;
+  */
+  X1 = Xo*CosAngle - Yo*SinAngle + dX;
+  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
+  X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2;
+  Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2;
   gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
   gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
 …
   gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
   gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
   PC11 = pc11*CosAngle - pc21*SinAngle;
   PC12 = pc12*CosAngle - pc22*SinAngle;
   PC21 = pc21*CosAngle + pc11*SinAngle;
   PC22 = pc22*CosAngle + pc12*SinAngle;
+  PC11 = pc11*CosAngle + pc21*SinAngle;
+  PC12 = pc12*CosAngle + pc22*SinAngle;
+  PC21 = pc21*CosAngle - pc11*SinAngle;
+  PC22 = pc22*CosAngle - pc12*SinAngle;
   gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
   gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
 …
   gfits_create_matrix (&buf[0].header, &buf[0].matrix);
   gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
   out_buff = (float *)buf[0].matrix.buffer;
   for (j = 0; j < Ly; j++) {
     for (i = 0; i < Lx; i++, out_buff++) {
+      x = (i - dX)*CosAngle + (j - dY)*SinAngle;
+      y = (dX - i)*SinAngle + (j - dY)*CosAngle;
+      X = (int) x;
+      Y = (int) y;
+      if (X < 0) continue;
+      if (X >= NX - 1) continue;
+      if (Y < 0) continue;
+      if (Y >= NY - 1) continue;
+      c = &in_buff[X + NX*Y];
+      fx = x - X;
+      fy = y - Y;
+      float xo = (i - Lx/2);
+      float yo = (j - Ly/2);
+      x =  xo*CosAngle + yo*SinAngle;
+      y = -xo*SinAngle + yo*CosAngle;
+      int I = x + NX/2;
+      int J = y + NY/2;
+      if (I < 0) continue;
+      if (I >= NX - 1) continue;
+      if (J < 0) continue;
+      if (J >= NY - 1) continue;
+      c = &in_buff[I + NX*J];
+      int X = (int) x;
+      int Y = (int) y;
+      double fx = x - X;
+      double fy = y - Y;
       *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
+//    *out_buff = c[0];
+    }
+  }

trunk/Ohana/src/opihi/cmd.data/spline.c

-              r41340
+              r41341
 int spline_load (int argc, char **argv);
 int spline_save (int argc, char **argv);
+int spline_print (int argc, char **argv);
 int spline_delete (int argc, char **argv);
 int spline_rename (int argc, char **argv);
 …
   {1, "load",       spline_load,       "write a spline to a FITS file"},
   {1, "save",       spline_save,       "read a spline from a FITS file"},
+  {1, "print",      spline_print,      "print a spline to stdout"},
   {1, "delete",     spline_delete,     "delete a spline"},
   {1, "rename",     spline_rename,     "rename a spline"},
 …
   gprint (GP_ERR, "    spline load   (spline) (filename)            : load a spline from a FITS file\n");
   gprint (GP_ERR, "    spline save   (spline) (filename) [-append]  : save a spline in FITS format\n");
+  gprint (GP_ERR, "    spline print  (spline)                       : print a spline to stdout\n");
   return FALSE;

trunk/Ohana/src/opihi/cmd.data/spline_commands.c

-              r41340
+              r41341
   if (xvec->Nelements != yvec->Nelements) {
     gprint (GP_ERR, "x and y vectors must be the same length\n");
+    return FALSE;
+  }
+  if (xvec->Nelements < 3) {
+    gprint (GP_ERR, "cannot make a spline with fewer than 3 knots\n");
     return FALSE;
+  }
 …
+}
+int spline_print (int argc, char **argv) {
+  int i;
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: spline print (name)\n");
+    return FALSE;
+  }
+  Spline *myspline = FindSpline (argv[1]);
+  if (!myspline) {
+    gprint (GP_ERR, "spline %s not found\n", argv[1]);
+    return (FALSE);
+  }
+  for (i = 0; i < myspline->Nknots; i++) {
+    gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]);
+  }
+  return TRUE;
+}
 int spline_load (int argc, char **argv) {
 …
 int spline_delete (int argc, char **argv) {
+  int status;
+  Spline *spline;
+  int N, status;
+  int QUIET = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    QUIET = TRUE;
+    remove_argument (N, &argc, argv);
+  }
   if (argc != 2) {
 …
+  }
   spline = FindSpline (argv[1]);
+  Spline *spline = FindSpline (argv[1]);
   if (spline == NULL) {
+    if (QUIET) return TRUE;
     gprint (GP_ERR, "spline %s not found\n", argv[1]);
     return FALSE;
 …
 int spline_rename (int argc, char **argv) {
-  Spline *spline;
   if (argc != 3) {
     gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n");
 …
+  }
   spline = FindSpline (argv[1]);
+  Spline *spline = FindSpline (argv[1]);
   if (spline == NULL) {
     gprint (GP_ERR, "spline %s not found\n", argv[1]);
     return FALSE;
+  }
+  Spline *newname = FindSpline (argv[2]);
+  if (newname != NULL) {
+    DeleteSpline (newname);
+  }
 …
   return TRUE;
+}
-/*
-int spline_listspline (int argc, char **argv) {
-  Spline *spline;
-  if (argc != 2) {
-    gprint (GP_ERR, "USAGE: spline listspline (spline)\n");
-    return FALSE;
+  }
-  spline = FindSpline (argv[1]);
-  if (spline == NULL) {
-    gprint (GP_ERR, "spline %s not found\n", argv[1]);
-    return FALSE;
+  }
-  ListPages (spline);
-  return TRUE;
+}
-*/

trunk/Ohana/src/opihi/cmd.data/test/cut.sh

-              r41340
+              r41341
  memtest1
 end
+macro test2
+  mgaussdev z 50 100 0.0 1.0
+  cut z x y_sum    x 0 0 z[][0] z[0][]
+  cut z x y_median x 0 0 z[][0] z[0][] -median
+  cut z x y_mean   x 0 0 z[][0] z[0][] -mean
+  cut z x y_inner  x 0 0 z[][0] z[0][] -inner
+end
+macro checkrange
+  if ($0 != 2)
+    echo "USAGE: checkrange (Npts)"
+    break
+  end
+  $Npts = $1
+  if ($Npts % 2)
+    $Ncenter = int($Npts / 2)
+    $Nquarter = int(0.25 * $Npts)
+    $Ns = $Ncenter - $Nquarter
+    $Ne = $Ncenter + $Nquarter
+  else
+    $Ncenter = int($Npts / 2) - 1
+    $Nquarter = int(0.25 * $Npts)
+    $Ns = $Ncenter - $Nquarter
+    $Ne = $Ncenter + $Nquarter + 1
+  end
+  if ($Ns < 0)
+    echo "error: Ns < 0: $Ns $Ncenter $Nquarter"
+  end
+  if ($Ne >= $Npts)
+    echo "error: Ne >= Npts: $Ne $Ncenter $Nquarter"
+  end
+  echo "$Npts : $Ncenter : $Nquarter : $Ns $Ne"
+  for i 0 $Npts
+    echo $i {($i >= $Ns) && ($i <= $Ne)}
+  end
+end
 # Test if cut works
 …
  end
 end
+macro memtest1
+ $i = 0
+ mcreate tim 100 100
+ cut tim xdir imx X 40 4 60 6
+ # do one to set up memory that should stay used
+ memory check
+ for i 0 100
+   cut tim xdir imx X 40 4 60 6
+ end
+ memory check
+ for i 0 100
+   cut tim xdir imx y 40 4 60 6
+ end
+ memory check
+ for i 0 100
+   cut -median tim xdir imx X 40 4 60 6
+ end
+ memory check
+ for i 0 100
+   cut -median tim xdir imx y 40 4 60 6
+ end
+ memory check
+ for i 0 100
+   cut -mean tim xdir imx X 40 4 60 6
+ end
+ memory check
+ for i 0 100
+   cut -mean tim xdir imx y 40 4 60 6
+ end
+ memory check
+end

trunk/Ohana/src/opihi/cmd.data/test/delete.sh

-              r41340
+              r41341
 end
+# Test delete
+macro test2
+ $PASS = 1
+ for i 0 1000
+   create v 0 200
+   delete v
+ end
+end

trunk/Ohana/src/opihi/cmd.data/test/medimage.sh

-              r41340
+              r41341
+macro go
+ mcreate a 30 30
+ for i 0 40
+  set a$i = zero(a) + $i
+ end
+ for i 0 40
+  medimage add t1 a$i
+ end
+end
+macro test.mean
+ $Nsample = 16
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 1.0
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ medimage calc t1 T -mean
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "expect {1/sqrt($Nsample)} : $C1"
+ plot -c red -x line x yf
+end
+macro test.median
+ $Nsample = 16
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 1.0
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+ end
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)"
+ plot -c red -x line x yf
+end
+macro test.wtmean
+ $Nsample = 8
+ $sig1 = 1.0
+ $sig2 = 3.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig1
+  set v = $sig1^2 + zero(t)
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig2
+  set v = $sig2^2 + zero(t)
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -wtmean
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
+ echo "expect {1/sqrt($S1)} : $C1"
+ plot -c red -x line x yf
+end
+macro test.irls
+ medimage delete -q t1
+ $Nsample = 16
+ $sig = 1.0
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  set v = $sig^2 + zero(t)
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ plot -c red -x line x yf
+ # get stats for irls
+ medimage calc t1 Ti -irls
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ plot -c red -x line x yf
+end
+###################33
+macro test.mean.var
+ $Nsample = 64
+ $sig = 2.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 100 100 0.0 $sig
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ medimage calc t1 T -mean -variance Tv
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+ imhist Tv xv yv -range -1 4 -delta 0.1
+ lim -n 2 xv yv; clear; box; plot xv yv -x hist
+ stat Tv
+ echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
+end
+macro test.median.var
+ $Nsample = 64
+ $sig = 2.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  medimage add t1 t
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -variance Tv
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+ imhist Tv xv yv -range -1 4 -delta 0.1
+ lim -n 2 xv yv; clear; box; plot xv yv -x hist
+ stat Tv
+ echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
+end
+macro test.wtmean.var
+ $Nsample = 32
+ $sig1 = 1.0
+ $sig2 = 1.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig1
+  set v = $sig1^2 + zero(t)
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig2
+  set v = $sig2^2 + zero(t)
+  medimage add t1 t -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
+ # (somewhat higher scatter)
+ medimage calc t1 T -wtmean -variance Tv
+ imhist -q T x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ plot -c red -x line x yf
+ stat -q Tv
+ $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
+ echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)}
+end
+macro test.irls.var
+ $Nsample = 16
+ $sig = 1.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 50 50 0.0 $sig
+  set v = $sig^2 + zero(t)
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ plot -c red -x line x yf
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ plot -c red -x line x yf
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
+end
+macro test.irls.boot.var
+ $Nsample = 64
+ $sig = 1.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 200 200 0.0 $sig
+  set v = $sig^2 + zero(t)
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+  mgaussdev noise 200 200 0.0 0.5
+  set ts = ts + noise
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  # echo $C1
+ end
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+ imhist -q Tm x y -range -10 10 -delta 0.02
+ lim -n 1 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig / sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from straight stdev: $C1"
+ # stats Tm
+ plot -c red -x line x yf
+ # get stats for irls
+ date
+ medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100
+ date
+ imhist -q Ti x y -range -10 10 -delta 0.02
+ lim -n 2 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig / sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss -q x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
+ # stats Ti
+ plot -c red -x line x yf
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range 0 {5*$sig/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
+end
+##############################
+macro test.irls.boot.test
+ $Nsample = 100
+ $sig1 = 1.0
+ medimage delete -q t1
+ for i 0 $Nsample
+  mgaussdev t 100 100 0.0 $sig1
+  set v = $sig1^2 + zero(t)
+  medimage add t1 t -variance v
+ end
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv -bootstrap
+ imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01
+ lim -n 2 x y; clear; box; plot -x hist x y
+ peak -q x y
+ $C0 = $peakpos
+ $C1 = 1.5*$sig1/sqrt($Nsample)
+ $C2 = $peakval
+ $C3 = 0
+ vgauss x y con yf
+ echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})"
+ # stats Ti
+ plot -c red -x line x yf
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range 0 {5*$sig1/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+ stat -q irls_npt
+ $Npix = $MEAN
+ stat -q Tv
+ echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})"
+end
+macro test.irls.range.var
+ medimage delete -q t1
+ for i 0 8
+  mgaussdev t 50 50 0.0 1.0
+  set v = 1.0 + zero(t)
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  echo $C1
+ end
+ for i 0 8
+  mgaussdev t 50 50 0.0 3.0
+  set v = 3.0 + zero(t)
+  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
+  set ts = t + bad
+  medimage add t1 ts -variance v
+  imhist -q t x y -range -10 10 -delta 0.1
+  $C0 = 0
+  $C1 = 1.5
+  $C2 = 400
+  $C3 = 0
+  vgauss -q x y con yf
+  echo $C1
+ end
+ # get stats for straight mean:
+ medimage calc t1 Tm -mean -variance Tv
+ imhist -q Tm x y -range -10 10 -delta 0.1
+ lim -n 1 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo $C1
+ stats Tm
+ plot -c red -x line x yf
+ # get stats for irls
+ medimage calc t1 Ti -irls -variance Tv
+ imhist -q Ti x y -range -10 10 -delta 0.1
+ lim -n 2 x y; clear; box; plot -x hist x y
+ $C0 = 0
+ $C1 = 1.5
+ $C2 = 400
+ $C3 = 0
+ vgauss -q x y con yf
+ echo $C1
+ stats Ti
+ plot -c red -x line x yf
+ stat -q Tv
+ echo $C1 vs {sqrt($MEDIAN)}
+ set dTv = sqrt(Tv)
+ imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
+end

trunk/Ohana/src/opihi/cmd.data/threshold.c

-              r41340
+              r41341
     remove_argument (N, &argc, argv);
     BinMax = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float ValMin = NAN;
+  float ValMax = NAN;
+  if ((N = get_argument (argc, argv, "-vrange"))) {
+    remove_argument (N, &argc, argv);
+    ValMin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ValMax = atof (argv[N]);
     remove_argument (N, &argc, argv);
+  }
 …
+  }
+  // ValMin, ValMax not compatible with BinMin, BinMax
+  if (isfinite(ValMin)) {
+    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMin);
+    BinMin = (vecx[0].elements.Flt[binTest+1] == ValMin) ? binTest : binTest + 1;
+    // XXX dangerous: check that this does not run off the end, etc
+  }
+  if (isfinite(ValMax)) {
+    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMax);
+    BinMax = (vecx[0].elements.Flt[binTest+1] == ValMax) ? binTest : binTest + 1;
+    // XXX dangerous: check that this does not run off the end, etc
+  }
   // use bisection to find the value
 …
   if (!REVERSE) {
+    // this algorithm assumes vecy[BinMax] > threshold, vecy[BinMin] < threshold
+    if (vecy[0].elements.Flt[BinMin] > value) {
+      gprint (GP_ERR, "ERROR: all values above threshold\n");
+      return FALSE;
+    }
+    if (vecy[0].elements.Flt[BinMax] < value) {
+      gprint (GP_ERR, "ERROR: all values below threshold\n");
+      return FALSE;
+    }
     while (Nhi - Nlo > 10) {
       N = 0.5*(Nlo + Nhi);
 …
+    }
   } else {
+    // this algorithm assumes vecy[BinMin] > threshold, vecy[BinMax] < threshold
+    if (vecy[0].elements.Flt[BinMin] < value) {
+      gprint (GP_ERR, "ERROR: all values below threshold\n");
+      return FALSE;
+    }
+    if (vecy[0].elements.Flt[BinMax] > value) {
+      gprint (GP_ERR, "ERROR: all values above threshold\n");
+      return FALSE;
+    }
     while (Nhi - Nlo > 10) {
       N = 0.5*(Nlo + Nhi);

trunk/Ohana/src/opihi/cmd.data/tvcolors.c

-              r41340
+              r41341
     gprint (GP_ERR, "USAGE: tvcolors (colormap) [-nan red green blue]\n");
     gprint (GP_ERR, " colormap options : greyscale, -greyscale, rainbow, heat, fullcolor, ruffcolor (also grayscale, -grayscale)\n");
+    gprint (GP_ERR, " colormap file name of the for WORD:PATH where WORD = (file,csvf,lgcy,cetf)\n");
+    gprint (GP_ERR, " lgcy:\n");
+    gprint (GP_ERR, "   the colormap file must contain 4 columns: f R B G; each line defines a color transition.\n");
+    gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
+    gprint (GP_ERR, "   R,B,G: 0 - 1 define the value of the color at the transition point\n");
+    gprint (GP_ERR, " file:\n");
+    gprint (GP_ERR, " colormap file: if colormap name is given as file:path/to/file, the colormap is read from the file\n");
     gprint (GP_ERR, "   the colormap file must contain 4 columns: f R G B; each line defines a color transition.\n");
     gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
     gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
-    gprint (GP_ERR, " csvf:\n");
-    gprint (GP_ERR, "   the colormap file must contain 4 columns: f,R,G,B separate by commas; each line defines a color transition.\n");
-    gprint (GP_ERR, "   f: 0 - 1 defines the scale value for the transition\n");
-    gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
-    gprint (GP_ERR, " cetf:\n");
-    gprint (GP_ERR, "   the colormap file must 256 rows contain 3 columns: R,G,B separate by commas; each line defines a color transition.\n");
-    gprint (GP_ERR, "   R,G,B: 0 - 1 define the value of the color at the transition point\n");
-    gprint (GP_ERR, "   CETF-format files can be found at : https://peterkovesi.com/projects/colourmaps\n");
     return (FALSE);
+  }

trunk/Ohana/src/opihi/cmd.data/vgauss.c

-              r41340
+              r41341
 /* local private functions */
 opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
+int vgauss_apply (int argc, char **argv);
 # define GET_VAR(V,A) \
 …
   Vector *xvec, *yvec, *svec, *ovec;
   char *c, name[16];
+  if ((N = get_argument (argc, argv, "-apply"))) {
+    remove_argument (N, &argc, argv);
+    int status = vgauss_apply (argc, argv);
+    return status;
+  }
   Quiet = FALSE;
 …
+  }
   if (!Quiet) gprint (GP_ERR, "%d iterations\n", i);
+  set_variable ("Chisq", chisq);
+  set_variable ("ChisqNu", chisq / (float) (Npts - Npar));
   for (i = 0; i < Npts; i++) {
 …
+}
+int vgauss_apply (int argc, char **argv) {
+  opihi_flt par[4];
+  Vector *xvec, *ovec;
+  char *c;
+  // NOTE: -apply has already been stripped by vgauss call
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: vgauss -apply <x> <out>\n");
+    return (FALSE);
+  }
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  CastVector (xvec, OPIHI_FLT);
+  GET_VAR (par[0], "C0");
+  GET_VAR (par[1], "C1");
+  GET_VAR (par[2], "C2");
+  GET_VAR (par[3], "C3");
+  int Npar = 4;
+  int Npts = xvec[0].Nelements;
+  ResetVector (ovec, OPIHI_FLT, Npts);
+  for (int i = 0; i < Npts; i++) {
+    ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, NULL);
+  }
+  ovec[0].Nelements = Npts;
+  return TRUE;
+}
 /* pars: x_o, sigma, I, back */
 opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
 …
   f = par[2]*r + par[3];
+  dpar[0] = par[2]*r*z/par[1];
+  dpar[1] = par[2]*r*z*z/par[1];
+  dpar[2] = r;
+  dpar[3] = 1;
+  if (dpar) {
+    dpar[0] = par[2]*r*z/par[1];
+    dpar[1] = par[2]*r*z*z/par[1];
+    dpar[2] = r;
+    dpar[3] = 1;
+  }
   return (f);

trunk/Ohana/src/opihi/cmd.data/write_vectors.c

-              r41340
+              r41341
+  }
+  /* option generate a FITS output table */
+  // option generate a FITS output table (FITS holds the filename)
+  // in FITS output context, -header is interpretted as a buffer containing
+  // header keywords to supplement the FITS table header
   char *FITS = NULL;
+  Header *fitsheader = NULL;
+  Buffer *headbuffer = NULL;
   if ((N = get_argument (argc, argv, "-fits"))) {
     remove_argument (N, &argc, argv);
     FITS = strcreate (argv[N]);
     remove_argument (N, &argc, argv);
+  }
+    if ((N = get_argument (argc, argv, "-header"))) {
+      remove_argument (N, &argc, argv);
+      if ((headbuffer = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+      fitsheader = &headbuffer->header;
+      remove_argument (N, &argc, argv);
+    }
+  }
   /* option generate a FITS output table */
   int CSV = FALSE;
 …
   if (FITS) {
     int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile);
+    int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile);
     free (vec);
     return status;
 …
   if (ADD_HEADER) {
     for (j = 0; j < Nvec; j++) {
+      if (j == 0) fprintf (f, "# ");
       if (CSV) {
         fprintf (f, "%s,", vec[j][0].name);
       } else {
-        if (j == 0) fprintf (f, "# ");
         fprintf (f, "%s ", vec[j][0].name);
+      }
 …
   /* default output format */
   if (format == (char *) NULL) {
+    char padChar = CSV ? ',' : ' ';
     for (i = 0; i < vec[0][0].Nelements; i++) {
       for (j = 0; j < Nvec; j++) {
+        if (vec[j][0].type == OPIHI_FLT) {
+          if (CSV) {
+            fprintf (f, "%.12g,", vec[j][0].elements.Flt[i]);
+          } else {
+            fprintf (f, "%.12g ", vec[j][0].elements.Flt[i]);
+          }
+        } else {
+          if (CSV) {
+            fprintf (f, OPIHI_INT_FMT",", vec[j][0].elements.Int[i]);
+          } else {
+            fprintf (f, OPIHI_INT_FMT" ", vec[j][0].elements.Int[i]);
+          }
+        switch (vec[j][0].type) {
+          case OPIHI_FLT:
+            fprintf (f, "%.12g%c", vec[j][0].elements.Flt[i], padChar);
+            break;
+          case OPIHI_INT:
+            fprintf (f, OPIHI_INT_FMT"%c", vec[j][0].elements.Int[i], padChar);
+            break;
+          case OPIHI_STR:
+            fprintf (f, "%s%c", vec[j][0].elements.Str[i], padChar);
+            break;
+        }
+      }
 …
         fmttype[j] = 'd';
         break;
+      case 's':
+        fmttype[j] = 's';
+        break;
       default:
         gprint (GP_ERR, "syntax error in format (only e,f,d,c,x allowed)\n");
+        gprint (GP_ERR, "syntax error in format (only e,f,d,c,x,s allowed)\n");
         return (FALSE);
+    }
 …
   strcat (fmtlist[Nvec-1], p0);
+  // check format types against vector types:
+  for (j = 0; j < Nvec; j++) {
+    switch (vec[j][0].type) {
+      case OPIHI_FLT:
+      case OPIHI_INT:
+        if (fmttype[j] == 's') {
+          gprint (GP_ERR, "mismatch between string format and numerical vector for %s\n", vec[j][0].name);
+          return FALSE;
+        }
+        break;
+      case OPIHI_STR:
+        if (fmttype[j] != 's') {
+          gprint (GP_ERR, "mismatch between numerical format and string vector for %s\n", vec[j][0].name);
+          return FALSE;
+        }
+        break;
+    }
+  }
   for (i = 0; i < vec[0][0].Nelements; i++) {
     for (j = 0; j < Nvec; j++) {
 …
+        }
+      }
+      if (fmttype[j] == 's') {
+        fprintf (f, fmtlist[j], vec[j][0].elements.Str[i]);
+      }
+    }
     fprintf (f, "\n");
 …
     gprint (GP_ERR, "  -append : write to the end of the existing file\n");
     gprint (GP_ERR, "  -fits NAME : write a fits table (extention name is NAME, column names match vector names)\n");
+    gprint (GP_ERR, "     in FITS output context, -header takes an additional argument which is interpretted\n");
+    gprint (GP_ERR, "     as a buffer containing header keywords to supplement the FITS table header\n");
     gprint (GP_ERR, "  -csv : write a comma-separated values file (eg, to read in excel)\n");
     gprint (GP_ERR, "  -f \"format\" : provide formatting codes for output:\n");

trunk/Ohana/src/opihi/dvo/Makefile

r41340	r41341
67	67	$(SRC)/gstar.$(ARCH).o \
68	68	$(SRC)/hosts.$(ARCH).o \
	69	$(SRC)/cimages.$(ARCH).o \
69	70	$(SRC)/images.$(ARCH).o \
70	71	$(SRC)/imbox.$(ARCH).o \

trunk/Ohana/src/opihi/dvo/avextract.c

r41340	r41341
272	272	// write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
273	273	if (RESULT_FILE && !SKIP_RESULTS) {
274		int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, FALSE, NULL, 0);
	274	int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nreturn, FALSE, FALSE, NULL, 0);
275	275	if (!status) {
276	276	goto escape;

trunk/Ohana/src/opihi/dvo/avmatch.c

-              r41340
+              r41341
       CoordsFile = abspath("coords.fits", 1024);
       int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
+      int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
       if (!status) goto escape;
+    }
 …
       vec[i][0].Nelements = Nfound;
+    }
     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields + 1, FALSE, FALSE, NULL, 0);
+    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nfields + 1, FALSE, FALSE, NULL, 0);
     free (vec[Nfields]->elements.Int);
     free (vec[Nfields]);

trunk/Ohana/src/opihi/dvo/avperiodomatch.c

r41340	r41341
99	99
100	100	CoordsFile = abspath("coords.fits", 1024);
101		int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
	101	int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
102	102	if (!status) goto escape;
103	103	free (vec);

trunk/Ohana/src/opihi/dvo/avselect.c

r41340	r41341
132	132
133	133	CoordsFile = abspath("coords.fits", 1024);
134		int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
	134	int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
135	135	if (!status) goto escape;
136	136	}
…	…
327	327	vec[NfieldsOut-1] = RADvec;
328	328	}
329		int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, NfieldsOut, FALSE, FALSE, NULL, 0);
	329	int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, NfieldsOut, FALSE, FALSE, NULL, 0);
330	330	if (!status) goto escape;
331	331	}

trunk/Ohana/src/opihi/dvo/dvo_host_utils.c

r41340	r41341
233	233	// write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
234	234	if (ResultFile) {
235		int status = WriteVectorTableFITS (ResultFile, "RESULT", vec, Nvec, FALSE, FALSE, NULL, 0);
	235	int status = WriteVectorTableFITS (ResultFile, "RESULT", NULL, vec, Nvec, FALSE, FALSE, NULL, 0);
236	236	if (!status) {
237	237	gprint (GP_ERR, "failed to write result file %s\n", ResultFile);

trunk/Ohana/src/opihi/dvo/images.c

-              r41340
+              r41341
   off_t i, Nimage, Nmosaic;
+  int j, status, InPic, leftside, *plist, TimeSelect, ByName;
+  int j, status, InPic, leftside, TimeSelect, ByName;
+  // int *plist, n, npts; -- what is/was plist?
   int WITH_MOSAIC, SOLO_MOSAIC;
   time_t tzero, tend;
   int N, NPTS, n, npts, Npts, kapa, *foundMosaic;
+  int n, N, NPTS, Npts, kapa, *foundMosaic;
   Vector Xvec, Yvec;
   double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius;
 …
   int TrianglePts  = wordhash ("TRI-");
+  npts = NPTS = 200;
+  // npts = 200
+  NPTS = 200;
   SetVector (&Xvec, OPIHI_FLT, NPTS);
   SetVector (&Yvec, OPIHI_FLT, NPTS);
 …
   Image *mosaic = NULL;
+  ALLOCATE (plist, int, NPTS);
+  n = N = 0;
+  /******** stage 1 : find the images to plot *********/
+  // ALLOCATE (plist, int, NPTS);
+  n = 0; // tracks number of images
+  N = 0; // tracks number of vertices
   for (i = 0; i < Nimage; i++) {
     if (ByName && strncmp (image[i].name, name, strlen(name))) continue;
 …
     if (!InPic) continue;
+    plist[n] = i;
+    n++;
+    if (n > npts - 1) {
+      npts += 200;
+      REALLOCATE (plist, int, npts);
+    }
+    N+=2*Npts;
+    // plist[n] = i;
+    // n++;
+    // if (n > npts - 1) {
+    //   npts += 200;
+    //   REALLOCATE (plist, int, npts);
+    // }
+    n ++;
+    N += 2*Npts;
     if (N + 16 >= NPTS) {  /* need to leave room for 8 point image */
       NPTS += 400;

trunk/Ohana/src/opihi/dvo/init.c

-              r41340
+              r41341
 int catdir_define   PROTO((int, char **));
 int ccd             PROTO((int, char **));
+int cimages         PROTO((int, char **));
 int cmatch          PROTO((int, char **));
 int cmd             PROTO((int, char **));
 …
   {1, "catdir",      catdir_define,"re-define CATDIR"},
 //  {1, "ccd",         ccd,          "plot color-color diagram"},
+  {1, "cimages",     cimages,      "fill image boxes with a color"},
   {1, "cmatch",      cmatch,       "match two catalogs"},
 //  {1, "cmd",         cmd,          "plot cmd of stars in current region"},

trunk/Ohana/src/opihi/dvo/mextract.c

r41340	r41341
352	352	// write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
353	353	if (RESULT_FILE && !SKIP_RESULTS) {
354		int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, FALSE, NULL, 0);
	354	int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nreturn, FALSE, FALSE, NULL, 0);
355	355	if (!status) goto escape;
356	356	}

trunk/Ohana/src/opihi/dvo/mmatch.c

r41340	r41341
148	148	// XXX this is now set for both cases...
149	149	CoordsFile = abspath("coords.fits", 1024);
150		int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
	150	int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
151	151	if (!status) goto escape;
152	152	}
…	…
349	349	vec[Nfields-1] = IDXvec;
350	350	}
351		int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields, FALSE, FALSE, NULL, 0);
	351	int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nfields, FALSE, FALSE, NULL, 0);
352	352	if (!status) goto escape;
353	353	}

trunk/Ohana/src/opihi/dvo/skyregion.c

-              r41340
+              r41341
   int N;
+  // alternate call form: skyregion -box Rcenter Dcenter radius
+  if ((N = get_argument (argc, argv, "-box"))) {
+    if ((N != 1) && (argc != 5)) {
+      gprint (GP_ERR, "USAGE: skyregion -box (Rcenter) (Dcenter) (radius)\n");
+      gprint (GP_ERR, "  Rcenter, Dcenter, radius in decimal degrees\n");
+      return (FALSE);
+    }
+    // argv[1] = -box
+    double Rcenter = atof(argv[2]);
+    double Dcenter = atof(argv[3]);
+    double Radius  = atof(argv[4]);
+    double Rmin = Rcenter - Radius/cos(DEG_RAD*Dcenter);
+    double Rmax = Rcenter + Radius/cos(DEG_RAD*Dcenter);
+    double Dmin = Dcenter - Radius;
+    double Dmax = Dcenter + Radius;
+    set_skyregion (Rmin, Rmax, Dmin, Dmax);
+    return TRUE;
+  }
   // dvo_client should have 2 standard arguments: -hostID and -hostdir
 …
       gprint (GP_ERR, "current skyregion: %f - %f : %f - %f\n", Rmin, Rmax, Dmin, Dmax);
       gprint (GP_ERR, "USAGE:  skyregion (min RA) (max RA) (min DEC) (max DEC)\n");
+      gprint (GP_ERR, "        skyregion -box (RA) (DEC) (radius) [sets Rmin,Rmax & Dmin,Dmax]\n");
       return FALSE;
+    }

trunk/Ohana/src/opihi/include/data.h

-              r41340
+              r41341
   int Nx;
   int Ny;
+  float **buffers;
+  float **flx;
+  float **var;
 } MedImageType;

trunk/Ohana/src/opihi/include/deimos.h

-              r41340
+              r41341
 } DeimosResult;
 float     *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row);
+float     *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row);
 void       deimos_make_kernel (float stilt, int Nx);
 void       deimos_free_kernel ();
 …
 // internal to make_model:
 float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row);
+float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row);
 float     *deimos_apply_tilt (float *input, int Nx, int Ny);
 void       deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny);
 void       deimos_add_background (opihi_flt *backgnd, float *out, int Nx, int Ny, int row);
 float     *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row);
+float     *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row);

trunk/Ohana/src/opihi/include/dvomath.h

-              r41340
+              r41341
 # define OPIHI_NAME_SIZE 1024
+/* OPIHI_FLT, OPIHI_INT and related are defined in libdvo/include/dvodb.h */
 # define NCHARS 256
-// # define opihi_flt double
-// # define opihi_int int
-// NOTE: if opihi_int is changed to unsigned, all subtraction and negation operations
-// need to result in a float value (or 3 - 5 will yield the unexpected value 2^32 - 2)
 # define REQUIRE_VECTOR_FLT(VECT,RVAL) { \
 …
 enum {ANYVECTOR, NEWVECTOR, OLDVECTOR};
 enum {ANYBUFFER, NEWBUFFER, OLDBUFFER};
-#ifdef NOT_MOVED_TO_LIBDVO
-enum {OPIHI_FLT, OPIHI_INT};
-#endif
 typedef struct {                        /* representation of a variable (0-D) */
 …
   union {
     void      *Ptr;
+    char     **Str;
     opihi_flt *Flt;
     opihi_int *Int;
 …
 /* vector IO functions */
 int           WriteVectorTableFITS  PROTO((char *filename, char *extname, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile));
+int           WriteVectorTableFITS  PROTO((char *filename, char *extname, Header *extraheader, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile));
 Vector      **ReadVectorTableFITS   PROTO((char *filename, char *extname, int *Nvec));

trunk/Ohana/src/opihi/lib.data/MedImageOps.c

-              r41340
+              r41341
   free (medimage[0].name);
   for (i = 0; i < medimage[0].Ninput; i++) {
+    free (medimage[0].buffers[i]);
+    free (medimage[0].flx[i]);
+    FREE (medimage[0].var[i]);
+  }
+  free (medimage[0].buffers);
+  free (medimage[0].flx);
+  free (medimage[0].var);
   free (medimage);
+}
 …
   medimage->Nx = Nx;
   medimage->Ny = Ny;
+  ALLOCATE (medimage->buffers, float *, 1);
+  ALLOCATE (medimage->flx, float *, 1);
+  ALLOCATE (medimage->var, float *, 1);
   medimages[N] = medimage;

trunk/Ohana/src/opihi/lib.data/SplineOps.c

-              r41340
+              r41341
     ALLOCATE (spline[0].y2, opihi_flt, spline[0].Nknots);
+  }
+  memset (spline[0].xk, 0, spline[0].Nknots * sizeof(opihi_flt));
+  memset (spline[0].yk, 0, spline[0].Nknots * sizeof(opihi_flt));
+  memset (spline[0].y2, 0, spline[0].Nknots * sizeof(opihi_flt));
+  if (spline[0].Nknots) {
+    memset (spline[0].xk, 0, spline[0].Nknots * sizeof(opihi_flt));
+    memset (spline[0].yk, 0, spline[0].Nknots * sizeof(opihi_flt));
+    memset (spline[0].y2, 0, spline[0].Nknots * sizeof(opihi_flt));
+  }
+}
 …
   // NOTE: if we want to compress the output table, use native byte order here (last element)
   gfits_set_bintable_column_reformat (&theader, &ftable, "X_KNOT", "double", myspline->xk, myspline->Nknots, FALSE);
   gfits_set_bintable_column_reformat (&theader, &ftable, "Y_KNOT", "double", myspline->yk, myspline->Nknots, FALSE);
   gfits_set_bintable_column_reformat (&theader, &ftable, "DY2_DX", "double", myspline->y2, myspline->Nknots, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "X_KNOT", "double", myspline->xk, myspline->Nknots, 0, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "Y_KNOT", "double", myspline->yk, myspline->Nknots, 0, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "DY2_DX", "double", myspline->y2, myspline->Nknots, 0, FALSE);
   if (!append) {

trunk/Ohana/src/opihi/lib.data/open_kapa.c

-              r41340
+              r41341
 int close_kapa (char *name) {
+  struct timespec request, remain;
   int N;
 …
+  }
   DelKapaDevice (name);
+  // avoid blocking on waitpid, test every 1 msec, up to 500 msec
+  request.tv_sec = 0;
+  request.tv_nsec = 1000000;
+  // try to harvest the child PID
+  int waitstatus = 0;
+  int result = waitpid (-1, &waitstatus, WNOHANG);
+  for (int i = 0; (i < 500) && (result == 0); i++) {
+    nanosleep (&request, &remain);
+    result = waitpid (-1, &waitstatus, WNOHANG);
+  }
+  if ((result == -1) && (errno != ECHILD)) {
+    fprintf (stderr, "unexpected error from waitpid (%d): programming error\n", errno);
+  }
+  if (result == 0) {
+    fprintf (stderr, "child did not exit (close_kapa), timeout");
+  }
   return (TRUE);
+}

trunk/Ohana/src/opihi/lib.data/spline.c

-              r41340
+              r41341
   float dy, dx, *tmp;
+  // spline is not valid with < 3 points
+  if (N < 3) return;
   ALLOCATE (tmp, float, N);
 …
   opihi_flt dy, dx, *tmp;
+  // spline is not valid with < 3 points
+  if (N < 3) return;
   ALLOCATE (tmp, opihi_flt, N);

trunk/Ohana/src/opihi/lib.data/starfuncs.c

r41340	r41341
75	75	FWHMx = 2.355sqrt (fabs(x2 / I - xx));
76	76	FWHMy = 2.355sqrt (fabs(y2 / I - yy));
	77
	78	// fprintf (stderr, "Mxx, Myy: %f, %f\n", x2/I - xx, y2/I - yy);
77	79	Sxy = xy / I - x*y;
78	80	mag = -2.5*log10(I);

trunk/Ohana/src/opihi/lib.data/style_args.c

-              r41340
+              r41341
+  }
+  if ((graphmode[0].style == KAPA_PLOT_POLYGON) || (graphmode[0].style == KAPA_PLOT_POLYFILL)) {
+    if ((N = get_argument (*argc, argv, "-npoint"))) {
+      remove_argument (N, argc, argv);
+      graphmode[0].ptype = atoi (argv[N]);
+      remove_argument (N, argc, argv);
+    } else {
+      gprint (GP_ERR, "polygon & polyfill styles require number of points argument: -npoint N\n");
+      return FALSE;
+    }
+  }
   return (TRUE);
+}

trunk/Ohana/src/opihi/lib.shell/VectorIO.c

-              r41340
+              r41341
 void gfits_compress_timing ();
+static int VectorGetMaxStringLength (Vector *vec) {
+  int MaxLength = 0;
+  if (vec[0].type != OPIHI_STR) return 0;
+  for (int i = 0; i < vec[0].Nelements; i++) {
+    MaxLength = MAX (MaxLength, strlen(vec[0].elements.Str[i]));
+  }
+  return MaxLength;
+}
 // write a set of vectors to a FITS FTable structure (vectors names become fits column names)
 static int WriteVectorTable (FTable *ftable, char *extname, Vector **vec, int Nvec, char *format, char nativeOrder) {
 …
   int j;
-  char *tformat = NULL;
   Header *theader = ftable->header;
   gfits_create_table_header (theader, "BINTABLE", extname);
+  ALLOCATE (tformat, char, 2*Nvec);
+  // allocate an array of strings to represent the format for each output field
+  // formats include single column formats (BIJKDE) and multi-column formats (e.g., 2I)
+  // we will have no more than Nvec fields (but we can have fewer)
+  int Nfield = 0;
+  ALLOCATE_PTR (tformat, char *, Nvec);
+  ALLOCATE_PTR (Nelement, int, Nvec);
   if (format) {
+    // the bintable format string can only define the byte-width of each field.  valid output columns are currently:
+    // the bintable format string can defines the byte-width of each field and number of elements (columns per field).
+    // valid output columns are currently:
     // B (char), I (16 bit short), J (32 bit int), E (32 bit float), D (64 bit double).
+    // the format string is just the sequence of types, eg: LIIJEED
+    // validate the format string
+    // the format string is just the sequence of types, eg: LIIJEED.
+    // it may have spaces or integer element counts:
+    // "2D 4I EEJ"
+    // *** validate the format string
+    // as I parse each elements, if it is a digit, I need parse that value
     char *ptr = format;
     for (j = 0; j < Nvec; j++) {
+    for (j = 0; j < Nvec; ) {
       while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++;
       if (*ptr == 0) {
 …
         goto escape;
+      }
+      if ((*ptr != 'B') && (*ptr != 'I') && (*ptr != 'J') && (*ptr != 'K') && (*ptr != 'D') && (*ptr != 'E')) {
+      // is there a leading integer?
+      char *endptr;
+      Nelement[Nfield] = strtol (ptr, &endptr, 10);
+      if (endptr == ptr) {
+        Nelement[Nfield] = 1;
+      }
+      ptr = endptr; // this should now point at the letter that is the format type
+      if ((*ptr != 'B') && (*ptr != 'I') && (*ptr != 'J') && (*ptr != 'K') && (*ptr != 'D') && (*ptr != 'E') && (*ptr != 'A')) {
         gprint (GP_ERR, "error in binary table format %s: invalid character %c\n", format, *ptr);
         goto escape;
+      }
+      tformat[2*j + 0] = *ptr;
+      tformat[2*j + 1] = 0; // a bit sleazy : use a 2xN string to store N 1-byte strings
+      int Nchar = snprintf (tformat[Nfield], 0, "%d%c", Nelement[Nfield], *ptr);
+      ALLOCATE (tformat[Nfield], char, Nchar + 1);
+      int Nout = snprintf (tformat[Nfield], Nchar + 1, "%d%c", Nelement[Nfield], *ptr);
+      myAssert (Nout <= Nchar, "oops");
+      // tformat[2*j + 0] = *ptr;
+      // tformat[2*j + 1] = 0; // a bit sleazy : use a 2xN string to store N 1-byte strings
+      // For numeric multi-valued fields, the number of elements from the format defines the number of
+      // vectors which go into that field.  For string-type vectors, the format specifies the maximum number of characters
+      // that go into the field.
+      // for example, a format code of 3E should match a list of three numeric-type vectors while a format code of 15A should
+      // match a single string-type vector which will supply up to 15 chars per row.
+      if (*ptr != 'A') {
+        j += Nelement[Nfield]; // advance past Nelement vectors
+      } else {
+        if (vec[j][0].type != OPIHI_STR) {
+          gprint (GP_ERR, "error in binary table format %s (mismatch between string format and numeric vector %s)\n", format, vec[j][0].name);
+          goto escape;
+        }
+        j ++; // advance past a single string-type vector (validate that?)
+      }
+      if (j > Nvec) {
+        gprint (GP_ERR, "error in binary table format %s (too few vectors for listed field)\n", format);
+        goto escape;
+      }
       ptr ++;
+      Nfield ++;
+    }
     while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++;
 …
   } else {
     for (j = 0; j < Nvec; j++) {
+      // if the format is not defined, just use the native byte-widths
+      tformat[2*j + 0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'K'; // this depends on opihi_int == int64_t for Int
+      tformat[2*j + 1] = 0;
+    }
+      switch (vec[j][0].type) {
+        case OPIHI_FLT:
+        case OPIHI_INT:
+          // if the format is not defined, just use the native byte-widths
+          ALLOCATE (tformat[j], char, 2);
+          tformat[j][0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'K';
+          tformat[j][1] = 0;
+          Nelement[j] = 1;
+          break;
+        case OPIHI_STR:
+          // we need to examine the vector to determine the length
+          Nelement[j] = VectorGetMaxStringLength(vec[j]);
+          int Nchar = snprintf (tformat[j], 0, "%d%c", Nelement[j], 'A');
+          ALLOCATE (tformat[j], char, Nchar + 1);
+          int Nout = snprintf (tformat[j], Nchar + 1, "%d%c", Nelement[j], 'A');
+          myAssert (Nout <= Nchar, "oops");
+          break;
+      }
+    }
+    Nfield = Nvec;
+  }
 …
   // output table (because the data goes to the named column below).  need to enforce
   // this somehow
+  for (j = 0; j < Nvec; j++) {
+    gfits_define_bintable_column (theader, &tformat[2*j], vec[j][0].name, NULL, NULL, 1.0, 0.0);
+  int ivec = 0;
+  for (j = 0; j < Nfield; j++) {
+    // XX need to loop over fields, and skip the additional vectors that are part of a field
+    // this call supported multiple columns per named field
+    gfits_define_bintable_column (theader, tformat[j], vec[ivec][0].name, NULL, NULL, 1.0, 0.0);
+    if (vec[ivec][0].type == OPIHI_STR) {
+      ivec ++;
+    } else {
+      ivec += Nelement[j];
+    }
+  }
+  // need to free the array
+  for (j = 0; j < Nfield; j++) {
+    free (tformat[j]);
+  }
   free (tformat);
 …
   gfits_create_table (theader, ftable);
+  // I need to add each vector in order, but I need to
+  // track which field it corresponds to.
   // add the vectors to the output array
+  for (j = 0; j < Nvec; j++) {
+    if (vec[j][0].type == OPIHI_FLT) {
+      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "double",  vec[j][0].elements.Flt, vec[j][0].Nelements, nativeOrder);
+  for (ivec = 0, j = 0; j < Nfield; j++) {
+    // the first vector provides the name for the field
+    char *fieldname = vec[ivec][0].name;
+    if (vec[ivec][0].type == OPIHI_STR) {
+      // string-type vectors need to be copied into a contiguous buffer with the right dimensions:
+      Vector *thisvec = vec[ivec];
+      ALLOCATE_PTR (strbuffer, char, Nelement[j]*thisvec->Nelements);
+      for (int i = 0; i < thisvec->Nelements; i++) {
+        int nChar = MIN (strlen (thisvec->elements.Str[i]), Nelement[j]);
+        // fprintf (stderr, "%d %d %d : %d : %d : %s\n", ivec, j, i, Nelement[j], nChar, thisvec->elements.Str[i]);
+        memcpy (&strbuffer[i*Nelement[j]], thisvec->elements.Str[i], nChar);
+      }
+      gfits_set_bintable_column (theader, ftable, fieldname, strbuffer, thisvec->Nelements);
+      free (strbuffer);
+      ivec ++;
     } else {
+//    gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int",     vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
+      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int64_t", vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
+      for (int k = 0; k < Nelement[j]; k++, ivec++) {
+        Vector *thisvec = vec[ivec];
+        switch (thisvec->type) {
+          case OPIHI_FLT:
+            gfits_set_bintable_column_reformat (theader, ftable, fieldname, "double",  thisvec->elements.Flt, thisvec->Nelements, k, nativeOrder);
+            break;
+          case OPIHI_INT:
+            gfits_set_bintable_column_reformat (theader, ftable, fieldname, "int64_t", thisvec->elements.Int, thisvec->Nelements, k, nativeOrder);
+            break;
+        }
+      }
+    }
+  }
 …
+}
+// insert a new header line, return end pointer
+char *gfits_insert_header_line (Header *header, char *endptr, char *newline) {
+  if (0) {
+    char temp1[32], temp2[32];
+    memset (temp1, 0, 32);
+    memset (temp2, 0, 32);
+    if (endptr) {
+      memcpy (temp1, endptr, 8);
+    }
+    memcpy (temp2, newline, 8);
+    fprintf (stderr, "insert: %s : %s\n", temp1, temp2);
+  }
+  /* find the END of the header, if not supplied */
+  if (!endptr) {
+    endptr = gfits_header_field (header, "END", 1);
+    if (endptr == NULL) return NULL;
+  }
+  /* is there enough space for 1 more line? */
+  if (header[0].datasize - (endptr - (header[0].buffer)) < 2*FT_LINE_LENGTH) {
+    // no, so expand by a full header block (FT_RECORD_SIZE = 32*80 = 2880)
+    header[0].datasize += FT_RECORD_SIZE;
+    REALLOCATE (header[0].buffer, char, header[0].datasize);
+    // re-find the "END" marker, in case new memory block is used
+    endptr = gfits_header_field (header, "END", 1);
+    if (endptr == NULL) return NULL;
+    memset (endptr + FT_LINE_LENGTH, ' ', FT_RECORD_SIZE);
+  }
+  /* push END line back 1 */
+  memmove ((endptr + FT_LINE_LENGTH), endptr, FT_LINE_LENGTH);
+  memset (endptr, ' ', FT_LINE_LENGTH);
+  strncpy (endptr, newline, 80);
+  endptr += FT_LINE_LENGTH;
+  return endptr;
+}
+static char *rawkeywords[] = {"SIMPLE", "BITPIX", "NAXIS", "PCOUNT", "GCOUNT", "EXTEND", "EXTNAME", "BSCALE", "BZERO", "TFIELDS", "TFORM", "TZERO", "TSCAL", "        ", NULL};
 // write a set of vectors to a FITS file (vectors names become fits column names)
 int WriteVectorTableFITS (char *filename, char *extname, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile) {
+int WriteVectorTableFITS (char *filename, char *extname, Header *extraheader, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile) {
   Header rawheader;
 …
   rawtable.header = &rawheader;
   if (!WriteVectorTable (&rawtable, extname, vec, Nvec, format, (compress != NULL))) goto escape;
   // NOTE : for compression, the table is constructed in native by order
+  // NOTE : for compression, the table is constructed in native byte order
   FTable *outtable = &rawtable;
 …
   if (!append) {
+    // generate a blank PHU header
     Header header;
     Matrix matrix;
 …
+  }
+  if (extraheader) {
+    // copy keywords which are not the standard or table keywords
+    char *buf = extraheader->buffer;
+    char *endptr = NULL;
+    for (int i = 0; i < extraheader->datasize; i+= FT_LINE_LENGTH, buf += FT_LINE_LENGTH) {
+      if (0) {
+        char temp1[32];
+        memset (temp1, 0, 32);
+        memcpy (temp1, buf, 8);
+        fprintf (stderr, "buffer: %s\n", temp1);
+      }
+      for (int j = 0; rawkeywords[j] != NULL; j++) {
+        if (!strncmp (buf, rawkeywords[j], strlen(rawkeywords[j]))) goto skip_insert;
+      }
+      endptr = gfits_insert_header_line (outheader, endptr, buf);
+      if (!endptr) {
+        gprint (GP_ERR, "failed to update FITS header with extra keywords\n");
+        return (FALSE);
+      }
+    skip_insert:
+      continue;
+    }
+  }
   // write the actual table data
   gfits_fwrite_Theader (f, outheader);

trunk/Ohana/src/opihi/lib.shell/VectorOps.c

-              r41340
+              r41341
   if (mode == NEWVECTOR) goto error;
   if (vectors[i]) {
+    // XXX warning: this will be a leak if vector is type OPIHI_STR
     if (vectors[i][0].elements.Ptr) free (vectors[i][0].elements.Ptr);
     free (vectors[i]);
 …
   if (i == Nvectors) return (FALSE);
+  // XXX warning: this will be a leak if vector is type OPIHI_STR
   if (vectors[i][0].elements.Ptr) free (vectors[i][0].elements.Ptr);
   free (vectors[i]);
 …
   if (i == Nvectors) return (FALSE);
+  // XXX warning: this will be a leak if vector is type OPIHI_STR
   if (vectors[i][0].elements.Ptr) free (vectors[i][0].elements.Ptr);
   free (vectors[i]);
 …
 int CopyVector (Vector *out, Vector *in) {
+  // XXX warning: this will be a leak if vector is type OPIHI_STR
   if (out[0].elements.Ptr) free (out[0].elements.Ptr);
   out[0].Nelements = in[0].Nelements;
 …
 int MatchVector(Vector *out, Vector *in, char type) {
+  // XXX warning: this will be a leak if vector is type OPIHI_STR
   if (out[0].elements.Ptr) free (out[0].elements.Ptr);
   out[0].Nelements = in[0].Nelements;
 …
 int ResetVector (Vector *vec, char type, int Nelements) {
+  // if the supplied vector is a string but the output is not a string, we need to free
+  // the unused elements
+  if ((vec[0].type == OPIHI_STR) && (vec[0].type != type)) {
+    for (int i = 0; i < vec[0].Nelements; i++) {
+      FREE (vec[0].elements.Str[i]);
+    }
+  }
   // a vector can only have >= 0 elements
   vec[0].Nelements = MAX(Nelements,0);
+  if (type == OPIHI_FLT) {
+    REALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1, Nelements));
+    vec[0].type = OPIHI_FLT;
+  } else {
+    REALLOCATE (vec[0].elements.Int, opihi_int, MAX(1, Nelements));
+    vec[0].type = OPIHI_INT;
+  switch (type) {
+    case OPIHI_FLT:
+      REALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1, Nelements));
+      vec[0].type = OPIHI_FLT;
+      break;
+    case OPIHI_INT:
+      REALLOCATE (vec[0].elements.Int, opihi_int, MAX(1, Nelements));
+      vec[0].type = OPIHI_INT;
+      break;
+    case OPIHI_STR:
+      REALLOCATE (vec[0].elements.Str, char *, MAX(1, Nelements));
+      vec[0].type = OPIHI_STR;
+      break;
+  }
   return TRUE;
 …
 // SetVector (vecx, OPIHI_FLT, MAX (Npts, 1));
+// Use this for an unallocated vector (e.g., static variable)
 int SetVector (Vector *vec, char type, int Nelements) {
   vec[0].Nelements = MAX(Nelements,0);
+  if (type == OPIHI_FLT) {
+    ALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1,Nelements));
+    vec[0].type = OPIHI_FLT;
+  } else {
+    ALLOCATE (vec[0].elements.Int, opihi_int, MAX(1,Nelements));
+    vec[0].type = OPIHI_INT;
+  switch (type) {
+    case OPIHI_FLT:
+      ALLOCATE (vec[0].elements.Flt, opihi_flt, MAX(1,Nelements));
+      vec[0].type = OPIHI_FLT;
+      break;
+    case OPIHI_INT:
+      ALLOCATE (vec[0].elements.Int, opihi_int, MAX(1,Nelements));
+      vec[0].type = OPIHI_INT;
+      break;
+    case OPIHI_STR:
+      ALLOCATE (vec[0].elements.Ptr, char *, MAX(1,Nelements));
+      vec[0].type = OPIHI_STR;
+      break;
+  }
   return TRUE;
 …
   if (vec[0].type == type) return TRUE;
+  if (type == OPIHI_FLT) {
+    opihi_flt *temp;
+    ALLOCATE (temp, opihi_flt, vec[0].Nelements);
+    opihi_flt *vo = temp;
+    opihi_int *vi = vec[0].elements.Int;
+    for (i = 0; i < vec[0].Nelements; i++, vo++, vi++) {
+      *vo = *vi;
+    }
+    free (vec[0].elements.Int);
+    vec[0].elements.Flt = temp;
+    vec[0].type = OPIHI_FLT;
+  } else {
+    opihi_int *temp;
+    ALLOCATE (temp, opihi_int, vec[0].Nelements);
+    opihi_int *vo = temp;
+    opihi_flt *vi = vec[0].elements.Flt;
+    for (i = 0; i < vec[0].Nelements; i++, vo++, vi++) {
+      *vo = *vi;
+    }
+    free (vec[0].elements.Flt);
+    vec[0].elements.Int = temp;
+    vec[0].type = OPIHI_INT;
+  switch (type) {
+    case OPIHI_FLT: {
+      opihi_flt *temp;
+      ALLOCATE (temp, opihi_flt, vec[0].Nelements);
+      opihi_flt *vo = temp;
+      opihi_int *vi = vec[0].elements.Int;
+      for (i = 0; i < vec[0].Nelements; i++, vo++, vi++) {
+        *vo = *vi;
+      }
+      free (vec[0].elements.Int);
+      vec[0].elements.Flt = temp;
+      vec[0].type = OPIHI_FLT;
+      break;
+    }
+    case OPIHI_INT: {
+      opihi_int *temp;
+      ALLOCATE (temp, opihi_int, vec[0].Nelements);
+      opihi_int *vo = temp;
+      opihi_flt *vi = vec[0].elements.Flt;
+      for (i = 0; i < vec[0].Nelements; i++, vo++, vi++) {
+        *vo = *vi;
+      }
+      free (vec[0].elements.Flt);
+      vec[0].elements.Int = temp;
+      vec[0].type = OPIHI_INT;
+      break;
+    }
+    case OPIHI_STR:
+    default:
+      // does it make sense to cast an int/flt vector to string?
+      break;
+  }
   return TRUE;
 …
   gprint (GP_LOG, "    N       name      size\n");
   for (i = 0; i < Nvectors; i++) {
+    if (vectors[i][0].type == OPIHI_FLT) {
+      gprint (GP_LOG, "%5d %10s %10d (FLT)\n", i, vectors[i][0].name, vectors[i][0].Nelements);
+    } else {
+      gprint (GP_LOG, "%5d %10s %10d (INT)\n", i, vectors[i][0].name, vectors[i][0].Nelements);
+    switch (vectors[i][0].type) {
+      case OPIHI_FLT:
+        gprint (GP_LOG, "%5d %10s %10d (FLT)\n", i, vectors[i][0].name, vectors[i][0].Nelements);
+        break;
+      case OPIHI_INT:
+        gprint (GP_LOG, "%5d %10s %10d (INT)\n", i, vectors[i][0].name, vectors[i][0].Nelements);
+        break;
+      case OPIHI_STR:
+        gprint (GP_LOG, "%5d %10s %10d (STR)\n", i, vectors[i][0].name, vectors[i][0].Nelements);
+        break;
+    }
+  }

trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c

-              r41340
+              r41341
     if (!strcmp (argv[i], "atan2"))  { type = ST_BINARY; strcpy (argv[i], "a"); goto gotit; }
     if (!strcmp (argv[i], "datan2")) { type = ST_BINARY; strcpy (argv[i], "d"); goto gotit; }
+    if (!strcmp (argv[i], "<~"))     { type = ST_BINARY; strcpy (argv[i], "l"); goto gotit; }
+    if (!strcmp (argv[i], "~>"))     { type = ST_BINARY; strcpy (argv[i], "r"); goto gotit; }
     if (!strcmp (argv[i], ","))      { type = ST_COMMA; goto gotit; }
 …
     if (!strcmp (argv[i], ">>"))     { type = ST_LOGIC; strcpy (argv[i], "U"); goto gotit; }
     if (!strcmp (argv[i], "<<"))     { type = ST_LOGIC; strcpy (argv[i], "D"); goto gotit; }
+    /* XXX I would like to change the syntax to allow << and >> to mean bitshifts
+       but that means breaking these older values which means MIN and MAX */
+    // for now, use <- and -> to mean bitshift
     if (!strcmp (argv[i], "&&"))     { type = ST_AND; strcpy (argv[i], "A"); goto gotit; }

trunk/Ohana/src/opihi/lib.shell/multicommand.c

-              r41340
+              r41341
   p = line;
+  int Nline = strlen(line);
   done = FALSE;
   status = TRUE;
 …
     tmpline = strncreate (p, q - p);
     stripwhite (tmpline);
     if (*tmpline) {
+    myAssert (tmpline, "oops");
+      if (bufferPending) {
+        // flush old messages
+        ExpectMessage (server, 0.2, &message);
+        bufferPending = FALSE;
+    // empty command, free and continue
+    if (*tmpline == 0) {
+      free (tmpline);
+      if (q == line + Nline) {
+        done = TRUE;
+      } else {
+        p = q + 1;
+        myAssert (p - line <= Nline, "oops");
+      }
+      continue;
+    }
+      status = command (tmpline, &outline, verbose);
+    if (bufferPending) {
+      // flush old messages
+      ExpectMessage (server, 0.2, &message);
+      bufferPending = FALSE;
+    }
+      if (status == -1) {
+        if (server) {
+          // send the command to the server instead
+          if (!SendMessage (server, "%s", outline)) {
+            switch (errno) {
+              case EPIPE:
+                gprint (GP_ERR, "server connection has died\n");
+                exit (32);
+              default:
+                gprint (GP_ERR, "server is busy...32\n");
+                bufferPending = TRUE;
+                goto escape;
+            }
+    // input tmpline is freed by command
+    status = command (tmpline, &outline, verbose);
+    if (status == -1) {
+      if (server) {
+        // send the command to the server instead
+        if (!SendMessage (server, "%s", outline)) {
+          switch (errno) {
+            case EPIPE:
+              gprint (GP_ERR, "server connection has died\n");
+              exit (32);
+            default:
+              gprint (GP_ERR, "server is busy...32\n");
+              bufferPending = TRUE;
+              goto escape;
+          }
+        }
+          // receive the command exit status
+          if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+            switch (errno) {
+              case EPIPE:
+                gprint (GP_ERR, "server connection has died\n");
+                exit (33);
+              default:
+                gprint (GP_ERR, "server is busy...33\n");
+                bufferPending = TRUE;
+                goto escape;
+            }
+          } else {
+            sscanf (message.buffer, "STATUS %d", &status);
+          }
+          // receive the resulting stderr
+          if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+            switch (errno) {
+              case EPIPE:
+                gprint (GP_ERR, "server connection has died\n");
+                exit (34);
+              default:
+                gprint (GP_ERR, "server is busy...34\n");
+                bufferPending = TRUE;
+                goto escape;
+            }
+          } else {
+            gwrite (message.buffer, 1, message.Nbuffer, GP_ERR);
+          }
+          // receive the resulting stdout
+          if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+            switch (errno) {
+              case EPIPE:
+                gprint (GP_ERR, "server connection has died\n");
+                exit (35);
+              default:
+                gprint (GP_ERR, "server is busy...35\n");
+                bufferPending = TRUE;
+                goto escape;
+            }
+          } else {
+            gwrite (message.buffer, 1, message.Nbuffer, GP_LOG);
+        // receive the command exit status
+        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+          switch (errno) {
+            case EPIPE:
+              gprint (GP_ERR, "server connection has died\n");
+              exit (33);
+            default:
+              gprint (GP_ERR, "server is busy...33\n");
+              bufferPending = TRUE;
+              goto escape;
+          }
         } else {
+          // if no server is defined, we treat an unknown command as an error
+          status = FALSE;
+        }
+      }
+    escape:
+      if (outline != NULL) free (outline);
+      if (!status && auto_break) done = TRUE;
+          sscanf (message.buffer, "STATUS %d", &status);
+        }
+        // receive the resulting stderr
+        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+          switch (errno) {
+            case EPIPE:
+              gprint (GP_ERR, "server connection has died\n");
+              exit (34);
+            default:
+              gprint (GP_ERR, "server is busy...34\n");
+              bufferPending = TRUE;
+              goto escape;
+          }
+        } else {
+          gwrite (message.buffer, 1, message.Nbuffer, GP_ERR);
+        }
+        // receive the resulting stdout
+        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
+          switch (errno) {
+            case EPIPE:
+              gprint (GP_ERR, "server connection has died\n");
+              exit (35);
+            default:
+              gprint (GP_ERR, "server is busy...35\n");
+              bufferPending = TRUE;
+              goto escape;
+          }
+        } else {
+          gwrite (message.buffer, 1, message.Nbuffer, GP_LOG);
+        }
+      } else {
+        // if no server is defined, we treat an unknown command as an error
+        status = FALSE;
+      }
+    }
+  escape:
+    if (outline != NULL) free (outline);
+    if (!status && auto_break) done = TRUE;
     p = q + 1;
+  }

trunk/Ohana/src/opihi/lib.shell/parse.c

-              r41340
+              r41341
         /* dvomath returns a new string, or NULL, with the result of the expression */
         val = dvomath (1, &V1, &size, 0);
         if (val == NULL) goto error;
+        if (val == NULL) { print_error (); goto error; }
         fval += atof(val);
         // save the result
 …
         /* dvomath returns a new string, or NULL, with the result of the expression */
         val = dvomath (1, &V1, &size, 0);
         if (val == NULL) goto error;
+        if (val == NULL) { print_error (); goto error; }
         fval -= atof(val);
         // save the result

trunk/Ohana/src/opihi/lib.shell/stack_math.c

-              r41340
+              r41341
     case '?': SSS_FUNC(M1 ? M2: M3);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as scalar trinary op!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case '?': VVV_FUNC(*M1 ? *M2: *M3);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as vector trinary op!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case '?': MMM_FUNC(*M1 ? *M2: *M3);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as matrix trinary op!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'A': VV_FUNC(ST_SCALAR_INT, (*M1 && *M2) ? 1 : 0);
     case 'O': VV_FUNC(ST_SCALAR_INT, (*M1 || *M2) ? 1 : 0);
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->vector->type == OPIHI_FLT)) {
+        // bitshift is not valid with float valuess
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      // I could just do this for all types and bitshift regardless...
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
+        *out = *M1 << *M2;
+      }
+      break;
+    }
+    case 'r': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->vector->type == OPIHI_FLT)) {
+        // bitshift is not valid with float valuess
+        CopyVector (OUT[0].vector, V1[0].vector);
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      // I could just do this for all types and bitshift regardless...
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
+        *out = *M1 >> *M2;
+      }
+      break;
+    }
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'A': SV_FUNC(ST_SCALAR_INT, (M1 && *M2) ? 1 : 0);
     case 'O': SV_FUNC(ST_SCALAR_INT, (M1 || *M2) ? 1 : 0);
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      CopyVector (OUT[0].vector, V2[0].vector);
+      if ((V1->type == ST_SCALAR_FLT) || (V2->vector->type == OPIHI_FLT)) {
+        // bitshift is not valid with float valuess
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      opihi_int  M1  =  V1[0].IntValue;                                 \
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M2++) {
+        *out = M1 << *M2;
+      }
+      break;
+    }
+    case 'r': {
+      CopyVector (OUT[0].vector, V2[0].vector);
+      if ((V1->type == ST_SCALAR_FLT) || (V2->vector->type == OPIHI_FLT)) {
+        // bitshift is not valid with float valuess
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      opihi_int  M1  =  V1[0].IntValue;                                 \
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M2++) {
+        *out = M1 >> *M2;
+      }
+      break;
+    }
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'A': VS_FUNC(ST_SCALAR_INT, (*M1 && M2) ? 1 : 0);
     case 'O': VS_FUNC(ST_SCALAR_INT, (*M1 || M2) ? 1 : 0);
+    case 'l': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->type == ST_SCALAR_FLT)) {
+        // bitshift is not valid with float valuess
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int  M2  =  V2[0].IntValue;                                 \
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++) {
+        *out = *M1 << M2;
+      }
+      break;
+    }
+    case 'r': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->type == ST_SCALAR_FLT)) {
+        // bitshift is not valid with float valuess
+        CopyVector (OUT[0].vector, V1[0].vector);
+        for (i = 0; i < Nx; i++) {
+          OUT[0].vector[0].elements.Flt[i] = NAN;
+        }
+        break;
+      }
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int  M2  =  V2[0].IntValue;                                 \
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++) {
+        *out = *M1 >> M2;
+      }
+      break;
+    }
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'A': MV_FUNC((*M1 && *M2) ? 1 : 0);
     case 'O': MV_FUNC((*M1 || *M2) ? 1 : 0);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'O': VM_FUNC((*M1 || *M2) ? 1 : 0);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'O': MM_FUNC((*M1 || *M2) ? 1 : 0);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'O': MS_FUNC((*M1 || M2) ? 1 : 0);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'O': SM_FUNC((M1 || *M2) ? 1 : 0);
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
 …
     case 'A': SS_FUNC(ST_SCALAR_INT, (M1 && M2) ? 1 : 0);
     case 'O': SS_FUNC(ST_SCALAR_INT, (M1 || M2) ? 1 : 0);
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      if ((V1->type == ST_SCALAR_FLT) || (V2->type == ST_SCALAR_FLT)) {
+        // bitshift is not valid with float valuess
+        OUT[0].type = ST_SCALAR_FLT;
+        OUT[0].FltValue = NAN;
+        break;
+      }
+      opihi_int M1 = V1[0].IntValue;
+      opihi_int M2 = V2[0].IntValue;
+      OUT[0].type = ST_SCALAR_INT;
+      OUT[0].IntValue = M1 << M2;
+      break;
+    }
+    case 'r': {
+      if ((V1->type == ST_SCALAR_FLT) || (V2->type == ST_SCALAR_FLT)) {
+        // bitshift is not valid with float valuess
+        OUT[0].type = ST_SCALAR_FLT;
+        OUT[0].FltValue = NAN;
+        break;
+      }
+      opihi_int M1 = V1[0].IntValue;
+      opihi_int M2 = V2[0].IntValue;
+      OUT[0].type = ST_SCALAR_INT;
+      OUT[0].IntValue = M1 >> M2;
+      break;
+    }
     default:
       snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
 …
 int S_unary (StackVar *OUT, StackVar *V1, char *op) {
-  char line[512]; // this is only used to report an error
 # define S_FUNC(OP,FTYPE) {                                             \
     if (V1->type == ST_SCALAR_FLT) {                                    \
 …
   if (!strcmp (op, "asinh"))  S_FUNC(asinh (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "acosh"))  S_FUNC(acosh (M1), ST_SCALAR_FLT);
-  if (!strcmp (op, "lgamma")) S_FUNC(lgamma (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "sin"))    S_FUNC(sin (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "cos"))    S_FUNC(cos (M1), ST_SCALAR_FLT);
 …
   if (!strcmp (op, "dacos"))  S_FUNC(acos (M1)*DEG_RAD, ST_SCALAR_FLT);
   if (!strcmp (op, "datan"))  S_FUNC(atan (M1)*DEG_RAD, ST_SCALAR_FLT);
+  if (!strcmp (op, "lgamma")) S_FUNC(lgamma (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "rnd"))    S_FUNC(0.0*M1 + drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "drnd"))   S_FUNC(0.0*M1 + drand48(), ST_SCALAR_FLT);
 …
   clear_stack (V1);
+  snprintf (line, 512, "error: op %s not defined!", op);
+  char line[512]; // this is only used to report an error
+  snprintf (line, 512, "error: op %s not defined as unary scalar op!", op);
   push_error (line);
   return (FALSE);
+}
 …
   if (!strcmp (op, "asinh"))  V_FUNC(asinh(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "acosh"))  V_FUNC(acosh(*M1), ST_SCALAR_FLT);
-  if (!strcmp (op, "lgamma")) V_FUNC(lgamma(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "sin"))    V_FUNC(sin(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "cos"))    V_FUNC(cos(*M1), ST_SCALAR_FLT);
 …
   if (!strcmp (op, "dacos"))  V_FUNC(acos(*M1)*DEG_RAD, ST_SCALAR_FLT);
   if (!strcmp (op, "datan"))  V_FUNC(atan(*M1)*DEG_RAD, ST_SCALAR_FLT);
+  if (!strcmp (op, "lgamma")) V_FUNC(lgamma(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "rnd"))    V_FUNC(drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "drnd"))   V_FUNC(drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "lrnd"))   V_FUNC(lrand48(), ST_SCALAR_INT);
   if (!strcmp (op, "mrnd"))   V_FUNC(mrand48(), ST_SCALAR_INT);
-  if (!strcmp (op, "ramp"))   V_FUNC(i, ST_SCALAR_INT);
-  if (!strcmp (op, "zero"))   V_FUNC(0, ST_SCALAR_INT);
   if (!strcmp (op, "not"))    V_FUNC(!(*M1), ST_SCALAR_INT);
   if (!strcmp (op, "--"))     V_FUNC(-1*(*M1), ST_SCALAR_INT); // NOTE: opihi_int is signed
   if (!strcmp (op, "isinf"))  V_FUNC(!finite(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "isnan"))  V_FUNC(isnan((opihi_flt)(*M1)), ST_SCALAR_FLT);
+  if (!strcmp (op, "ramp"))   V_FUNC(i, ST_SCALAR_INT);
   if (!strcmp (op, "xramp"))  V_FUNC(i, ST_SCALAR_INT);
   if (!strcmp (op, "yramp"))  V_FUNC(0, ST_SCALAR_INT);
   if (!strcmp (op, "zramp"))  V_FUNC(0, ST_SCALAR_INT);
+  if (!strcmp (op, "zero"))   V_FUNC(0, ST_SCALAR_INT);
   /* xramp, yramp, zramp above only make sense for matrices. for vectors, xramp = ramp, yramp = zero */
 # undef V_FUNC
+escape:
+  // free the temp vector if needed
   if (V1[0].type == ST_VECTOR_TMP) {
     free (V1[0].vector[0].elements.Ptr);
 …
   clear_stack (V1);
+  char line[512]; // this is only used to report an error
+  snprintf (line, 512, "error: op %s not defined as unary vector op!", op);
+  push_error (line);
+  return (FALSE);
+escape:
+  if (V1[0].type == ST_VECTOR_TMP) {
+    free (V1[0].vector[0].elements.Ptr);
+    free (V1[0].vector);
+    V1[0].vector = NULL;
+  }
+  clear_stack (V1);
   return (TRUE);
+}
 # define M_FUNC(OP) { for (i = 0; i < Npix; i++, out++, M1++) { *out = (OP); }}
+# define M_FUNC(OP) { for (i = 0; i < Npix; i++, out++, M1++) { *out = (OP); } goto escape; }
 int M_unary (StackVar *OUT, StackVar *V1, char *op) {
 …
 // if (!strcmp (op, "rint"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = nearbyint (*M1); }}
   if (!strcmp (op, "="))     { }
+  if (!strcmp (op, "="))      { goto escape; }
   if (!strcmp (op, "abs"))    M_FUNC(fabs(*M1));
   if (!strcmp (op, "int"))    M_FUNC((opihi_flt)(long long)(*M1));
 …
   if (!strcmp (op, "asinh"))  M_FUNC(asinh(*M1));
   if (!strcmp (op, "acosh"))  M_FUNC(acosh(*M1));
-  if (!strcmp (op, "lgamma")) M_FUNC(lgamma(*M1));
   if (!strcmp (op, "sin"))    M_FUNC(sin(*M1));
   if (!strcmp (op, "cos"))    M_FUNC(cos(*M1));
 …
   if (!strcmp (op, "dacos"))  M_FUNC(acos(*M1)*DEG_RAD);
   if (!strcmp (op, "datan"))  M_FUNC(atan(*M1)*DEG_RAD);
+  if (!strcmp (op, "not"))    M_FUNC(!(*M1));
+  if (!strcmp (op, "--"))     M_FUNC(-(*M1));
+  if (!strcmp (op, "lgamma")) M_FUNC(lgamma(*M1));
   if (!strcmp (op, "rnd"))    M_FUNC(drand48());
   if (!strcmp (op, "drnd"))   M_FUNC(drand48());
   if (!strcmp (op, "lrnd"))   M_FUNC(lrand48());
   if (!strcmp (op, "mrnd"))   M_FUNC(mrand48());
+  if (!strcmp (op, "not"))    M_FUNC(!(*M1));
+  if (!strcmp (op, "--"))     M_FUNC(-(*M1));
   if (!strcmp (op, "ramp"))   M_FUNC(i);
-  if (!strcmp (op, "zero"))   M_FUNC(0);
   if (!strcmp (op, "isinf"))  M_FUNC(!finite(*M1));
   if (!strcmp (op, "isnan"))  M_FUNC(isnan(*M1));
+  if (!strcmp (op, "zero"))   M_FUNC(0);
   /* xrm and yrm only make sense for 2D matrices. see special meaning for vectors */
 …
+      }
+    }
+    goto escape;
+  }
   if (!strcmp (op, "yramp")) {
 …
+      }
+    }
+    goto escape;
+  }
   if (!strcmp (op, "zramp")) {
 …
+      }
+    }
+  }
+    goto escape;
+  }
+  if (V1[0].type == ST_MATRIX_TMP) {
+    free (V1[0].buffer[0].header.buffer);
+    free (V1[0].buffer[0].matrix.buffer);
+    free (V1[0].buffer);
+  }
+  clear_stack (V1);
+  char line[512]; // this is only used to report an error
+  snprintf (line, 512, "error: op %s not defined as unary matrix op!", op);
+  push_error (line);
+  return (FALSE);
+ escape:
   if (V1[0].type == ST_MATRIX_TMP) {
     free (V1[0].buffer[0].header.buffer);

trunk/Ohana/src/opihi/lib.shell/startup.c

-              r41340
+              r41341
       LOAD_RC = FALSE;
+    }
+    if ((N = get_argument (*argc, argv, "--no-rc"))) {
+      remove_argument (N, argc, argv);
+      LOAD_RC = FALSE;
+    }
+    if ((N = get_argument (*argc, argv, "-norc"))) {
+      remove_argument (N, argc, argv);
+      LOAD_RC = FALSE;
+    }
+    if ((N = get_argument (*argc, argv, "-no-rc"))) {
+      remove_argument (N, argc, argv);
+      LOAD_RC = FALSE;
+    }
     ONLY_INPUT = FALSE;

trunk/Ohana/src/opihi/mana/Makefile

-              r41340
+              r41341
 $(SRC)/deimos_getobj.$(ARCH).o \
 $(SRC)/deimos_fitobj.$(ARCH).o \
-$(SRC)/deimos_fitalt.$(ARCH).o \
-$(SRC)/deimos_mkalt.$(ARCH).o \
 $(SRC)/deimos_mkmodel.$(ARCH).o \
 $(SRC)/deimos_arclines.$(ARCH).o \
+$(SRC)/deimos_fitarc.$(ARCH).o \
+$(SRC)/deimos_fitprofile.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o

trunk/Ohana/src/opihi/mana/deimos.c

-              r41340
+              r41341
 int deimos_fitslit (int argc, char **argv);
 int deimos_mkobj (int argc, char **argv);
-int deimos_mkalt (int argc, char **argv);
 int deimos_getobj (int argc, char **argv);
+// int deimos_getalt (int argc, char **argv);
 int deimos_fitobj (int argc, char **argv);
 int deimos_fitalt (int argc, char **argv);
+// int deimos_fitalt (int argc, char **argv);
 int deimos_arclines (int argc, char **argv);
+int deimos_fitarc (int argc, char **argv);
+int deimos_fitprofile (int argc, char **argv);
 static Command deimos_commands[] = {
+  {1, "fitobj", deimos_fitobj, "fit for object parameters"},
+  {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"},
+  {1, "getobj", deimos_getobj, "determine crude object parameters"},
+  {1, "mkobj", deimos_mkobj, "make a full object image"},
+  {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"},
+  {1, "mkslit", deimos_mkslit, "make a slit image"},
+  {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
+  {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"},
+  {1, "fitobj",     deimos_fitobj,     "fit for object parameters using LMM"},
+  {1, "fitalt",     deimos_fitobj,     "fit for object parameters using LMM"},
+  {1, "getobj",     deimos_getobj,     "determine crude object parameters"},
+  {1, "getalt",     deimos_getobj,     "determine crude object parameters"},
+  {1, "mkobj",      deimos_mkobj,      "make a full object image"},
+  {1, "mkslit",     deimos_mkslit,     "make a slit image"},
+  {1, "fitslit",    deimos_fitslit,    "fit slit image to observed slit flux"},
+  {1, "arclines",   deimos_arclines,   "detect arclines using LSF and STILT"},
+  {1, "fitarc",     deimos_fitarc,     "fit arclamp lines"},
+  {1, "fitprofile", deimos_fitprofile, "fit slit profile"},
 };

trunk/Ohana/src/opihi/mana/deimos_arclines.c

-              r41340
+              r41341
   deimos_make_LSF_kernel (LSF, stilt, Nx);
+  // FILE *Fout = NULL;
   // loop over the rows
   for (int iy = 0; iy < Ny; iy++) {
     float Fsum = 0;
+    float Ksum = 0;
+    // if (iy == 5802) Fout = fopen ("test.5802.dat", "w");
+    // if (iy == 5815) Fout = fopen ("test.5815.dat", "w");
     // cross-correlation of buffer values and kernel values
     for (int ix = 0; ix < Nx; ix++) {
 …
         if ((iy + ky) < 0) continue;
         if ((iy + ky) >= Ny) continue;
+        if (!isfinite(Fin[ix + (iy + ky)*Nx])) continue;
+        if (!isfinite(kernel[ix + ko*Nx])) continue;
         Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
+        Ksum += kernel[ix + ko*Nx];
+        // if ((iy == 5815) || (iy == 5802)) {
+        //   fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]);
+        // }
+      }
+    }
+    // if ((iy == 5815) || (iy == 5802)) fclose (Fout);
     coord->elements.Flt[iy] = iy;
     flux->elements.Flt[iy] = Fsum;

trunk/Ohana/src/opihi/mana/deimos_fitslit.c

-              r41340
+              r41341
 # include "data.h"
+# define IRLS_TOLERANCE 1e-4
+float weight_cauchy_square_flt (float x2);
+static void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky);
 int deimos_fitslit (int argc, char **argv) {
+  // fitslit (window) (slit) (flux) (sky)
+  // int N;
+  // fitslit (window) (vari) (slit) (flux) (sky)
   Vector *flux    = NULL;
 …
   Buffer *wind    = NULL;
+  Buffer *vari    = NULL;
   Buffer *slit    = NULL;
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: deimos fitslit (window) (slit) (flux) (sky)\n");
+    gprint (GP_ERR, "  inputs:  window (observed 2D flux), slit (model 2D flux)\n");
+  int N;
+  int Niter = 10;
+  if ((N = get_argument (argc, argv, "-irls-iter"))) {
+    remove_argument (N, &argc, argv);
+    Niter = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  Buffer *model = NULL;
+  if ((N = get_argument (argc, argv, "-model"))) {
+    remove_argument (N, &argc, argv);
+    if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: deimos fitslit (window) (variance) (slit) (flux) (sky)\n");
+    gprint (GP_ERR, "  inputs:  window (observed 2D flux), variance (on 2D flux), slit (model 2D flux)\n");
     gprint (GP_ERR, "  outputs: flux (best-fit 1D flux), sky (best-fit 1D background)\n");
     return FALSE;
 …
   // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
   if ((wind = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((slit = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((flux = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sky  = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  // XXX confirm slit and wind have same dimensions
+  if ((vari = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((slit = SelectBuffer (argv[3], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((flux = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sky  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
   // define the output window
 …
   int Ny = wind[0].matrix.Naxis[1];
+  // confirm wind, vari, slit have same dimensions
   if (Nx != slit[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
   if (Ny != slit[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
+  if (Nx != vari[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
+  if (Ny != vari[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
   ResetVector (flux, OPIHI_FLT, Ny);
   ResetVector (sky,  OPIHI_FLT, Ny);
+  if (model) {
+    if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; }
+  }
   float *Fwind = (float *) wind[0].matrix.buffer;
+  float *Fvari = (float *) vari[0].matrix.buffer;
   float *Fprof = (float *) slit[0].matrix.buffer;
+  float *FmodOut = model ? (float *) model[0].matrix.buffer : NULL;
+  ALLOCATE_PTR (weight, float, Nx);
+  ALLOCATE_PTR (rawwgt, float, Nx);
   // loop over the rows
   for (int iy = 0; iy < Ny; iy++) {
     // calculate elements of the chi-square for this row:
+    double F = 0.0, S = 0.0;
+    double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
+    double wt = 1.0;
+    // set weight based on the variance
     for (int ix = 0; ix < Nx; ix++) {
+      float fxy = Fwind[ix + iy*Nx];
+      float pxy = Fprof[ix + iy*Nx];
+      R  += wt;
+      P  += pxy*wt;
+      P2 += pxy*pxy*wt;
+      f  += fxy*wt;
+      fp += fxy*pxy*wt;
+      float vari = Fvari[ix + iy*Nx];
+      weight[ix] = !isfinite(vari) || (fabs(vari) < 1e-6) ? 1.0 : 1.0 / vari;
+      rawwgt[ix] = weight[ix];
+    }
+    double det = R*P2 - P*P;
+    double S = (P2*f - P*fp) / det;
+    double F = (R*fp - P*f) / det;
+    // calculate elements of the chi-square for this row:
+    fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
+    // do an IRLS loop here to reject outliers
+    int converged = FALSE;
+    for (int iter = 0; (iter < Niter) && !converged; iter++) {
+      double Flast = F;
+      double Slast = S;
+      // calculate weight modification based on distances (squared).
+      // use modifier to calculate new weighted mean
+      for (int ix = 0; ix < Nx; ix++) {
+        // F & S are model parameters; compare the model and observed values for each pixel
+        float Fmodel = S + F*Fprof[ix + iy*Nx];
+        float dV = (Fwind[ix + iy*Nx] - Fmodel);
+        float d2 = SQ(dV) * rawwgt[ix];
+        float Mod = weight_cauchy_square_flt (d2);
+        weight[ix] = Mod * rawwgt[ix];
+      }
+      fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
+      float dF = fabs(F - Flast);
+      float dS = fabs(S - Slast);
+      if ((dF < F * IRLS_TOLERANCE) && (dS < S * IRLS_TOLERANCE)) converged = TRUE;
+    }
+    if (FmodOut) {
+      for (int ix = 0; ix < Nx; ix++) {
+        // F & S are model parameters; compare the model and observed values for each pixel
+        float Fmodel = S + F*Fprof[ix + iy*Nx];
+        float dV = (Fwind[ix + iy*Nx] - Fmodel);
+        float d2 = SQ(dV) * rawwgt[ix];
+        float Mod = weight_cauchy_square_flt (d2);
+        if (Mod < 0.1) {
+          FmodOut[ix + iy*Nx] = NAN;
+        } else {
+          FmodOut[ix + iy*Nx] = Fmodel;
+        }
+      }
+    }
     flux->elements.Flt[iy] = F;
     sky->elements.Flt[iy] = S;
+  }
   return TRUE;
+}
+// weight is 1 / variance
+// Fwind, Fprof are pointers to the start of a single row
+// Fwind is the observed flux in the window
+// Fprof is the slit profile
+void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky) {
+  // calculate elements of the chi-square fit for this row:
+  double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
+  for (int ix = 0; ix < Nx; ix++) {
+    double fxy = Fwind[ix];
+    double pxy = Fprof[ix];
+    // skip NAN / Inf pixels in either window or slit
+    if (!isfinite(fxy)) continue;
+    if (!isfinite(pxy)) continue;
+    double wt = (weight == NULL) ? 1.0 : weight[ix];
+    R  += wt;
+    P  += pxy*wt;
+    P2 += pxy*pxy*wt;
+    f  += fxy*wt;
+    fp += fxy*pxy*wt;
+  }
+  double det = R*P2 - P*P;
+  *sky = (P2*f - P*fp) / det;
+  *flux = (R*fp - P*f) / det;
+}

trunk/Ohana/src/opihi/mana/deimos_mkalt.c

-              r41340
+              r41341
 # include "data.h"
 # include "deimos.h"
+XXX: deprecated
 // Use this to test the deimos_make_object function used by deimos_fitobj
 …
   // generate the model based on the test values
   float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, Nx, Ny, 0);
+  float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, NULL, Nx, Ny, 0);
   ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0);

trunk/Ohana/src/opihi/mana/deimos_mkmodel.c

-              r41340
+              r41341
 // generate a subimage of the full model for the range row - row + Ny
 // obj,sky,bck are subset vectors corresponding to the range row to row+Ny
+float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) {
+  float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row);
+// the slit_trace_* and psf_trace refer to the full image starting at row
+float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row) {
+  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
   float *new = deimos_apply_tilt (raw, Nx, Ny);
   deimos_apply_profile (profile, new, Nx, Ny);
   deimos_add_background (bck, new, Nx, Ny, row);
   float *out = deimos_apply_trace (trace, new, Nx, Ny, row);
+  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
   free (raw);
 …
 // *** generate a raw image with no slit tilt and no trace offsets.  apply the PSF to the
 // object flux, add in the sky flux
 float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) {
+float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) {
   OHANA_UNUSED_PARAM(row);
 …
   // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
   int Xoff = (int)(Xref - Npsf);
+  int XoffRef = (int)(Xref - Npsf);
   // Nx : width of output window
 …
     // we are generating the image for just the row range row to row+Ny
+    // note that obj & sky are subset vectors for just this range of pixels
     opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
     opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
+    // integral and fractional pixel offsets of PSF (none, by default)
+    int dxi = 0;
+    float dxf = 0.0, dxr = 1.0;
+    int Xoff = XoffRef;
+    if (psf_trace) {
+      // the psf_trace is referenced against the full image, so we need to add row to iy:
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
+      dxi = floor(dx);
+      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+      dxr = 1 - dxf;
+      Xoff += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
+    }
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     // flux = obj * PSF + sky
     for (int ix = 0; ix < Nx; ix++) {
 …
       // only add in the flux if we are in range of the PSF
       if ((n >= 0) && (n < NpsfFull)) {
+        value += objVy * psfV[n];
+        if (doInterp) {
+          if ((n > 0) && (n < NpsfFull)) {
+            if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+            }
+          }
+          if (n == 0) {
+            if (isfinite(psfV[n])) {
+              value += objVy*psfV[n]*dxr;
+            }
+          }
+        } else {
+          if (isfinite(psfV[n])) {
+            value += objVy * psfV[n];
+          }
+        }
+      }
       out[ix + iy*Nx] = value;
 …
 // shift pixels in x based on the trace
 float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) {
+float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
   ALLOCATE_PTR (output, float, Nx*Ny);
 …
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     // NOTE: spline is evaluated at full dispersion coordinate
+    float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row);
+    int iy_full = iy + row;
+    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
+    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
     // extract the integer pixel offset and the fractional offset

trunk/Ohana/src/opihi/mana/deimos_mkobj.c

-              r41340
+              r41341
   // input parameters:
+  // * trace       : spline fit of slit central x pos vs y-coord
+  // * profile     : slit window profile (vector)
+  // * object      : vector of object flux vs y-coord
+  // * sky         : vector of local sky signal vs y-coord
+  // * background  : vector of extra-slit background flux vs y-coord
+  // * PSF         : point-spread function vector (flux normalized, x-dir)
+  // * LSF         : sline-spread function vector (flux normalized, y-dir)
+  // * stilt       : slit tilt response : 2D kernel?
+  // * slit_trace_red : spline fit of slit central x pos vs y-coord for red chip
+  // * slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip
+  // * psf_trace      : spline fit of psf central x pos vs y-coord
+  // * profile        : slit window profile (vector)
+  // * object         : vector of object flux vs y-coord
+  // * sky            : vector of local sky signal vs y-coord
+  // * backgnd        : vector of extra-slit background flux vs y-coord
+  // * PSF            : point-spread function vector (flux normalized, x-dir)
+  // * LSF            : sline-spread function vector (flux normalized, y-dir)
+  // * stilt          : slit tilt response : 2D kernel?
+  // * redlimit       : pixel break between red and blu chips
+  // * Nwave          : number of wavelength pixels in output image
   int N;
   // if any of these are not defined, they will have assumed identity values
+  Spline *trace   = NULL;
+  Vector *profile = NULL;
+  Vector *object  = NULL;
+  Vector *sky     = NULL;
+  Vector *backgnd = NULL;
+  Vector *psf     = NULL;
+  Vector *lsf     = NULL;
+  Spline *slit_trace_red = NULL;
+  Spline *slit_trace_blu = NULL;
+  Spline *psf_trace  = NULL;
+  Vector *profile    = NULL;
+  Vector *object     = NULL;
+  Vector *sky        = NULL;
+  Vector *backgnd    = NULL;
+  Vector *psf        = NULL;
+  Vector *lsf        = NULL;
   // add tilt later
   Buffer *output  = NULL;
+  Buffer *output     = NULL;
   float stilt = 0.0; // angle of the slit
+  if ((N = get_argument (argc, argv, "-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
+  // for a red vs blu spline, we need to specify the split point
+  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
+  // this more generic (need to define the ranges somewhere)
+  int redlimit = 4096;
+  if ((N = get_argument (argc, argv, "-redlimit"))) {
+    remove_argument (N, &argc, argv);
+    redlimit = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-psf-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
+  }
 …
+  }
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
   if (argc != 3) {
+    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
+    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-object vector] [-sky vector] [-backgnd vector] [-slit-trace-red spline] [-slit-trace-blu spline] [-psf-trace spline] [-profile vector] [-psf vector] [-lsf vector] [-stilt angle] [-Nwave Npixel] [-redlimit pixel]\n");
+    gprint (GP_ERR, "  slit_trace_red : spline fit of slit central x pos vs y-coord for red chip\n");
+    gprint (GP_ERR, "  slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip\n");
+    gprint (GP_ERR, "  psf_trace      : spline fit of psf central x pos vs y-coord\n");
+    gprint (GP_ERR, "  profile        : slit window profile (vector)\n");
+    gprint (GP_ERR, "  object         : vector of object flux vs y-coord \n");
+    gprint (GP_ERR, "  sky            : vector of local sky signal vs y-coord\n");
+    gprint (GP_ERR, "  backgnd        : vector of extra-slit background flux vs y-coord\n");
+    gprint (GP_ERR, "  PSF            : point-spread function vector (flux normalized, x-dir)\n");
+    gprint (GP_ERR, "  LSF            : sline-spread function vector (flux normalized, y-dir)\n");
+    gprint (GP_ERR, "  stilt          : slit tilt response : 2D kernel? \n");
+    gprint (GP_ERR, "  redlimit       : pixel break between red and blu chips\n");
+    gprint (GP_ERR, "  Nwave          : number of wavelength pixels in output image\n");
     return FALSE;
+  }
 …
   // if we are appying a trace offset spline, we need to generate an output window which
   // is Nx + the full swing of the trace, then window back down
   if (trace) {
+  if (slit_trace_red && slit_trace_blu) {
     float dXmin = +1000;
     float dXmax = -1000;
     for (int iy = 0; iy < Ny; iy++) {
       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
+      float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
+      float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       dXmin = MIN (dx, dXmin);
       dXmax = MAX (dx, dXmax);
 …
     if (psf) Npof -= Npsf;
+    // integral and fractional pixel offsets of PSF (none, by default)
+    int dxi = 0;
+    float dxf = 0.0, dxr = 1.0;
+    if (psf_trace) {
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
+      dxi = floor(dx);
+      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+      dxr = 1 - dxf;
+      Npof += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
+      if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
+    }
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     // flux = obj * PSF + sky
     for (int ix = 0; ix < Nx; ix++) {
 …
         // only add in the flux if we are in range of the PSF
         if ((n >= 0) && (n < NpsfFull)) {
+          value += objVy * psfV[n];
+          if (doInterp) {
+            if ((n > 0) && (n < NpsfFull)) {
+              if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+                value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+              }
+            }
+            if (n == 0) {
+              if (isfinite(psfV[n])) {
+                value += objVy*psfV[n]*dxr;
+              }
+            }
+          } else {
+            if (isfinite(psfV[n])) {
+              value += objVy*psfV[n];
+            }
+          }
+        }
       } else {
 …
   // shift pixels in x based on the trace
   if (trace) {
+  if (slit_trace_red && slit_trace_blu) {
     ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
     float *outTrace = (float *) outTraceBuffer;
 …
       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
+      float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
+      float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       // extract the integer pixel offset and the fractional offset

trunk/Ohana/src/opihi/mana/deimos_mkslit.c

-              r41340
+              r41341
   // input parameters:
   //   slit profile response : vector of fractional flux vs x-coord
+  //   trace                 : spline fit of slit central x pos vs y-coord
+  //   flux                  : vector of input signal vs y-coord
+  //   trace_ref,trace_blu   : spline fit of slit central x pos vs y-coord
+  //   flux (optional)       : vector of input signal vs y-coord
+  //   sky  (optional)       : vector of input signal vs y-coord
+  //   fluxbins (optional)   : size of output image in y-direction
+  //   NOTE: one of flux, sky, or fluxbins must be defined
+  // output : buff (an image with Nx defined by profile and Ny defined
+  //                by flux, sky, or fluxbins)
   int N;
+  Vector *profile = NULL;
+  Vector *flux    = NULL;
+  Vector *sky     = NULL;
+  Spline *trace   = NULL;
+  Buffer *buff    = NULL;
+  Vector *profile   = NULL;
+  Vector *flux      = NULL;
+  Vector *sky       = NULL;
+  Spline *trace_red = NULL;
+  Spline *trace_blu = NULL;
+  Buffer *buff      = NULL;
   // user-specified flux value (otherwise assumed to be 1)
 …
+  }
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
+  // either flux or fluxbins must be specified to define output size
+  int profilebins = 0;
+  if ((N = get_argument (argc, argv, "-profilebins"))) {
+    remove_argument (N, &argc, argv);
+    profilebins = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  // for a red vs blu spline, we need to specify the split point
+  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
+  // this more generic (need to define the ranges somewhere)
+  int redlimit = 4096;
+  if ((N = get_argument (argc, argv, "-redlimit"))) {
+    remove_argument (N, &argc, argv);
+    redlimit = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace_red) (trace_blu) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
     return FALSE;
+  }
 …
   // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
+  if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
+  if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((profile   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((trace_red = FindSpline (argv[2])) == NULL) return (FALSE);
+  if ((trace_blu = FindSpline (argv[3])) == NULL) return (FALSE);
+  if ((buff      = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
   // define the output window
+  int Nx = 2*profile[0].Nelements + 1;
+  // if no profile half-width is specified, use supplied profile
+  int Nx = profilebins ? 2*profilebins + 1 : 2*profile[0].Nelements + 1;
   int Ny = flux ? flux[0].Nelements : fluxbins;
 …
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
+    Spline *trace = (iy < redlimit) ? trace_red : trace_blu;
     float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
 …
       // set the sky in the entire window
       out[ix + iy*Nx] = sky ? sky->elements.Flt[iy] : 0.0;
+      out[ix + iy*Nx] = NAN;
       // equivalent coord in the profile:
 …
       // a default value:
+      float vsky = sky ? sky->elements.Flt[iy] : 0.0;
       float vout = NAN;
 …
       if (flux) vout *= flux->elements.Flt[iy];
       out[ix + iy*Nx] += vout;
+      out[ix + iy*Nx] = vsky + vout;
+    }
+  }

trunk/Ohana/src/opihi/pantasks/Makefile

-              r41340
+              r41341
 $(SRC)/task_command.$(ARCH).o \
 $(SRC)/task_options.$(ARCH).o \
+$(SRC)/version.$(ARCH).o \
+$(SRC)/ipptool2book.$(ARCH).o
+$(SRC)/version.$(ARCH).o
 libs = \

trunk/Ohana/src/opihi/pantasks/init.c

-              r41340
+              r41341
 int verbose         PROTO((int, char **));
 int version         PROTO((int, char **));
-int ipptool2book    PROTO((int, char **));
 static Command cmds[] = {
 …
   {1, "host",       task_host,     "define host machine for a task"},
   {1, "nice",       task_nice,     "set nice level for a task"},
-  {1, "ipptool2book", ipptool2book, "convert queue with ipptool output to book"},
   {1, "kill",       kill_job,      "kill job"},
   {1, "nmax",       task_nmax,     "define maximum number of jobs for a task"},

trunk/Ohana/src/opihi/pantasks/init_server.c

-              r41340
+              r41341
 int version         PROTO((int, char **));
 int server          PROTO((int, char **));
-int ipptool2book    PROTO((int, char **));
 int server_run      PROTO((int, char **));
 …
   {1, "host",       task_host,     "define host machine for a task"},
   {1, "nice",       task_nice,     "set nice level for a task"},
-  {1, "ipptool2book", ipptool2book, "convert queue with ipptool output to book"},
   {1, "kill",       kill_job,      "kill job"},
   {1, "nmax",       task_nmax,     "define maximum number of jobs for a task"},

Context Navigation

Legend:

Download in other formats: