Index: trunk/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 41341)
@@ -39,4 +39,11 @@
   }
 
+  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"))) {
@@ -80,6 +87,9 @@
   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);
Index: trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c	(revision 41341)
@@ -13,4 +13,11 @@
   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;
@@ -79,6 +86,9 @@
   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);
@@ -90,6 +100,6 @@
   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;
Index: trunk/Ohana/src/opihi/cmd.astro/cgrid.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cgrid.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/cgrid.c	(revision 41341)
@@ -83,4 +83,7 @@
   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"))) {
@@ -115,6 +118,6 @@
   }
 
-  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);
@@ -151,8 +154,5 @@
   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? */
@@ -182,4 +182,5 @@
     }
   }
+
   
   /* set spacings for DEC */
@@ -345,3 +346,24 @@
   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);
 }
+
Index: trunk/Ohana/src/opihi/cmd.astro/csystem.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 41341)
@@ -12,4 +12,14 @@
   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
@@ -84,7 +94,45 @@
 
     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);
@@ -154,4 +202,9 @@
   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;
 }
Index: trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c	(revision 41341)
@@ -11,4 +11,17 @@
     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);
   }
@@ -101,4 +114,13 @@
       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; }
   }
 
@@ -169,4 +191,20 @@
 	}
       }
+    }
+  }
+  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];
     }
   }
Index: trunk/Ohana/src/opihi/cmd.astro/mkgauss.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/mkgauss.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/mkgauss.c	(revision 41341)
@@ -10,4 +10,10 @@
   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
@@ -60,4 +66,6 @@
   /* 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++) {
@@ -67,5 +75,5 @@
       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;
     }
Index: trunk/Ohana/src/opihi/cmd.astro/region.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 41340)
+++ trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 41341)
@@ -3,5 +3,5 @@
 int region (int argc, char **argv) {
   
-  double Ra, Dec, Radius;
+  double Ra, Dec;
   float dx, dy;
   int N, kapa, NoClear, dXpix, dYpix;
@@ -54,4 +54,26 @@
   }
 
+  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);
@@ -81,6 +103,11 @@
   }
 
-  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", 
@@ -90,28 +117,51 @@
 	     &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;
@@ -126,5 +176,6 @@
   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);
@@ -132,15 +183,38 @@
 
   /* 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);
