Index: trunk/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 39233)
+++ trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 39294)
@@ -11,5 +11,5 @@
   Buffer *bf;
   Vector *vr, *vd;
-  opihi_flt *r, *d, x, y;
+  opihi_flt x, y;
   int kapa;
   Graphdata graphmode;
@@ -24,4 +24,11 @@
   double Rmin = graphmode.coords.crval1 - 182.0;
   double Rmax = graphmode.coords.crval1 + 182.0;
+
+  Vector *vv = NULL;
+  if ((N = get_argument (argc, argv, "-value"))) {
+    remove_argument (N, &argc, argv);
+    if ((vv = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
 
   float scale = 0.0;
@@ -45,4 +52,8 @@
     gprint (GP_ERR, "USAGE: cdensify buffer R D\n");
     gprint (GP_ERR, " option: -psf [dot] (circle) (square) (gauss)\n");
+    gprint (GP_ERR, " other options:\n");
+    gprint (GP_ERR, "   -raw : do not renormalize PSF\n");
+    gprint (GP_ERR, "   -scale : spatial scale factor for PSF (radius, half-width, or sigma)\n");
+    gprint (GP_ERR, "   -value (vector) : multiply sum by the given vector\n");
     return (FALSE);
   }
@@ -53,4 +64,11 @@
 
   if (vr[0].Nelements != vd[0].Nelements) return (FALSE);
+
+  if (vv) {
+    if (vv[0].Nelements != vr[0].Nelements) {
+      gprint (GP_ERR, "mis-match in vector lengths\n");
+      return FALSE;
+    }
+  }
 
   REQUIRE_VECTOR_FLT (vr, FALSE); 
@@ -99,7 +117,13 @@
   InitCoords (&coords, "DEC--TAN");
 
-  r = vr[0].elements.Flt;
-  d = vd[0].elements.Flt;
+  opihi_flt *r = vr[0].elements.Flt;
+  opihi_flt *d = vd[0].elements.Flt;
+
+  opihi_flt *Fs = vv ? vv[0].elements.Flt : NULL;
+  opihi_int *Is = vv ? vv[0].elements.Int : NULL;
+  int isFloatScale = (vv && vv[0].type == OPIHI_FLT);
+
   val = (float *)bf[0].matrix.buffer;
+
   for (i = 0; i < vr[0].Nelements; i++, r++, d++) {
     double rn = ohana_normalize_angle (*r);
@@ -108,4 +132,9 @@
     coords.crval1 = rn;
     coords.crval2 = *d;
+
+    float F = 1.0;
+    if (vv) {
+      F = isFloatScale ? Fs[i] : Is[i];
+    }
 
     switch (PSFTYPE) {
@@ -118,5 +147,9 @@
 	if (Xb < 0) continue;
 	if (Yb < 0) continue;
-	val[Xb + Yb*Nx] ++;
+	if (vv) {
+	  val[Xb + Yb*Nx] += F;
+	} else {
+	  val[Xb + Yb*Nx] ++;
+	}
 	break;
       case IS_SQUARE:
@@ -134,5 +167,9 @@
 	    if (Xb < 0) continue;
 	    if (Yb < 0) continue;
-	    val[Xb + Yb*Nx] += Normalize ? fSquare : 1.0;
+	    if (vv) {
+	      val[Xb + Yb*Nx] += Normalize ? fSquare*F : F;
+	    } else {
+	      val[Xb + Yb*Nx] += Normalize ? fSquare : 1.0;
+	    }
 	  }
 	}
@@ -154,5 +191,9 @@
 	    if (Xb < 0) continue;
 	    if (Yb < 0) continue;
-	    val[Xb + Yb*Nx] += Normalize ? fCircle : 1.0;
+	    if (vv) {
+	      val[Xb + Yb*Nx] += Normalize ? fCircle*F : F;
+	    } else {
+	      val[Xb + Yb*Nx] += Normalize ? fCircle : 1.0;
+	    }
 	  }
 	}
@@ -173,5 +214,9 @@
 	    if (Xb < 0) continue;
 	    if (Yb < 0) continue;
-	    val[Xb + Yb*Nx] += fGauss*exp(-fSigma*r2);
+	    if (vv) {
+	      val[Xb + Yb*Nx] += F*fGauss*exp(-fSigma*r2);
+	    } else {
+	      val[Xb + Yb*Nx] += fGauss*exp(-fSigma*r2);
+	    }
 	  }
 	}