Index: /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/Makefile	(revision 36386)
+++ /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/Makefile	(revision 36387)
@@ -28,4 +28,5 @@
 $(SRC)/cval.$(ARCH).o		   \
 $(SRC)/czplot.$(ARCH).o	   \
+$(SRC)/cdensify.$(ARCH).o	   \
 $(SRC)/drizzle.$(ARCH).o	   \
 $(SRC)/flux.$(ARCH).o		   \
Index: /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 36387)
+++ /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 36387)
@@ -0,0 +1,181 @@
+# 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;
+  
+  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)");
+  
+  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, "RA---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;
+	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: /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/init.c	(revision 36386)
+++ /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.astro/init.c	(revision 36387)
@@ -13,4 +13,5 @@
 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 **));
@@ -73,4 +74,5 @@
   {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"},
Index: /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.data/densify.c
===================================================================
--- /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.data/densify.c	(revision 36386)
+++ /branches/eam_branches/ipp-20131211/Ohana/src/opihi/cmd.data/densify.c	(revision 36387)
@@ -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);