Index: /branches/eam_branches/ipp-20230313/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- /branches/eam_branches/ipp-20230313/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 42583)
+++ /branches/eam_branches/ipp-20230313/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 42584)
@@ -1,3 +1,4 @@
 # include "data.h"
+# define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
 
 # define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
@@ -116,4 +117,7 @@
   PutCoords (&newcoords, &bf[0].header);
   
+  // use the mask to prevent double-counting
+  ALLOCATE_PTR (mask, char, Nx*Ny); 
+
   float scalescale = scale*scale;
   float scale2 = (scale + 1.0) * (scale + 1.0);
@@ -142,7 +146,13 @@
     coords.crval1 = rn;
     coords.crval2 = *d;
+    
+    // XXX do not oversample by more than a factor of 10
+    float dXn = dX * MAX(fabs(dCOS(*d)), 0.1);
 
     float F = 1.0;
     if (vv) { F = isFloatScale ? Fs[i] : Is[i]; }
+
+    // reset the mask so we do not double-count
+    memset (mask, 0, Nx*Ny);
 
     switch (PSFTYPE) {
@@ -162,5 +172,5 @@
 	break;
       case IS_SQUARE:
-	for (ix = -scale; ix <= scale; ix += dX) {
+	for (ix = -scale; ix <= scale; ix += dXn) {
 	  for (iy = -scale; iy <= scale; iy += dY) {
 	    double rp, dp;
@@ -184,6 +194,6 @@
 	break;
       case IS_CIRCLE:
-	for (ix = -scale; ix <= scale; ix += dX) {
-	  for (iy = -scale; iy <= scale; iy += dY) {
+	for (ix = -scale - dXn; ix <= scale + dXn; ix += dXn) {
+	  for (iy = -scale - dY; iy <= scale + dY; iy += dY) {
 	    float r2 = ix*ix + iy*iy;
 	    double rp, dp;
@@ -199,4 +209,6 @@
 	    if (Xb < 0) continue;
 	    if (Yb < 0) continue;
+	    if (mask[Xb + Yb*Nx]) continue;
+	    mask[Xb + Yb*Nx] = 1;
 	    if (vv) {
 	      val[Xb + Yb*Nx] += Normalize ? fCircle*F : F;
@@ -208,5 +220,5 @@
 	break;
       case IS_GAUSS:
-	for (ix = -3.0*scale; ix <= 3.0*scale; ix += dX) {
+	for (ix = -3.0*scale; ix <= 3.0*scale; ix += dXn) {
 	  for (iy = -3.0*scale; iy <= 3.0*scale; iy += dY) {
 	    float r2 = ix*ix + iy*iy;
@@ -232,4 +244,6 @@
     }
   }
+
+  free (mask);
   return (TRUE);
 }