Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 36489)
@@ -28,7 +28,9 @@
 $(SRC)/cval.$(ARCH).o		   \
 $(SRC)/czplot.$(ARCH).o	   \
+$(SRC)/cdensify.$(ARCH).o	   \
 $(SRC)/drizzle.$(ARCH).o	   \
 $(SRC)/flux.$(ARCH).o		   \
 $(SRC)/fitplx.$(ARCH).o	   \
+$(SRC)/fitpm.$(ARCH).o	   \
 $(SRC)/fixwrap.$(ARCH).o	   \
 $(SRC)/fixcols.$(ARCH).o	   \
Index: /trunk/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 36489)
+++ /trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 36489)
@@ -0,0 +1,192 @@
+# include "data.h"
+
+# define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
+enum {IS_DOT, IS_SQUARE, IS_CIRCLE, IS_GAUSS};
+
+int cdensify (int argc, char **argv) {
+
+  int i, Nx, Ny, Xb, Yb, N, Xpix, Ypix;
+  double Xmin, Xmax, dX, Ymin, Ymax, dY, ix, iy;
+  float *val;
+  Buffer *bf;
+  Vector *vr, *vd;
+  opihi_flt *r, *d, x, y;
+  int kapa;
+  Graphdata graphmode;
+
+  int Normalize = TRUE;
+  if ((N = get_argument (argc, argv, "-raw"))) {
+    remove_argument (N, &argc, argv);
+    Normalize = FALSE;
+  }
+
+  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;
+  double Rmin = graphmode.coords.crval1 - 182.0;
+  double Rmax = graphmode.coords.crval1 + 182.0;
+
+  float scale = 0.0;
+  if ((N = get_argument (argc, argv, "-scale"))) {
+    remove_argument (N, &argc, argv);
+    scale = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int PSFTYPE = IS_DOT;
+  if ((N = get_argument (argc, argv, "-psf"))) {
+    remove_argument (N, &argc, argv);
+    if (!strcasecmp(argv[N], "dot"))    PSFTYPE = IS_DOT;
+    if (!strcasecmp(argv[N], "square")) PSFTYPE = IS_SQUARE;
+    if (!strcasecmp(argv[N], "circle")) PSFTYPE = IS_CIRCLE;
+    if (!strcasecmp(argv[N], "gauss"))  PSFTYPE = IS_GAUSS;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: cdensify buffer R D\n");
+    gprint (GP_ERR, " option: -psf [dot] (circle) (square) (gauss)\n");
+    return (FALSE);
+  }
+  
+  if ((bf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((vr = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vd = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if (vr[0].Nelements != vd[0].Nelements) return (FALSE);
+
+  REQUIRE_VECTOR_FLT (vr, FALSE); 
+  REQUIRE_VECTOR_FLT (vd, FALSE); 
+
+  KapaGetImageRange (kapa, &Xmin, &Xmax, &Ymax, &Ymin, &Xpix, &Ypix);
+  Xmax = graphmode.xmax;
+  Xmin = graphmode.xmin;
+  Ymax = graphmode.ymax;
+  Ymin = graphmode.ymin;
+  dX = (Xmax - Xmin) / (Xpix - 1);
+  dY = (Ymax - Ymin) / (Ypix - 1);
+
+  CHECKVAL(Xmin);
+  CHECKVAL(Xmax);
+  CHECKVAL(dX);
+
+  CHECKVAL(Ymin);
+  CHECKVAL(Ymax);
+  CHECKVAL(dY);
+
+  Nx = (Xmax - Xmin) / dX + 1;
+  Ny = (Ymax - Ymin) / dY + 1;
+  
+  Coords newcoords = graphmode.coords;
+  newcoords.cdelt1 *= dX;
+  newcoords.cdelt2 *= dY;
+  newcoords.crpix1 = (newcoords.crpix1 - Xmin) / dX;
+  newcoords.crpix2 = (newcoords.crpix2 - Ymin) / dY;
+
+  gfits_free_matrix (&bf[0].matrix);
+  gfits_free_header (&bf[0].header);
+  CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
+  strcpy (bf[0].file, "(empty)");
+  PutCoords (&newcoords, &bf[0].header);
+  
+  float scalescale = scale*scale;
+  float scale2 = (scale + 1.0) * (scale + 1.0);
+  float fSquare = 1.0 / scale2;
+  float fCircle = 1.0 / (3.141592 * scale2);
+  float fSigma  = 0.5 / scale2;
+  float fGauss  = 1.0 / (2.0 * 3.141592 * scale2);
+
+  // generate the PSF in a local tangent plane
+  Coords coords;
+  coords.crpix1 = coords.crpix2 = 0.0;
+  coords.crval1 = coords.crval2 = 0.0;
+  coords.cdelt1 = coords.cdelt2 = 1.0;
+  coords.pc1_1  = coords.pc2_2  = 1.0;
+  coords.pc1_2  = coords.pc2_1  = 0.0;
+  coords.Npolyterms = 0;
+  strcpy (coords.ctype, "DEC--TAN");
+
+  r = vr[0].elements.Flt;
+  d = vd[0].elements.Flt;
+  val = (float *)bf[0].matrix.buffer;
+  for (i = 0; i < vr[0].Nelements; i++, r++, d++) {
+    double rn = ohana_normalize_angle (*r);
+    while (rn < Rmin) rn += 360.0;
+    while (rn > Rmax) rn -= 360.0;
+    coords.crval1 = rn;
+    coords.crval2 = *d;
+
+    switch (PSFTYPE) {
+      case IS_DOT:
+	RD_to_XY (&x, &y, rn, *d, &graphmode.coords);
+	Xb = (x - Xmin) / dX;
+	Yb = (y - Ymin) / dY;
+
+	RD_to_XY (&x, &y, rn, *d, &newcoords);
+	Xb = (int) x;
+	Yb = (int) y;
+	if (Xb >= Nx) continue;
+	if (Yb >= Ny) continue;
+	if (Xb < 0) continue;
+	if (Yb < 0) continue;
+	val[Xb + Yb*Nx] ++;
+	break;
+      case IS_SQUARE:
+	for (ix = -scale; ix <= scale; ix += dX) {
+	  for (iy = -scale; iy <= scale; iy += dY) {
+	    double rp, dp;
+	    XY_to_RD (&rp, &dp, ix, iy, &coords);
+	    while (rp < Rmin) rp += 360.0;
+	    while (rp > Rmax) rp -= 360.0;
+	    RD_to_XY (&x, &y, rp, dp, &graphmode.coords);
+	    Xb = (x - Xmin) / dX;
+	    Yb = (y - Ymin) / dY;
+	    if (Xb >= Nx) continue;
+	    if (Yb >= Ny) continue;
+	    if (Xb < 0) continue;
+	    if (Yb < 0) continue;
+	    val[Xb + Yb*Nx] += Normalize ? fSquare : 1.0;
+	  }
+	}
+	break;
+      case IS_CIRCLE:
+	for (ix = -scale; ix <= scale; ix += dX) {
+	  for (iy = -scale; iy <= scale; iy += dY) {
+	    float r2 = ix*ix + iy*iy;
+	    double rp, dp;
+	    if (r2 > scalescale) continue;
+	    XY_to_RD (&rp, &dp, ix, iy, &coords);
+	    while (rp < Rmin) rp += 360.0;
+	    while (rp > Rmax) rp -= 360.0;
+	    RD_to_XY (&x, &y, rp, dp, &graphmode.coords);
+	    Xb = (x - Xmin) / dX;
+	    Yb = (y - Ymin) / dY;
+	    if (Xb >= Nx) continue;
+	    if (Yb >= Ny) continue;
+	    if (Xb < 0) continue;
+	    if (Yb < 0) continue;
+	    val[Xb + Yb*Nx] += Normalize ? fCircle : 1.0;
+	  }
+	}
+	break;
+      case IS_GAUSS:
+	for (ix = -3.0*scale; ix <= 3.0*scale; ix += dX) {
+	  for (iy = -3.0*scale; iy <= 3.0*scale; iy += dY) {
+	    float r2 = ix*ix + iy*iy;
+	    double rp, dp;
+	    XY_to_RD (&rp, &dp, ix, iy, &coords);
+	    while (rp < Rmin) rp += 360.0;
+	    while (rp > Rmax) rp -= 360.0;
+	    RD_to_XY (&x, &y, rp, dp, &graphmode.coords);
+	    Xb = (x - Xmin) / dX;
+	    Yb = (y - Ymin) / dY;
+	    if (Xb >= Nx) continue;
+	    if (Yb >= Ny) continue;
+	    if (Xb < 0) continue;
+	    if (Yb < 0) continue;
+	    val[Xb + Yb*Nx] += fGauss*exp(-fSigma*r2);
+	  }
+	}
+	break;
+    }
+  }
+  return (TRUE);
+}
Index: /trunk/Ohana/src/opihi/cmd.astro/fitplx.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 36489)
@@ -112,5 +112,5 @@
   coords.pc1_2  = coords.pc2_1 = 0.0;
   coords.Npolyterms = 1;
-  strcpy (coords.ctype, "RA---SIN");
+  strcpy (coords.ctype, "DEC--SIN");
 
   double *X, *Y, *t, *pX, *pY, *dX, *dY;
Index: /trunk/Ohana/src/opihi/cmd.astro/fitpm.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/fitpm.c	(revision 36489)
+++ /trunk/Ohana/src/opihi/cmd.astro/fitpm.c	(revision 36489)
@@ -0,0 +1,262 @@
+# include "astro.h"
+# define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
+
+# define ESCAPE(MSG,...) {			\
+    gprint (GP_ERR, MSG, __VA_ARGS__);		\
+    return FALSE; }
+
+typedef struct {
+  double Ro, dRo;
+  double Do, dDo;
+
+  double uR, duR;
+  double uD, duD;
+
+  double chisq;
+  int Nfit;
+} PMFit;
+
+int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
+
+int fitpm (int argc, char **argv) {
+  
+  int i, N;
+
+  Vector *rvec, *dvec, *tvec, *dRvec, *dDvec;
+
+  Vector *mvec = NULL; // mask vector
+  if ((N = get_argument (argc, argv, "-mask"))) {
+    remove_argument (N, &argc, argv);
+    if ((mvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+    CastVector (mvec, OPIHI_INT);
+  }
+
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+  if ((N = get_argument (argc, argv, "-vv"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = 2;
+  }
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
+    // what about the errors?
+    return (FALSE);
+  }
+
+  /* select input / output buffers */
+  // if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  // Nx = buf[0].header.Naxis[0];
+  // Ny = buf[0].header.Naxis[1];
+  
+  if ((rvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[1]);
+  if ((dvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[3]);
+  if ((tvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[5]);
+
+  if ((dRvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[2]);
+  if ((dDvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[4]);
+
+  double *R = rvec->elements.Flt;
+  double *D = dvec->elements.Flt;
+  double *T = tvec->elements.Flt;
+  
+  double *dR = dRvec->elements.Flt;
+  double *dD = dDvec->elements.Flt;
+
+  int *mask = NULL;
+  if (mvec) {
+    mask = mvec->elements.Int;
+  }
+
+  N = tvec->Nelements; // XXX check other lengths
+
+  // find mean values to remove
+  double Npts = 0;
+  double Tmean = 0;
+  double Rmean = 0;
+  double Dmean = 0;
+  double Tmin = +1000000;
+  double Tmax = -1000000;
+  for (i = 0; i < N; i++) {
+    if (mask && !mask[i]) continue;
+    Rmean += R[i];
+    Dmean += D[i];
+    Tmean += T[i];
+    Tmin = MIN(Tmin, T[i]);
+    Tmax = MAX(Tmax, T[i]);
+    Npts += 1.0;
+  }
+  Rmean /= Npts;
+  Dmean /= Npts;
+  Tmean /= Npts;
+
+  float Trange = Tmax - Tmin;
+  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax);
+
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  Coords coords;
+  coords.crval1 = Rmean;
+  coords.crval2 = Dmean;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "DEC--SIN");
+
+  double *X, *Y, *t, *dX, *dY;
+  ALLOCATE (X, double, N);
+  ALLOCATE (Y, double, N);
+  ALLOCATE (dX, double, N);
+  ALLOCATE (dY, double, N);
+  ALLOCATE (t, double, N);
+
+  int n = 0;
+  for (i = 0; i < N; i++) {
+    if (mask && !mask[i]) continue;
+    RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords);
+    dX[n] = dR[i];
+    dY[n] = dD[i];
+    t[n] = (T[i] - Tmean) / 365.25;
+    n++;
+  }
+
+  PMFit fit;
+  if (!FitPMonly (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
+    return FALSE;
+  }
+
+  // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
+  
+  XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords);
+  if (VERBOSE) {
+    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo);
+    fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD);
+    fprintf (stderr, "chisq: %f Nfit %d\n", fit.chisq, fit.Nfit);
+  }
+
+  set_variable ("RA",   Rmean);
+  set_variable ("DEC",  Dmean);
+  set_variable ("dR",   fit.dRo);
+  set_variable ("dD",   fit.dDo);
+  set_variable ("uR",   fit.uR);
+  set_variable ("uD",   fit.uD);
+  set_variable ("duR",   fit.duR);
+  set_variable ("duD",   fit.duD);
+  set_variable ("plx",  0.0);
+  set_variable ("dplx", 0.0);
+  
+  set_variable ("Tmean",  Tmean);
+  set_variable ("Trange", Trange);
+  set_variable ("Prange", 0.0);
+
+  set_variable ("chisq", fit.chisq);
+  set_variable ("Nfit",  fit.Nfit);
+
+  return (TRUE);
+}
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
+
+  int i;
+
+  static double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  double chisq, Xf, Yf;
+
+  /* do I need to do this as 2 2x2 matrix equations? */
+  if (A == NULL) {
+    ALLOCATE (A, double *, 4);
+    ALLOCATE (B, double *, 4);
+    for (i = 0; i < 4; i++) {
+      ALLOCATE (A[i], double, 4);
+      ALLOCATE (B[i], double, 1);
+      memset (A[i], 0, 4*sizeof(double));
+      memset (B[i], 0, 1*sizeof(double));
+    }
+  }
+
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+  for (i = 0; i < Npts; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(dX[i]);
+    wy = 1.0 / SQ(dY[i]);
+
+    Wx += wx;
+    Wy += wy;
+
+    Tx += T[i]*wx;
+    Ty += T[i]*wy;
+    
+    Tx2 += SQ(T[i])*wx;
+    Ty2 += SQ(T[i])*wy;
+    
+    Xs += X[i]*wx;
+    Ys += Y[i]*wy;
+
+    XT += X[i]*T[i]*wx;
+    YT += Y[i]*T[i]*wy;
+  }
+
+  A[0][0] = Wx;
+  A[0][1] = Tx;
+
+  A[1][0] = Tx;
+  A[1][1] = Tx2;
+
+  A[2][2] = Wy;
+  A[2][3] = Ty;
+
+  A[3][2] = Ty;
+  A[3][3] = Ty2;
+
+  B[0][0] = Xs;
+  B[1][0] = XT;
+  B[2][0] = Ys;
+  B[3][0] = YT;
+
+  if (!dgaussjordan ((double **)A, (double **)B, 4, 1)) {
+    if (VERBOSE) fprintf (stderr, "error in fit\n");
+    if (VERBOSE == 2) {
+      int j;
+      for (i = 0; i < 4; i++) {
+	for (j = 0; j < 4; j++) {
+	  fprintf (stderr, "%e ", A[i][j]);
+	}
+	fprintf (stderr, " : %e\n", A[i][0]);
+      }
+    }
+    return FALSE;
+  }
+
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].duR = sqrt(A[1][1]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].duD = sqrt(A[3][3]);
+  
+  // add up the chi square for the fit
+  chisq = 0.0;
+  for (i = 0; i < Npts; i++) {
+    Xf = fit[0].Ro + fit[0].uR*T[i];
+    Yf = fit[0].Do + fit[0].uD*T[i];
+    chisq += SQ(X[i] - Xf) / SQ(dX[i]);
+    chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+    // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+  }
+  fit[0].Nfit = Npts;
+
+  // the reduced chisq is divided by (Ndof = 2*Npts - 4)
+  fit[0].chisq = chisq / (2.0*Npts - 4.0);
+  return (TRUE);
+}
Index: /trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 36489)
@@ -13,7 +13,9 @@
 int czplot                  PROTO((int, char **));
 int czcplot                 PROTO((int, char **));
+int cdensify                PROTO((int, char **));
 int drizzle                 PROTO((int, char **));
 int flux                    PROTO((int, char **));
 int fitplx                  PROTO((int, char **));
+int fitpm                   PROTO((int, char **));
 int fixwrap                 PROTO((int, char **));
 int fiximage                PROTO((int, char **));
@@ -73,7 +75,9 @@
   {1, "czplot",      czplot,       "plot scaled vectors in sky coordinates"},
   {1, "czcplot",     czcplot,      "plot color-scaled vectors in sky coordinates"},
+  {1, "cdensify",    cdensify,      "vectors to density history on projection"},
   {1, "drizzle",     drizzle,      "transform image to image"},
   {1, "flux",        flux,         "flux in a convex contour"},
   {1, "fitplx",      fitplx,       "fit proper motion and parallax"},
+  {1, "fitpm",       fitpm,        "fit proper motion only"},
   {1, "fixwrap",     fixwrap,      "fix megacam over-wrapped pixels"},
   {1, "fiximage",    fiximage,     "fix pixels in an image by interpolation"},
Index: /trunk/Ohana/src/opihi/cmd.astro/region.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.astro/region.c	(revision 36489)
@@ -80,22 +80,22 @@
   if (!ohana_str_to_radec (&Ra, &Dec, argv[1], argv[2])) return (FALSE);
   Radius = atof (argv[3]);
-  strcpy (graphmode.coords.ctype, "RA---TAN");
+  strcpy (graphmode.coords.ctype, "DEC--TAN");
   if (argc == 5) {
     if (!strcasecmp (argv[4], "TAN")) 
-      strcpy (graphmode.coords.ctype, "RA---TAN");
+      strcpy (graphmode.coords.ctype, "DEC--TAN");
     if (!strcasecmp (argv[4], "SIN")) 
-      strcpy (graphmode.coords.ctype, "RA---SIN");
+      strcpy (graphmode.coords.ctype, "DEC--SIN");
     if (!strcasecmp (argv[4], "ARC")) 
-      strcpy (graphmode.coords.ctype, "RA---ARC");
+      strcpy (graphmode.coords.ctype, "DEC--ARC");
     if (!strcasecmp (argv[4], "STG")) 
-      strcpy (graphmode.coords.ctype, "RA---STG");
+      strcpy (graphmode.coords.ctype, "DEC--STG");
     if (!strcasecmp (argv[4], "ZEA"))
-      strcpy (graphmode.coords.ctype, "RA---ZEA");
+      strcpy (graphmode.coords.ctype, "DEC--ZEA");
     if (!strcasecmp (argv[4], "AIT")) 
-      strcpy (graphmode.coords.ctype, "RA---AIT");
+      strcpy (graphmode.coords.ctype, "DEC--AIT");
     if (!strcasecmp (argv[4], "GLS")) 
-      strcpy (graphmode.coords.ctype, "RA---GLS");
+      strcpy (graphmode.coords.ctype, "DEC--GLS");
     if (!strcasecmp (argv[4], "PAR")) 
-      strcpy (graphmode.coords.ctype, "RA---PAR");
+      strcpy (graphmode.coords.ctype, "DEC--PAR");
   }
 
Index: /trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 36489)
@@ -25,4 +25,5 @@
 $(SRC)/book.$(ARCH).o		\
 $(SRC)/book_commands.$(ARCH).o	\
+$(SRC)/cast.$(ARCH).o		\
 $(SRC)/center.$(ARCH).o	\
 $(SRC)/clear.$(ARCH).o		\
Index: /trunk/Ohana/src/opihi/cmd.data/cast.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/cast.c	(revision 36489)
+++ /trunk/Ohana/src/opihi/cmd.data/cast.c	(revision 36489)
@@ -0,0 +1,71 @@
+# include "data.h"
+
+int cast (int argc, char **argv) {
+  
+  int  i, valid;
+  Vector *ivec, *ovec;
+
+  ivec = ovec = NULL;
+
+  if (argc < 6) goto usage;
+
+  valid = TRUE;
+  valid &= !strcmp(argv[2], "=");
+  valid &= !strcmp(argv[4], "as");
+  if (!valid) goto usage;
+
+  // out type 
+  char outType = OPIHI_NOTYPE;
+  if (!strcasecmp(argv[5], "int")) outType = OPIHI_INT;
+  if (!strcasecmp(argv[5], "float")) outType = OPIHI_FLT;
+  if (outType == OPIHI_NOTYPE) goto usage;
+
+  if ((ovec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto error;
+  if ((ivec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto error;
+
+  // ovec matches ivec in type
+  ResetVector (ovec, outType, ivec[0].Nelements);
+
+  // we have four cases: (ivec == flt or int) and (tvec == flt or int)
+  if ((ivec->type == OPIHI_FLT) && (ovec->type == OPIHI_FLT)) {
+    opihi_flt *vi = ivec[0].elements.Flt;
+    opihi_flt *vo = ovec[0].elements.Flt;
+    for (i = 0; i < ovec[0].Nelements; i++, vi++, vo++) {
+      *vo = *vi;
+    }
+    return (TRUE);
+  }
+  if ((ivec->type == OPIHI_FLT) && (ovec->type == OPIHI_INT)) {
+    opihi_flt *vi = ivec[0].elements.Flt;
+    opihi_int *vo = ovec[0].elements.Int;
+    for (i = 0; i < ovec[0].Nelements; i++, vi++, vo++) {
+      *vo = *vi;
+    }
+    return (TRUE);
+  }
+  if ((ivec->type == OPIHI_INT) && (ovec->type == OPIHI_FLT)) {
+    opihi_int *vi = ivec[0].elements.Int;
+    opihi_flt *vo = ovec[0].elements.Flt;
+    for (i = 0; i < ovec[0].Nelements; i++, vi++, vo++) {
+      *vo = *vi;
+    }
+    return (TRUE);
+  }
+  if ((ivec->type == OPIHI_INT) && (ovec->type == OPIHI_INT)) {
+    opihi_int *vi = ivec[0].elements.Int;
+    opihi_int *vo = ovec[0].elements.Int;
+    for (i = 0; i < ovec[0].Nelements; i++, vi++, vo++) {
+      *vo = *vi;
+    }
+    return (TRUE);
+  }
+
+error:
+  DeleteVector (ovec);
+  return (FALSE);
+
+usage:
+  gprint (GP_ERR, "SYNTAX: cast vec = vec as (int/float)\n");
+  return (FALSE);
+}
+
Index: /trunk/Ohana/src/opihi/cmd.data/densify.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/densify.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.data/densify.c	(revision 36489)
@@ -2,8 +2,9 @@
 
 # define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
+enum {IS_DOT, IS_SQUARE, IS_CIRCLE, IS_GAUSS};
 
 int densify (int argc, char **argv) {
 
-  int i, Nx, Ny, Xb, Yb, N, Xpix, Ypix, good, UseGraph;
+  int i, Nx, Ny, Xb, Yb, ix, iy, N, Xpix, Ypix, good, UseGraph;
   double Xmin, Xmax, dX, Ymin, Ymax, dY;
   float *val;
@@ -24,8 +25,26 @@
   }
 
+  float scale = 0.0;
+  if ((N = get_argument (argc, argv, "-scale"))) {
+    remove_argument (N, &argc, argv);
+    scale = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int PSFTYPE = IS_DOT;
+  if ((N = get_argument (argc, argv, "-psf"))) {
+    remove_argument (N, &argc, argv);
+    if (!strcasecmp(argv[N], "dot"))    PSFTYPE = IS_DOT;
+    if (!strcasecmp(argv[N], "square")) PSFTYPE = IS_SQUARE;
+    if (!strcasecmp(argv[N], "circle")) PSFTYPE = IS_CIRCLE;
+    if (!strcasecmp(argv[N], "gauss"))  PSFTYPE = IS_GAUSS;
+    remove_argument (N, &argc, argv);
+  }
+
   good = UseGraph ? (argc == 4) : (argc == 10);
   if (!good) {
     gprint (GP_ERR, "USAGE: densify buffer x y Xmin Xmax dX Ymin Ymax dY\n");
     gprint (GP_ERR, "   OR: densify buffer x y -graph\n");
+    gprint (GP_ERR, " option: -psf [dot] (circle) (square) (gauss)\n");
     return (FALSE);
   }
@@ -69,4 +88,7 @@
   CHECKVAL(dY);
 
+  float scaleX = (scale > 0.0) ? scale / dX : 3.0;
+  float scaleY = (scale > 0.0) ? scale / dY : 3.0;
+
   Nx = (Xmax - Xmin) / dX + 1;
   Ny = (Ymax - Ymin) / dY + 1;
@@ -76,4 +98,10 @@
   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
   strcpy (bf[0].file, "(empty)");
+  
+  float scale2 = (scaleX + 1.0) * (scaleY + 1.0);
+  float fSquare = 1.0 / scale2;
+  float fCircle = 1.0 / (3.141592 * scale2);
+  float fSigma  = 0.5 / scale2;
+  float fGauss  = 1.0 / (2.0 * 3.141592 * scale2);
 
   x = vx[0].elements.Flt;
@@ -83,9 +111,53 @@
     Xb = (*x - Xmin) / dX;
     Yb = (*y - Ymin) / dY;
-    if (Xb >= Nx) continue;
-    if (Yb >= Ny) continue;
-    if (Xb < 0) continue;
-    if (Yb < 0) continue;
-    val[Xb + Yb*Nx] ++;
+    switch (PSFTYPE) {
+      case IS_DOT:
+	if (Xb >= Nx) continue;
+	if (Yb >= Ny) continue;
+	if (Xb < 0) continue;
+	if (Yb < 0) continue;
+	val[Xb + Yb*Nx] ++;
+	break;
+      case IS_SQUARE:
+	for (ix = Xb - scaleX; ix <= Xb + scaleX; ix++) {
+	  for (iy = Yb - scaleY; iy <= Yb + scaleY; iy++) {
+	    if (ix >= Nx) continue;
+	    if (iy >= Ny) continue;
+	    if (ix < 0) continue;
+	    if (iy < 0) continue;
+	    val[ix + iy*Nx] += fSquare;
+	  }
+	}
+	break;
+      case IS_CIRCLE:
+	for (ix = Xb - scaleX; ix <= Xb + scaleX; ix++) {
+	  float dX = ix - Xb;
+	  for (iy = Yb - scaleY; iy <= Yb + scaleY; iy++) {
+	    float dY = iy - Yb;
+	    float r2 = dX*dX + dY*dY;
+	    if (r2 > 9) continue;
+	    if (ix >= Nx) continue;
+	    if (iy >= Ny) continue;
+	    if (ix < 0) continue;
+	    if (iy < 0) continue;
+	    val[ix + iy*Nx] += fCircle;
+	  }
+	}
+	break;
+      case IS_GAUSS:
+	for (ix = Xb - scaleX; ix <= Xb + scaleX; ix++) {
+	  float dX = ix - Xb;
+	  for (iy = Yb - scaleY; iy <= Yb + scaleY; iy++) {
+	    float dY = iy - Yb;
+	    float r2 = dX*dX + dY*dY;
+	    if (ix >= Nx) continue;
+	    if (iy >= Ny) continue;
+	    if (ix < 0) continue;
+	    if (iy < 0) continue;
+	    val[ix + iy*Nx] += fGauss*exp(-fSigma*r2);
+	  }
+	}
+	break;
+    }
   }
   return (TRUE);
Index: /trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/init.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.data/init.c	(revision 36489)
@@ -10,4 +10,5 @@
 int center           PROTO((int, char **));
 int parity           PROTO((int, char **));
+int cast             PROTO((int, char **));
 int circstats        PROTO((int, char **));
 int clear            PROTO((int, char **));
@@ -164,4 +165,5 @@
   {1, "buffers",      list_buffers,     "list the currently allocated buffers (images)"},
   {1, "center",       center,           "center image on coords"},
+  {1, "cast",         cast,             "cast input vector to specified type"},
   {1, "circstats",    circstats,        "circular statistics"},
   {1, "clear",        clear,            "erase plot"},
Index: /trunk/Ohana/src/opihi/cmd.data/vgauss.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/vgauss.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/cmd.data/vgauss.c	(revision 36489)
@@ -48,4 +48,7 @@
   if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
+  CastVector (xvec, OPIHI_FLT);
+  CastVector (yvec, OPIHI_FLT);
+
   int Nsvec = strlen(argv[3]);
 
@@ -69,6 +72,4 @@
   }
 
-  CastVector (xvec, OPIHI_FLT);
-  CastVector (yvec, OPIHI_FLT);
   CastVector (svec, OPIHI_FLT);
   // XXX Cast is failing.
Index: /trunk/Ohana/src/opihi/dvo/find_matches.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/find_matches.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/dvo/find_matches.c	(revision 36489)
@@ -50,5 +50,5 @@
   tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
   tcoords.Npolyterms = 1;
-  strcpy (tcoords.ctype, "RA---ARC");
+  strcpy (tcoords.ctype, "DEC--ARC");
 
   // this region includes a boundary layer of size RADIUS
Index: /trunk/Ohana/src/opihi/dvo/gimages.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/gimages.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/dvo/gimages.c	(revision 36489)
@@ -136,5 +136,5 @@
       local.pc2_1 = local.pc1_2 = 0.0;
       local.Npolyterms = 1;
-      strcpy (local.ctype, "RA---TAN");
+      strcpy (local.ctype, "DEC--TAN");
 
       if (typehash == DistortImage) {
Index: /trunk/Ohana/src/opihi/dvo/showtile.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/showtile.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/dvo/showtile.c	(revision 36489)
@@ -34,5 +34,5 @@
   coords.pc1_2  = coords.pc2_1  = 0.0;
   coords.Npolyterms = 0;
-  strcpy (coords.ctype, "RA---TAN");
+  strcpy (coords.ctype, "DEC--TAN");
   
   /* fill in top-left region */
Index: /trunk/Ohana/src/opihi/dvo/simage.c
===================================================================
--- /trunk/Ohana/src/opihi/dvo/simage.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/dvo/simage.c	(revision 36489)
@@ -43,6 +43,6 @@
   /* get astrometry information */
   strcpy (coords.ctype, "NONE");
-  gfits_scan (&header, "CTYPE1",   "%s",  1, coords.ctype);
-  if (strcmp (coords.ctype, "RA---PLY")) {
+  gfits_scan (&header, "CTYPE2",   "%s",  1, coords.ctype);
+  if (strcmp (coords.ctype, "DEC--PLY")) {
     gprint (GP_ERR, "ERROR: wrong astrometric info in header\n");
     return (FALSE);
Index: /trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 36489)
@@ -30,4 +30,8 @@
     if (!strcmp (argv[i], "abs"))    { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "int"))    { type = ST_UNARY; goto gotit; }
+    if (!strcmp (argv[i], "floor"))  { type = ST_UNARY; goto gotit; }
+    if (!strcmp (argv[i], "round"))  { type = ST_UNARY; goto gotit; }
+    if (!strcmp (argv[i], "ceil"))   { type = ST_UNARY; goto gotit; }
+    if (!strcmp (argv[i], "rint"))   { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "exp"))    { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "ten"))    { type = ST_UNARY; goto gotit; }
Index: /trunk/Ohana/src/opihi/lib.shell/stack_math.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 36488)
+++ /trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 36489)
@@ -1123,4 +1123,7 @@
   if (!strcmp (op, "abs"))    S_FUNC(fabs(M1), ST_SCALAR_INT);
   if (!strcmp (op, "int"))    S_FUNC((long long)(M1), ST_SCALAR_INT);
+  if (!strcmp (op, "floor"))  S_FUNC(floor (M1), ST_SCALAR_FLT);
+  if (!strcmp (op, "ceil"))   S_FUNC(ceil (M1), ST_SCALAR_FLT);
+  if (!strcmp (op, "rint"))   S_FUNC(nearbyint (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "exp"))    S_FUNC(exp (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "ten"))    S_FUNC(pow (10.0,M1), ST_SCALAR_FLT);
@@ -1202,4 +1205,7 @@
   if (!strcmp (op, "abs"))    V_FUNC(fabs(*M1), ST_SCALAR_INT);
   if (!strcmp (op, "int"))    V_FUNC((long long)(*M1), ST_SCALAR_INT);
+  if (!strcmp (op, "floor"))  V_FUNC(floor (*M1), ST_SCALAR_FLT);
+  if (!strcmp (op, "ceil"))   V_FUNC(ceil (*M1), ST_SCALAR_FLT);
+  if (!strcmp (op, "rint"))   V_FUNC(nearbyint (*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "exp"))    V_FUNC(exp(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "ten"))    V_FUNC(pow(10.0,*M1), ST_SCALAR_FLT);
@@ -1273,4 +1279,9 @@
   if (!strcmp (op, "abs"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = fabs(*M1);         }}
   if (!strcmp (op, "int"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = (opihi_flt)(long long)(*M1); }}
+
+  if (!strcmp (op, "floor")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = floor (*M1); }}
+  if (!strcmp (op, "ceil"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = ceil (*M1); }}
+  if (!strcmp (op, "rint"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = nearbyint (*M1); }}
+
   if (!strcmp (op, "exp"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = exp(*M1);          }}
   if (!strcmp (op, "ten"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = pow(10.0,*M1);     }}