Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c	(revision 42392)
@@ -23,4 +23,5 @@
 /* Find all the burned patches and bright stars */
    err = star_detect(nx, ny, NX, NY, imbuf, mbuf, msbuf, cell, cellnum);
+
    err = burn_check(nx, ny, NX, NY, imbuf, mbuf, cell);
 
@@ -34,5 +35,5 @@
 //   fprintf(stderr, "Got through grow mask for burns\n");
    err = grow_mask(nx, ny, NX, mbuf, BMASK_GROW, RMASK_GROW, 
-		   MASK_STAR_HALO, cell->nstar, cell->star);
+         MASK_STAR_HALO, cell->nstar, cell->star);
 //   fprintf(stderr, "Got through grow mask for stars\n");
 
@@ -70,6 +71,7 @@
    printf("  #    cx  cy  max     sx  sy    ex  ey  midy  diff  S/B/E  slope     zero\n");
    for(k=0; k<cell->nburn; k++) {
-      ymid = (cell->burn[k].y0m + cell->burn[k].y1m +
-	      cell->burn[k].y0p + cell->burn[k].y1p + 2) / 4;
+      /*ymid = (cell->burn[k].y0m + cell->burn[k].y1m +
+	      cell->burn[k].y0p + cell->burn[k].y1p + 2) / 4;*/
+      ymid = (cell->burn[k].sy + cell->burn[k].ey + 1) / 2;
       i = cell->burn[k].nfit / 2;
       printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %5d %2d %1d %1d %6.3f %8d\n", 
@@ -88,6 +90,7 @@
    printf("  #    cx  cy  max     sx  sy    ex  ey  midy  S/P\n");
    for(k=0; k<cell->nstar; k++) {
-      ymid = (cell->star[k].y0m + cell->star[k].y1m +
-	      cell->star[k].y0p + cell->star[k].y1p + 2) / 4;
+      /*ymid = (cell->star[k].y0m + cell->star[k].y1m +
+	      cell->star[k].y0p + cell->star[k].y1p + 2) / 4;*/
+      ymid = (cell->star[k].sy + cell->star[k].ey + 1) / 2;
 
       printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %1d\n", 
@@ -142,5 +145,5 @@
 		       MTYPE *mask, CELL *cell)
 {
-   int err, k;
+   int err,k;
 
 /* Given these big stars, identify which ones really have a burn */
@@ -160,5 +163,5 @@
 	     k, cell->burn[k].sx, cell->burn[k].sy, 
 	     cell->burn[k].ex, cell->burn[k].ey, 
-	     cell->burn[k].y0m, cell->burn[k].y1m, cell->burn[k].y0p, cell->burn[k].y0p,
+	     cell->burn[k].y0m, cell->burn[k].y1m, cell->burn[k].y0p, cell->burn[k].y1p,
 	     cell->burn[k].diff, cell->burn[k].burned);
    }
@@ -242,9 +245,20 @@
    int y0, y1, ymid, xburn, xref, dx=0, dy=0;
 
+   /*introduce a maximum to the burn contrast check, to avoid a very large number in crowded regions*/
+   int burn_contrast = MIN(0.3*rms,5); 
+   /*introduce a maximum size for the burn contrast check box, instead of a random large 100 pixels */
+   /*for small stars you will lose signal after 10 pixels, so use the fitted box size as a good guide*/
+   int test_box = MAX(box->ey-box->sy,10); 
+
+   if(VERBOSE & VERB_BOXGROW) {
+	printf("burn contrast test using box of max size: %8d\n", test_box);
+   }
+
 /* FIXME: needs a test for flat-toppedness as well as trail... */
 /* FIXME: needs a test for pure, massive saturation */
 
 /* Center line */
-   ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+   //ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+   ymid = (box->sy+box->ey)/2;
 
 /* Irrelevant, too near the top to detect a burn anyway */
@@ -283,5 +297,5 @@
       y0 = ymid + dy;
       y1 = y0 + MIN(ny-1-y0, 2*ymid-y0-1);
-      if(y1-y0 > 100) y1 = y0 + 100;
+      if(y1-y0 > test_box) y1 = y0 + test_box;
    }
 
@@ -292,26 +306,30 @@
       nburn = nref = 0;
       for(i=box->sx; i<=box->ex; i++) {
-	 if(mask[i+j*NX] == MASK_NONE) {
-	    xburn += data[i+j*NX];
-	    nburn++;
-	 }
-	 if(dx != 0) {	/* Side to side */
-	    if(mask[i+dx+j*NX] == MASK_NONE) {
-	       xref += data[i+dx+j*NX];
-	       nref++;
-	    }
-	 } else {	/* top-bottom */
-	    if(mask[i+(2*ymid-j)*NX] == MASK_NONE) {
-	       xref += data[i+(2*ymid-j)*NX];
-	       nref++;
-	    }
-	 }
+	     if(mask[i+j*NX] == MASK_NONE) {
+	       xburn += data[i+j*NX];
+	       nburn++;
+	     }
+	     if(dx != 0) {	/* Side to side */
+	       if(mask[i+dx+j*NX] == MASK_NONE) {
+	        xref += data[i+dx+j*NX];
+	        nref++;
+	       }
+	     } else {	/* top-bottom */
+	       if(mask[i+(2*ymid-j)*NX] == MASK_NONE) {
+	        xref += data[i+(2*ymid-j)*NX];
+	        nref++;
+	       }
+	     }
       }
       if(nburn > 0) xburn /= nburn;
       if(nref > 0) xref /= nref;
       if(nref > 0 && nburn > 0) median_buf[nmed++] = xburn - xref;
+      if(VERBOSE & VERB_BOXGROW) {
+	printf("contrast: %8d %8d %8d %8d %8d\n", j, xburn, (2*ymid-j), xref, xburn - xref);
+      }
+
    }
    box->diff = int_median(nmed, median_buf);
-   box->burned = (box->diff > 0.3*rms);
+   box->burned = (box->diff > burn_contrast);
    return(0);
 }
@@ -335,5 +353,6 @@
 
 /* Center line */
-   ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+   //ymid = (box->y0m + box->y0p + box->y1m + box->y1p + 2) / 4;
+   ymid = (box->sy + box->ey + 1) / 2;
 
 /* Starting point offset */
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c	(revision 42392)
@@ -89,6 +89,6 @@
    PSF_THRESH   =  5000;	/* Threshold for a star to be a PSF */
 
-   BMASK_GROW =  5.0;		/* Growth of burned boxes in size */
-   RMASK_GROW =  3.0;		/* Growth of burn/star in diameter */
+   BMASK_GROW =  3.0;		/* Growth of burned boxes in pixels */
+   RMASK_GROW =  2.0;		/* Growth of burn/star in diameter */
 
    MIN_PSF_SIZE  =    3;	/* Min box size for stamp selection */
@@ -97,5 +97,5 @@
    MAX_PSF_PER_CELL = 20; 	/* Max number of PSF stars accepted per cell */
 
-   NEGLIGIBLE_TRAIL = 0.5;	/* Don't sweat less than this * sigma */
+   NEGLIGIBLE_TRAIL = 0.4;	/* Don't sweat less than this * sigma */
    EXPIRE_TRAIL_TIME = 2000;	/* Expire a persist after this [sec] */
 
@@ -564,13 +564,13 @@
 
 /* Does this cell look kosher? */
-	 if(OTA[cell].rms*OTA[cell].rms > 
-	    OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
-	    if(VERBOSE > 0) {
-	       fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
-	       fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n", 
-		       OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
-	    }
-	    continue;
-	 }
+//	 if(OTA[cell].rms*OTA[cell].rms > 
+//	    OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
+//	    if(VERBOSE > 0) {
+//	       fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
+//	       fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n", 
+//		       OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
+//	    }
+//	    continue;
+//	 }
 
 	 if(apply) {
@@ -601,12 +601,9 @@
 	    if(VERBOSE & VERB_NORM) burn_blab(OTA+cell);
 
-/* Fix up the streaks (don't bother if table only) */
-// Why the heck do I get to skip this if table only?  Makes no sense!
-//	    if(!tableonly) {
-	       persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, 
+	    persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, 
 			   OTA+cell);
 /* Tell us about it? */
-	       if(VERBOSE & VERB_NORM) persist_blab(OTA+cell);
-//	    }
+	    if(VERBOSE & VERB_NORM) persist_blab(OTA+cell);
+
 	 }
 
@@ -667,5 +664,5 @@
 {
 /* Make some space for medians */
-   if(2*SKY_MARG*ny > nmedian_buf) {
+/*   if(2*SKY_MARG*ny > nmedian_buf) {
       if(median_buf != NULL) free(median_buf);
       if( (median_buf = (int *)calloc(2*SKY_MARG*ny, sizeof(int))) == NULL) {
@@ -675,4 +672,15 @@
 
       nmedian_buf = 2*SKY_MARG*ny;
+   }
+*/
+
+/* Make some space for cell copy */
+   if(nx*ny > nmedian_buf) {
+      if(median_buf != NULL) free(median_buf);
+      if( (median_buf = (int *)calloc(nx*ny, sizeof(DTYPE))) == NULL) {
+       fprintf(stderr, "\rerror: failed to alloc cell copy\n");
+       exit(-668);
+      }
+      nmedian_buf = nx*ny;
    }
 
@@ -728,4 +736,6 @@
 /* Get sky stats */
    for(k=n=nsat=0; k<nx*ny; k+=((617*nx)/1000)) {
+//   for(k=n=nsat=0; k<nx*ny; k+=((100*nx)/1000)) {
+//   for(k=n=nsat=0; k<nx*ny; k++) {
       i = k % nx;
       j = k / nx;
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h	(revision 42392)
@@ -70,5 +70,5 @@
 #define STAR_RADIUS  4		/* Radius over which a star ctr must be max */
 #define SKY_MARG 3		/* Horiz offset from burn area for sky */
-#define FIT_EDGE 10		/* How far beyond saturation to start fit? */
+#define FIT_EDGE 5		/* How far beyond saturation to start fit? */
 #define Y_SCALE 0.01		/* Scale factor for y in fits */
 
@@ -102,5 +102,5 @@
 
 /* Fit parameters */
-#define FIT_MIN_SLOPE  -5.0	/* minimum slope which is a credible fit */
+#define FIT_MIN_SLOPE -10.0	/* minimum slope which is a credible fit */
 #define FIT_MAX_SLOPE   0.0	/* maximum slope which is a credible fit */
 
@@ -144,4 +144,5 @@
       int sxfit;	/* starting column for fits */
       int exfit;	/* ending column for fits */
+      int eyfit;  /* y-coord ending column for fits */
       int nfit;		/* how many columns were corrected? */
       IMTYPE *stamp;	/* postage stamp of this object */
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnutils.c	(revision 42392)
@@ -18,5 +18,5 @@
 STATIC int grow_box(int nx, int ny, int NX, DTYPE *data, int thr, OBJBOX *box)
 {
-   int i, j, is, ie, left, right;
+   int i, j,i2,j2, is, ie,jlo,jhi,ilo,ihi, left, right;
 
    left = right = 1;
@@ -27,57 +27,104 @@
 
    do {
-/* Push the right wall until it's clear */
+      /* Push the right wall until it's clear */
+      jlo = box->sy;
+      jhi = box->ey;
       if(right) {
-	 for(i=box->ex+1; i<nx; i++) {
-	    for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
-	    if(j > box->ey) break;
-	    if(VERBOSE & VERB_BOXGROW) {
-	       printf("right: %8d %8d\n", i, j);
-	    }
-	 }
-	 box->ex = i - 1;
-	 right = 0;
-      }
-
-/* Push the left wall until it's clear */
+	     for(i=box->ex+1; i<nx; i++) {
+	       for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
+          if(j > box->ey) break;
+          for(j2=box->ey; j2>=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
+          if(j2 < box->sy) break;
+
+          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+          /*If not, there is a discontinuity, and we might have overlapping sources*/
+          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
+
+	       if(VERBOSE & VERB_BOXGROW) {
+	         printf("right: %8d %8d %8d\n", i, j,data[i+j*NX]);
+	       }
+          jlo = j;
+          jhi = j2;
+	     }
+	     box->ex = i - 1;
+	     right = 0;
+      }
+
+      /* Push the left wall until it's clear */
+      jlo = box->sy;
+      jhi = box->ey;
       if(left) {
-	 for(i=box->sx-1; i>=0; i--) {
-	    for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
-	    if(j > box->ey) break;
-	    if(VERBOSE & VERB_BOXGROW) {
-	       printf("left: %8d %8d\n", i, j);
-	    }
-	 }
-	 box->sx = i + 1;
-	 left = 0;
-      }
-
-/* Push the top wall until it's clear, catching corners! */
+	     for(i=box->sx-1; i>=0; i--) {
+	       for(j=box->sy; j<=box->ey; j++) if(data[i+j*NX] >= thr) break;
+	       if(j > box->ey) break;
+          for(j2=box->ey; j2<=box->sy; j2--) if(data[i+j2*NX] >= thr) break;
+          if(j2 < box->sy) break;
+
+          /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+          /*If not, there is a discontinuity, and we might have overlapping sources*/
+          if(!(j >= jlo-0.1*(jhi-jlo) && j2 <= jhi+0.1*(jhi-jlo))) {break;}
+
+	       if(VERBOSE & VERB_BOXGROW) {
+	         printf("left: %8d %8d %8d\n", i, j,data[i+j*NX]);
+	       }
+          jlo = j;
+          jhi = j2;
+	     }
+	     box->sx = i + 1;
+	     left = 0;
+      }
+
+      /* Push the top wall until it's clear, catching corners! */
       is = MAX(0, box->sx-1);
       ie = MIN(nx-1, box->ex+1);
+      ilo = is;
+      ihi = ie;
       for(j=box->ey+1; j<ny; j++) {
-	 for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
-	 if(i == is && box->sx>0) left = 1;	/* Need to redo left */
-	 if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
-	 if(i > ie) break;			/* All clear */
-	 if(VERBOSE & VERB_BOXGROW) {
-	    printf("top: %8d %8d\n", i, j);
-	 }
+	     for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
+        if(i > ie) break;        /* All clear */
+        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
+        if(i2 < is) break;        /* All clear */
+
+	     if(i == is && box->sx>0) left = 1;	/* First point hits thr, so need to redo left */
+	     if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /*Last point hits thr, so redo right */
+
+        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+        /*If not, there is a discontinuity, and we might have overlapping sources*/
+        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) {break;}
+
+	     if(VERBOSE & VERB_BOXGROW) {
+	       printf("top: %8d %8d %8d\n", i, j,data[i+j*NX]);
+	     }
+        ilo = i;
+        ihi = i2;
       }
       box->ey = j - 1;
 
-/* Push the bottom wall until it's clear, catching corners! */
+      /* Push the bottom wall until it's clear, catching corners! */
       is = MAX(0,box->sx-1);
       ie = MIN(nx-1,box->ex+1);
+      ilo = is;
+      ihi = ie;
       for(j=box->sy-1; j>=0; j--) {
-	 for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
-	 if(i == is && box->sx>0) left = 1;	/* Need to redo left */
-	 if(data[is+1+j*NX] >= thr) box->y0m = j;
-	 if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
-	 if(data[ie-1+j*NX] >= thr) box->y1m = j;
-	 if(i > ie) break;			/* All clear */
-	 if(VERBOSE & VERB_BOXGROW) {
-	    printf("bottom: %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m);
-	 }
+	     for(i=is; i<=ie; i++) if(data[i+j*NX] >= thr) break;
+        if(i > ie) break;        /* All clear */
+        for(i2=ie; i2<=is; i--) if(data[i2+j*NX] >= thr) break;
+        if(i2 < is) break;        /* All clear */
+
+	     if(i == is && box->sx>0) left = 1;	/* Need to redo left */
+        if(data[ie+j*NX] >= thr && box->ex<nx-1) right = 1;  /* Redo right */
+
+	     if(data[is+1+j*NX] >= thr) box->y0m = j;
+	     if(data[ie-1+j*NX] >= thr) box->y1m = j;
+
+        /*Since stars are roughly diamonds, the new vertices should be contained within the previous ones*/
+        /*If not, there is a discontinuity, and we might have overlapping sources*/
+        if(!(i >= ilo-0.1*(ihi-ilo) && i2 <= ihi+0.1*(ihi-ilo))) break;
+
+	     if(VERBOSE & VERB_BOXGROW) {
+	       printf("bottom: %8d %8d %8d %8d %8d\n", i, j, box->y0m, box->y1m,data[i+j*NX]);
+	     }
+        ilo = i;
+        ihi = i2;
       }
       box->sy = j + 1;
@@ -126,5 +173,5 @@
 		     double rfac, int maskval, int nbox, OBJBOX *box)
 {
-   int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1;
+   int i, j, k, x0, y0, db, dr, r2, i0, i1, j0, j1,boxsize;
 
    for(k=0; k<nbox; k++) {
@@ -137,20 +184,22 @@
 
       for(j=j0; j<=j1; j++) {
-	 for(i=i0+1; i<i1; i++) {
-	    if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
-	    if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
-	 }
+	     for(i=i0+1; i<i1; i++) {
+	      if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
+	      if(mask[i+j*NX] == MASK_NONE) mask[i+j*NX] = maskval;
+	     }
       }
 
 /* Push out the center to a diameter rfac*xsize */
       x0 = (box[k].sx + box[k].ex + 1) / 2;
-      y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
-      dr = NINT((box[k].ex - box[k].sx + 1) * 0.5 * rfac);
+      //y0 = (box[k].y0p + box[k].y1p + box[k].y0m + box[k].y1m + 2) / 4;
+      y0 = (box[k].sy + box[k].ey + 1) / 2;
+      boxsize = MAX((box[k].ex - box[k].sx + 1),(box[k].ey - box[k].sy + 1));
+      dr = NINT(boxsize * 0.5 * rfac);
       for(j=MAX(0,y0-dr); j<=MIN(ny-1, y0+dr); j++) {
-	 for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
-	    r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
-	    if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
+	     for(i=MAX(0,x0-dr); i<=MIN(nx-1, x0+dr); i++) {
+	       r2 = (i-x0)*(i-x0) + (j-y0)*(j-y0);
+	       if(r2 < dr*dr && mask[i+j*NX] == MASK_NONE)
 	       mask[i+j*NX] = maskval;
-	 }
+	     }
       }
    }
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistfix.c	(revision 42392)
@@ -23,29 +23,29 @@
 
 /* Merge all overlapping streaks */
-   persist_merge(cell);
+  persist_merge(cell);
+
 
 /* Fix up all the persistence streaks */
    for(k=0; k<cell->npersist; k++) {
-
-/* Is this just a blasted area being carried for IPP? */
-      if((cell->persist)[k].func == BURN_BLASTED) {
-	 if(cell->time - (cell->persist)[k].time > EXPIRE_TRAIL_TIME) {
+     /* Is this just a blasted area being carried for IPP? */
+     if((cell->persist)[k].func == BURN_BLASTED) {
+	   if(cell->time - (cell->persist)[k].time > EXPIRE_TRAIL_TIME) {
 	    (cell->persist)[k].fiterr = FIT_EXPIRED;
-	 }
-	 continue;
-      }
+	   }
+	   continue;
+     }
 
 /* This had a significant positive slope, what do we do now? */
-      if( (cell->persist)[k].func == BURN_POSSLOPE) {
+     if( (cell->persist)[k].func == BURN_POSSLOPE) {
 /* try again, so let it slide through... */
-      }
+     }
 
 /* Fit the trail */
-      err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0, 
+     err = fit_trail(nx, ny, NX, imbuf, mbuf, cell->persist+k, 0, 
 	     cell->sky+cell->bias, cell->rms, BURN_EXP);
 
 /* Subtract out the fit */
-      if(!cell->persist[k].fiterr) 
-	 err = sub_fit(nx, ny, NX, buf, cell->persist+k, 1);
+     if(!cell->persist[k].fiterr) 
+	   err = sub_fit(nx, ny, NX, buf, cell->persist+k, 1);
    }
 
@@ -65,7 +65,8 @@
    for(k=0; k<cell->npersist; k++) {
       i = (cell->persist[k].ex - cell->persist[k].sx + 1) / 2;
-      ymid = (cell->persist[k].y0m + cell->persist[k].y1m +
-	      cell->persist[k].y0p + cell->persist[k].y1p + 2) / 4;
-      printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %8.5f %4d\n", 
+      /*ymid = (cell->persist[k].y0m + cell->persist[k].y1m +
+	      cell->persist[k].y0p + cell->persist[k].y1p + 2) / 4;*/
+        ymid = (cell->persist[k].sy + cell->persist[k].ey + 1) / 2;
+      printf("%3d %5d %3d %5d %5d %3d %5d %3d %5d %2d %8.5f %8.5f\n", 
 	     k, cell->persist[k].cx, cell->persist[k].cy,
 	     cell->persist[k].max, 
@@ -73,5 +74,5 @@
 	     cell->persist[k].ex, cell->persist[k].ey, ymid,
 	     cell->persist[k].func, 
-	     cell->persist[k].slope, NINT((cell->persist[k]).zero[i]));
+	     cell->persist[k].slope, (cell->persist[k]).zero[i]);
    }
    return(0);
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c	(revision 42392)
@@ -38,5 +38,5 @@
    while(fgets(line, 1024, fp) != NULL) {
       if(line[0] == '#') continue;
-      sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d %d",
+      sscanf(line, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %lf %d %d %d %d %d",
 	     &boxbuf[nbox].cell, &boxbuf[nbox].time,
 	     &boxbuf[nbox].cx, &boxbuf[nbox].cy,
@@ -51,15 +51,15 @@
 	     &boxbuf[nbox].slope, &boxbuf[nbox].nfit,
 	     &boxbuf[nbox].sxfit, &boxbuf[nbox].exfit,
-	     &boxbuf[nbox].fiterr);
+	     &boxbuf[nbox].fiterr, &boxbuf[nbox].eyfit);
       if(boxbuf[nbox].nfit > 0) {
-	 boxbuf[nbox].zero = (double *)calloc(boxbuf[nbox].nfit, sizeof(double));
-	 boxbuf[nbox].xfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
-	 boxbuf[nbox].yfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
-	 if(boxbuf[nbox].zero == NULL ||
-	    boxbuf[nbox].xfit == NULL ||
-	    boxbuf[nbox].yfit == NULL) {
-	    fprintf(stderr, "\rerror: failed to alloc boxbuf\n");
-	    exit(-673);
-	 }
+	     boxbuf[nbox].zero = (double *)calloc(boxbuf[nbox].nfit, sizeof(double));
+	     boxbuf[nbox].xfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
+	     boxbuf[nbox].yfit = (int *)calloc(boxbuf[nbox].nfit, sizeof(int));
+	     if(boxbuf[nbox].zero == NULL ||
+	      boxbuf[nbox].xfit == NULL ||
+	      boxbuf[nbox].yfit == NULL) {
+	     fprintf(stderr, "\rerror: failed to alloc boxbuf\n");
+	     exit(-673);
+	   }
 
 	 for(i=0; i<boxbuf[nbox].nfit; i++) {
@@ -123,5 +123,5 @@
 {
    int i, j, k, n, nbox, xs, xe;
-   int jp, kp;
+   int jp, kp,xdiff;
    int lapmax;
    OBJBOX *box;
@@ -137,65 +137,70 @@
    for(i=0; i<myMAXSIZE; i++) boxid[i] = -1;
    for(k=0; k<nbox; k++) {
-      if(box[k].exfit >= myMAXSIZE-1) continue;
-      for(i=box[k].sxfit; i<=box[k].exfit; i++) xusage[i] += 1;
+      if(box[k].ex >= myMAXSIZE-1) continue;
+      for(i=box[k].sx; i<=box[k].ex; i++) xusage[i] += 1;
    }
 
 /* Identify clusters */
    for(xs=0; xs<myMAXSIZE-1; xs++) {
-      if(xusage[xs] == 0) continue;
+    if(xusage[xs] == 0) continue;
       for(xe=xs, lapmax=0; xe<myMAXSIZE-1; xe++) {
-	 if(xusage[xe+1] == 0) break;
-	 lapmax = MAX(lapmax, xusage[xe]);
+	     if(xusage[xe+1] == 0) break;
+	     lapmax = MAX(lapmax, xusage[xe]);
       }
       if(lapmax == 1) {	/* No overlap?  No problem. */
-	 xs = xe + 1;	/* Hop to next gap */
-	 continue;
-      }
-
-/* Which boxes overlap this cluster? */
+	    xs = xe + 1;	/* Hop to next gap */
+	    continue;
+      }
+
+      /* Which boxes overlap this cluster? */
       for(k=n=0; k<nbox; k++) {
-	 if(box[k].sxfit > xe || box[k].exfit < xs) continue;
-	 boxid[n] = k;
-	 yctr[n] = (box[k].y0m+box[k].y1m+box[k].y0p+box[k].y0p+2) / 4;
-	 n++;
-      }
-
-/* Case 1: different y start => sever */
+        if(box[k].sx > xe || box[k].ex < xs) continue;
+	     boxid[n] = k;
+	     //yctr[n] = (box[k].y0m+box[k].y1m+box[k].y0p+box[k].y1p+2) / 4;
+	     yctr[n] = (box[k].sy+box[k].ey+1) / 2;
+	     n++;
+      }
+
       for(kp=0; kp<n; kp++) {
-	 k = boxid[kp];
-	 zk = 0.0;
-	 if(box[k].nfit > 0) zk = box[k].zero[box[k].nfit/2];
-	 for(jp=kp+1; jp<n; jp++) {
-	    j = boxid[jp];
-	    zj = 0.0;
-	    if(box[j].nfit > 0) zj = box[j].zero[box[j].nfit/2];
-	    if(ABS(yctr[jp]-yctr[kp]) > DIFFERENT_STREAK) {
-/* Trim back the feebler streak */
-	       if(zk > zj) {
-		  if(box[j].sxfit >= box[k].sxfit) 
-		     box[j].sxfit = MAX(box[j].sxfit, box[k].exfit+1);
-		  if(box[j].exfit <= box[k].exfit) 
-		     box[j].exfit = MIN(box[j].exfit, box[k].sxfit-1);
-	       } else {
-		  if(box[k].sxfit >= box[j].sxfit) 
-		     box[k].sxfit = MAX(box[k].sxfit, box[j].exfit+1);
-		  if(box[k].exfit <= box[j].exfit) 
-		     box[k].exfit = MIN(box[k].exfit, box[j].sxfit-1);
-	       }
-	    }
-	 }
-      }
-
-/* Case 2: pretty much the same y start => union */
-      for(kp=0; kp<n; kp++) {
-	 k = boxid[kp];
-	 for(jp=kp+1; jp<n; jp++) {
-	    j = boxid[jp];
-	    /* CZW: Since I added the initialization statement above, */
-	    /*      any box with boxid == -1 hasn't been set, and */
-	    /*      therefore needs to be dropped, I think. */
-	    if ((k < 0)||(j < 0)) { continue; }
-	    if(ABS(yctr[jp]-yctr[kp]) <= DIFFERENT_STREAK) {
-/* Merge k into j, trash k from further consideration */
+	    k = boxid[kp];
+	    zk = 0.0;
+	    if(box[k].nfit > 0) zk = box[k].zero[box[k].nfit/2];
+	    for(jp=kp+1; jp<n; jp++) {
+	     j = boxid[jp];
+	     zj = 0.0;
+	     if(box[j].nfit > 0) zj = box[j].zero[box[j].nfit/2];
+	     if(box[j].sx > box[k].sx) xdiff = box[j].sx - box[k].ex;
+	     if(box[j].sx <= box[k].sx) xdiff = box[k].sx - box[j].ex;
+
+	     /* CZW: Since I added the initialization statement above, */
+	     /*      any box with boxid == -1 hasn't been set, and */
+	     /*      therefore needs to be dropped, I think. */
+	     if ((k < 0)||(j < 0)) { continue; }
+
+        /*Tdb20220222: In the case of overlapping x-coords but not overlapping y-coords*/
+        /*I fail to see the need to change the x-range of any streaks. That just spells disaster*/
+        // Case 1: different y start => sever 
+	     if(ABS(yctr[jp]-yctr[kp]) > DIFFERENT_STREAK) {
+         // Trim back the feebler streak 
+	      if(zk > zj) {
+		    /*
+		    if(box[j].sxfit >= box[k].sxfit) 
+	       box[j].sxfit = MAX(box[j].sxfit, box[k].exfit+1);
+		    if(box[j].exfit <= box[k].exfit) 
+		    box[j].exfit = MIN(box[j].exfit, box[k].sxfit-1);
+		    */
+	      } else {
+	       /*	
+          if(box[k].sxfit >= box[j].sxfit) 
+		    box[k].sxfit = MAX(box[k].sxfit, box[j].exfit+1);
+		    if(box[k].exfit <= box[j].exfit) 
+		    box[k].exfit = MIN(box[k].exfit, box[j].sxfit-1);
+		    */
+	      }
+	     }
+
+	     /* Case 2: pretty much the same y start => union */
+	     if((ABS(yctr[jp]-yctr[kp]) <= DIFFERENT_STREAK) && (xdiff < 0) ) {
+          /* Merge k into j, trash k from further consideration */
 	       box[j].time = MAX(box[j].time, box[k].time);
 	       box[j].sx = MIN(box[j].sx, box[k].sx);
@@ -213,6 +218,7 @@
 	       box[j].sxfit = MIN(box[j].sxfit, box[k].sxfit);
 	       box[j].exfit = MAX(box[j].exfit, box[k].exfit);
-
-	       box[k].exfit = box[k].sxfit - 1;
+	       box[j].eyfit = MIN(box[j].eyfit, box[k].eyfit);
+
+	       box[k].exfit = -999;
 	       yctr[kp] = -2 * DIFFERENT_STREAK;
 	    }
@@ -226,8 +232,8 @@
 /* Excise the boxes which have been eradicated (sxfit > exfit) */
    for(k=0; k<nbox; k++) {
-      if(box[k].sxfit > box[k].exfit) {
-	 for(j=k; j<nbox-1; j++) memcpy(box+j, box+j+1, sizeof(OBJBOX));
-	 nbox--;
-	 k--;
+      if(box[k].exfit < 0) {
+	    for(j=k; j<nbox-1; j++) memcpy(box+j, box+j+1, sizeof(OBJBOX));
+	    nbox--;
+	    k--;
       }
    }
@@ -256,5 +262,5 @@
    fprintf(fp,  "# Cell: %d  sky= %d   rms= %d   bias= %d\n", 
 	   cell[0].cell, cell[0].sky, cell[0].rms, cell[0].bias);
-   fprintf(fp, "#Cell time    cx  cy  max   y0   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p F up    slope nfit sxf exf\n");
+   fprintf(fp, "#Cell time    cx  cy  max   y0   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p F up    slope nfit sxf exf fiterr eyf\n");
 
    for(j=0; j<MAXCELL; j++) {
@@ -276,5 +282,6 @@
 	    if(cell[j].persist[k].nfit <= 0) continue;
 	 }
-	 fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d\n",
+
+	 fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d %3d\n",
 		 j, cell[j].persist[k].time, 
 		 cell[j].persist[k].cx, cell[j].persist[k].cy, 
@@ -289,5 +296,5 @@
 		 cell[j].persist[k].slope, cell[j].persist[k].nfit,
 		 cell[j].persist[k].sxfit, cell[j].persist[k].exfit, 
-		 cell[j].persist[k].fiterr);
+		 cell[j].persist[k].fiterr, cell[j].persist[k].eyfit);
 	 for(i=0; i<cell[j].persist[k].nfit; i++) {
 	    fprintf(fp, "%3d %3d %8.4f\n", cell[j].persist[k].xfit[i], 
@@ -298,20 +305,20 @@
 /* Second: new burns */
       for(k=0; k<cell[j].nburn; k++) {
-	 if(!cell[j].burn[k].burned) continue;
-	 if(PERSIST_RETAIN) {
-/* Keep fits which have a dubious slope */
-	    if(cell[j].burn[k].fiterr != FIT_SLOPE_ERROR) {
+	    if(!cell[j].burn[k].burned) continue;
+	    if(PERSIST_RETAIN) {
+        /* Keep fits which have a dubious slope */
+	     if(cell[j].burn[k].fiterr != FIT_SLOPE_ERROR) {
 	       if(cell[j].burn[k].fiterr && 
-		  cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
+		    cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
 	       if(cell[j].burn[k].nfit <= 0) continue;
+	     }
+	    } else {
+	     if(cell[j].burn[k].fiterr && 
+	       cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
+	     if(cell[j].burn[k].nfit <= 0) continue;
 	    }
-	 } else {
-	    if(cell[j].burn[k].fiterr && 
-	       cell[j].burn[k].fiterr != FIT_TOP_ERROR) continue;
-	    if(cell[j].burn[k].nfit <= 0) continue;
-	 }
 
 	 i = (cell[j].burn[k].ex - cell[j].burn[k].sx + 1) / 2;
-	 fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d\n", 
+	 fprintf(fp, "%3d %7d  %3d %3d %5d %3d  %3d %3d  %3d %3d  %3d %3d %3d %3d %3d %3d %3d %3d  %1d %1d %9.6f %3d %3d %3d %d %3d\n", 
 		 j, cell[j].burn[k].time, 
 		 cell[j].burn[k].cx, cell[j].burn[k].cy,
@@ -326,5 +333,5 @@
 		 cell[j].burn[k].slope, cell[j].burn[k].nfit,
 		 cell[j].burn[k].sxfit, cell[j].burn[k].exfit, 
-		 cell[j].burn[k].fiterr);
+		 cell[j].burn[k].fiterr, cell[j].burn[k].eyfit);
 	 for(i=0; i<cell[j].burn[k].nfit; i++) {
 	    fprintf(fp, "%3d %3d %8.4f\n", cell[j].burn[k].xfit[i], 
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/stardetect.c	(revision 42392)
@@ -34,115 +34,118 @@
    starcut = STAR_THRESH/2 + cell->sky + cell->bias;
 
+   /*try a new value for thresh_lo which is 1 sigmas above the background, but with a safety limit of 1000 counts*/
+	trailthresh = cell->sky + (10*cell->rms);
+	if(trailthresh <= 1000) trailthresh = 1000;
+	trailthresh = trailthresh + cell->bias;
+
 /* Look at all the pixels which pass the burn threshold */
    for(j=0; j<ny; j++) {
       xon = -1;
       for(i=0; i<nx; i++) {
-/* Big enough?  Initialize a box and push it outwards to the trail level. */
-	 if(mask[i+j*NX] > MASK_NONE) continue;
-	 if(xon < 0) {
-	    if(data[i+j*NX] > burnthresh) {
-	       xon = i;
-	       thresh_hi = burnthresh;
-	       thresh_lo = trailthresh;
-	       if(VERBOSE & VERB_DETECT) {
-		  printf("Starting a burn at %d %d %d\n", i, j, data[i+j*NX]);
-	       }
-	    } else if(data[i+j*NX] > starthresh) {
-/* Local maximum? */
-	       if(!local_max(i, j, STAR_RADIUS, nx, ny, NX, data)) continue;
-	       if(data[i+j*NX] < veto[i+j*NX]) continue;
-	       xon = i;
-//	       thresh_hi = starthresh;
-//	       thresh_lo = starcut;
-	       thresh_hi = data[i+j*NX];
-	       thresh_lo = STAR_FRAC*(data[i+j*NX]-cell->sky-cell->bias) + 
-		  cell->sky + cell->bias;
-	       if(VERBOSE & VERB_DETECT) {
-		  printf("Starting a star at %d %d %d\n", i, j, data[i+j*NX]);
-	       }
-	    }
-	 }
-
-	 if(xon >= 0 && (data[i+j*NX] < thresh_hi || i == nx-1)) {
-
-//	 if(xon < 0 && data[i+j*NX] > burnthresh && mask[i+j*NX] == 0) xon = i;
-//	 if(xon > 0 && (data[i+j*NX] < burnthresh || i == nx-1)) {
-
-	    xoff = i<nx-1 ? i-1 : nx-1;
-	    if(nbox >= MAXBURN) {
-	       fprintf(stderr, "error: too many burn boxes\n");
-	       return(-1);
-	    }
-	    boxbuf[nbox].sx = xon;
-	    boxbuf[nbox].ex = xoff;
-	    boxbuf[nbox].sy = j;
-	    for(k=j; k<ny && data[(xon+xoff)/2+k*NX] > thresh_hi; k++);
-	    boxbuf[nbox].ey = k-1;
-	    grow_box(nx, ny, NX, data, thresh_lo, boxbuf+nbox);
-/* Fill in max and center info */
-	    boxbuf[nbox].max = 0;
-	    for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
-	       for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
-		  if(data[k+l*NX] > boxbuf[nbox].max) {
-		     boxbuf[nbox].cx = k;
-		     boxbuf[nbox].cy = l;
-		     boxbuf[nbox].max = data[k+l*NX];
-		  }
-	       }
-	    }
+         /* Big enough?  Initialize a box and push it outto trail level. */
+	      if(mask[i+j*NX] > MASK_NONE) continue;
+	      if(xon < 0) {
+	        if(data[i+j*NX] > burnthresh) {
+	          xon = i;
+	          thresh_hi = burnthresh;
+	          thresh_lo = trailthresh;
+	          if(VERBOSE & VERB_DETECT) {
+		         printf("Starting a burn at %d %d %d %d\n", i, j, data[i+j*NX],thresh_lo);
+	          }
+	        } else if(data[i+j*NX] > starthresh) {
+             /* Local maximum? */
+	          if(!local_max(i, j, STAR_RADIUS, nx, ny, NX, data)) continue;
+	          if(data[i+j*NX] < veto[i+j*NX]) continue;
+	          xon = i;
+//	          thresh_hi = starthresh;
+//	          thresh_lo = starcut;
+	          thresh_hi = data[i+j*NX];
+	          thresh_lo = STAR_FRAC*(data[i+j*NX]-cell->sky-cell->bias) + 
+		       cell->sky + cell->bias;
+	          if(VERBOSE & VERB_DETECT) {
+		         printf("Starting a star at %d %d %d %d\n", i, j, data[i+j*NX],thresh_lo);
+	          }
+	        }
+	      }
+
+	      if(xon >= 0 && (data[i+j*NX] < thresh_hi || i == nx-1)) {
+	        xoff = i<nx-1 ? i-1 : nx-1;
+	        if(nbox >= MAXBURN) {
+	           fprintf(stderr, "error: too many burn boxes\n");
+	           return(-1);
+	        }
+	        boxbuf[nbox].sx = xon;
+	        boxbuf[nbox].ex = xoff;
+	        boxbuf[nbox].sy = j;
+	        for(k=j; k<ny && data[(xon+xoff)/2+k*NX] > thresh_hi; k++);
+	        boxbuf[nbox].ey = k-1;
+	        grow_box(nx, ny, NX, data, thresh_lo, boxbuf+nbox);
+           /* Fill in max and center info */
+	        boxbuf[nbox].max = 0;
+	        for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
+	         for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
+		       if(data[k+l*NX] > boxbuf[nbox].max) {
+		        boxbuf[nbox].cx = k;
+		        boxbuf[nbox].cy = l;
+		        boxbuf[nbox].max = data[k+l*NX];
+		       }
+	         }
+	        }
+
 /* A box which triggered on starthresh cannot envelop burnthresh */
-	    if(thresh_hi < burnthresh && boxbuf[nbox].max > burnthresh) {
-	       if(VERBOSE & VERB_DETECT) {
-		  printf("Ditching box %d %d %d %d  max = %d > burnthresh\n",
-			 boxbuf[nbox].sx, boxbuf[nbox].sy,
-			 boxbuf[nbox].ex, boxbuf[nbox].ey,
-			 boxbuf[nbox].max);
-		  fflush(stdout);
-	       }
+	        if(thresh_hi < burnthresh && boxbuf[nbox].max > burnthresh) {
+	          if(VERBOSE & VERB_DETECT) {
+		        printf("Ditching box %d %d %d %d  max = %d > burnthresh\n",
+			     boxbuf[nbox].sx, boxbuf[nbox].sy,
+			     boxbuf[nbox].ex, boxbuf[nbox].ey,
+			     boxbuf[nbox].max);
+		        fflush(stdout);
+	          }
 /* Mark the veto mask so as not to trigger on this one again */
-	       for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
-		  for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
-		     veto[k+l*NX] = boxbuf[nbox].max;
-		  }
-	       }
-	       xon = -1;
-	       continue;
-	    }
+	          for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
+		         for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
+		          veto[k+l*NX] = boxbuf[nbox].max;
+		         }
+	          }
+	          xon = -1;
+	          continue;
+	        }
 
 /* Mask the pixels so as not to catch them again! */
-	    for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
-	       for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
-		  mask[k+l*NX] = MASK_CTR;
-	       }
-	    }
-	    boxbuf[nbox].sat = boxbuf[nbox].max >= burnthresh;
-	    boxbuf[nbox].max -= cell->sky + cell->bias;
+	        for(l=boxbuf[nbox].sy; l<=boxbuf[nbox].ey; l++) {
+	         for(k=boxbuf[nbox].sx; k<=boxbuf[nbox].ex; k++) {
+		       mask[k+l*NX] = MASK_CTR;
+	         }
+	        }
+	        boxbuf[nbox].sat = boxbuf[nbox].max >= burnthresh;
+	        boxbuf[nbox].max -= cell->sky + cell->bias;
 /* Burned only if it triggered on burnthresh, not if it enveloped it */
-	    boxbuf[nbox].burned = thresh_hi >= burnthresh;
+	        boxbuf[nbox].burned = thresh_hi >= burnthresh;
 /* But take a hard look at stars with really bright centers... */
-	    if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) {
-	       burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox);
-	    }
-	    boxbuf[nbox].time = cell->time;
-	    boxbuf[nbox].cell = cellnum;
-	    boxbuf[nbox].stamp = NULL;
-	    boxbuf[nbox].func = FUNC_NONE;
-	    boxbuf[nbox].sxfit = boxbuf[nbox].sx;
-	    boxbuf[nbox].exfit = boxbuf[nbox].ex;
-	    boxbuf[nbox].nfit = 0;
-	    boxbuf[nbox].zero = NULL;
-	    boxbuf[nbox].xfit = boxbuf[nbox].yfit = NULL;
-	    if(VERBOSE & VERB_DETECT) {
-	       printf("New box at %d %d %d %d", 
+	        if(!boxbuf[nbox].burned && boxbuf[nbox].max > MAX_THRESH) {
+	         burn_test(nx, ny, NX, data, cell->rms, mask, boxbuf+nbox);
+	        }
+	        boxbuf[nbox].time = cell->time;
+	        boxbuf[nbox].cell = cellnum;
+	        boxbuf[nbox].stamp = NULL;
+	        boxbuf[nbox].func = FUNC_NONE;
+ 
+	        boxbuf[nbox].sxfit = boxbuf[nbox].sx;
+	        boxbuf[nbox].exfit = boxbuf[nbox].ex;
+	        boxbuf[nbox].nfit = 0;
+	        boxbuf[nbox].zero = NULL;
+	        boxbuf[nbox].xfit = boxbuf[nbox].yfit = NULL;
+	        if(VERBOSE & VERB_DETECT) {
+	         printf("New box at %d %d %d %d", 
 		      boxbuf[nbox].sx, boxbuf[nbox].sy,boxbuf[nbox].ex, boxbuf[nbox].ey);
-	       printf("   %d %d %d %d", 
+	         printf("   %d %d %d %d", 
 		      boxbuf[nbox].x0m, boxbuf[nbox].x0p,boxbuf[nbox].x1m, boxbuf[nbox].x1p);
-	       printf("   %d %d %d %d\n", 
+	         printf("   %d %d %d %d\n", 
 		      boxbuf[nbox].y0m, boxbuf[nbox].y0p,boxbuf[nbox].y1m, boxbuf[nbox].y1p);
-	       fflush(stdout);
-	    }
-	    xon = -1;
-	    nbox++;
-	 }
+	         fflush(stdout);
+	        }
+	        xon = -1;
+	        nbox++;
+	      }
       }
    }
Index: /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c
===================================================================
--- /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c	(revision 42391)
+++ /trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c	(revision 42392)
@@ -14,5 +14,5 @@
 #include "burnparams.h"
 
-STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE];
+STATIC double ybuf[MAXSIZE], zbuf[MAXSIZE], wbuf[MAXSIZE], wbuf0[MAXSIZE];
 
 /****************************************************************/
@@ -86,5 +86,5 @@
 #define MAXWGT 2.0	/* Maximum allowed weight */
 #define MINWGT 0.1	/* Minimum allowed weight */
-#define MAXDEV 0.5	/* Maximum ln deviation above first pass fit */
+#define MAXDEV 1.	/* Maximum ln deviation above first pass fit */
 #define MAXRMS 4.0	/* Minimum factor of RMS to be chopped */
 
@@ -94,7 +94,7 @@
 		    int up, int bckgnd, int rms, int fitfunc)
 {
-   int i, j, k, err, nfit;
+   int i, j, k, err, nfit,firstfitpix;
    int xs, xe, y0, y1, y2, dy, yfit;
-   double slope, zero, zsum, wsum, trial;
+   double slope, zero, zsum, wsum, trial,zfun;
    char fooname[60];
    FILE *fp=stdout;
@@ -105,23 +105,27 @@
    }
 
-   xs = box->sxfit;
-   xe = box->exfit;
-
+   /*xs = box->sxfit; */
+   /*xe = box->exfit; */
+   /*TdB20220222: Do not use s/exfit here, which will fail in the case of positive slopes*/
+   xs = box->sx;
+   xe = box->ex;
+   
    if(up) {
-      y0 = (box->y0m + box->y1m + box->y0p + box->y0p + 2) / 4;
-/* Maybe a bit better?  Maybe not... */
-//      y0 = box->ey;
-      y1 = box->ey + FIT_EDGE;	/* Start of fit, budged away from the burn */
-      y2 = ny - 1;		/* End of fit */
+      //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
+      y0 = (box->sy + box->ey + 1) / 2;
+      y1 = box->ey + FIT_EDGE;   /* Start of fit, budged away from the burn */
+      y2 = ny - 1;      /* End of fit */
       dy = +1;
-      yfit = box->sy;		/* Start of validity, TBD later by col */
+      yfit = box->sy;      /* Start of validity, TBD later by col */
+      firstfitpix = MAX((box->ey - box->sy)*2,20) ;  /* Number of pixels to include in the initial fit*/
    } else {
-      y0 = (box->y0m + box->y1m + box->y0p + box->y0p + 2) / 4;
-      y1 = 0;			/* Start of fit */
-      y2 = box->sy - FIT_EDGE;	/* End of fit, budged away from the burn */
-//      y2 = box->ey - FIT_EDGE;	/* End of fit, budged away from the burn */
-//      y2 = y0;	/* End of fit, budged away from the burn */
+      //y0 = (box->y0m + box->y1m + box->y0p + box->y1p + 2) / 4;
+      y0 = (box->sy + box->ey + 1) / 2;
+      y1 = 0;        /* Start of fit */
+      y2 = y0 - FIT_EDGE;  /* End of fit, budged away from the burn */
       dy = -1;
-      yfit = box->ey;		/* Start of validity, TBD later by col */
+      yfit = box->ey;      /* Start of validity, TBD later by col */
+      firstfitpix = 100;  /* Number of pixels to include in the initial fit*/
+
    }
 
@@ -147,4 +151,5 @@
    box->up = up;
    box->y0 = y0;
+   box->eyfit = up?y2:y1;
    for(i=0; i<xe-xs+1; i++) {
       box->zero[i] = 0.0;
@@ -157,11 +162,11 @@
       box->slope = -1.0;
       box->fiterr = FIT_TOP_ERROR;
-/* No fit here because the burn is to top, but the persist may be fitable */
+      /* No fit here because the burn is to top, but the persist may be fitable */
       if(box->sy > FIT_EDGE) {
-	 box->nfit = xe - xs + 1;
-/* However, if burn extends to bottom, persist will not be able to fit */
+	    box->nfit = xe - xs + 1;
+       /* However, if burn extends to bottom, persist will not be able to fit */
       } else {
-	 box->nfit = 1;		/* Just enough to keep it alive */
-	 box->func = BURN_BLASTED;
+	    box->nfit = 1;		/* Just enough to keep it alive */
+	    box->func = BURN_BLASTED;
       }
       return(0);
@@ -174,4 +179,5 @@
 
 /* Accumulate points to fit */
+/* use wbuf0 to select only the first FIRSTFITPIX pixels for first fit */
    nfit = y2 - y1 + 1;
    for(j=y1; j<=y2; j++) {
@@ -179,24 +185,29 @@
       zbuf[j]=0.0;
       for(i=xs, k=0; i<=xe; i++) {
-//	 printf("%4d %4d %7d %7d\n", i, j, data[i+j*NX], mask[i+j*NX]);
-	 if(data[i+j*NX] != NODATA && mask[i+j*NX] == MASK_NONE) {
-//	 if(data[i+j*NX] != NODATA) {
-	    zbuf[j] += data[i+j*NX];
-	    k++;
-	 }
-      }
-//      printf("%4d %4d %7d\n", k, j, zbuf[j]);
+	     if(data[i+j*NX] != NODATA && (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO) ) {
+	      zbuf[j] += data[i+j*NX];
+	      k++;
+	     }
+        if(mask[i+j*NX] != MASK_NONE) {
+        }
+      }
       zbuf[j] = zbuf[j] / MAX(1,k) - bckgnd;
-/* Linearize for fit... */
+
+/*    Linearize for fit... */
       if(fitfunc == BURN_PWR) ybuf[j] = log(ybuf[j]);
       wbuf[j] = 0;
+      wbuf0[j] = 0;
       if(zbuf[j] > 1) {
-	 wbuf[j] = MIN(MAXWGT, (zbuf[j] / rms));  /* Quadratic not good */
-	 wbuf[j] = MAX(MINWGT, wbuf[j]);
-	 zbuf[j] = log(zbuf[j]);
+	     wbuf[j] = MIN(MAXWGT, (zbuf[j] / rms));  /* Quadratic not good */
+	     wbuf[j] = MAX(MINWGT, wbuf[j]);
+	     zbuf[j] = log(zbuf[j]);
+        wbuf0[j] = wbuf[j];
+        if(j < (y2-firstfitpix) && !up) wbuf0[j] = 0;
+        if(j > (y1+firstfitpix) && up) wbuf0[j] = 0;
+        //printf("fit: %4d %.3f %.3f %.3f %.3f\n", j, ybuf[j], zbuf[j], wbuf[j], wbuf0[j]);
       }
    }
 /* First pass fit */
-   err = wlinearfit(nfit, ybuf+y1, zbuf+y1, wbuf+y1, &slope, &zero);
+   err = wlinearfit(nfit, ybuf+y1, zbuf+y1, wbuf0+y1, &slope, &zero);
    if(VERBOSE & VERB_FIT) {
       printf("nfit= %d  err= %d  slope= %.3f  zero= %.3f\n", 
@@ -210,6 +221,7 @@
 /* Trim away any really big positive deviations */
    for(j=y1; j<=y2; j++) {
-      if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV && 
-	 zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
+//      if(zbuf[j]-(ybuf[j]*slope+zero) > MAXDEV && 
+//	 zbuf[j] > log(MAXRMS*rms)) wbuf[j] = 0;
+      if(abs(zbuf[j]-(ybuf[j]*slope+zero)) > MAXDEV) wbuf[j] = 0;
       if(VERBOSE & VERB_FITPROF) {
 	 fprintf(fp, "%9.4f %9.4f %9.4f %3d\n", ybuf[j], zbuf[j], wbuf[j], j);
@@ -239,4 +251,6 @@
 //      if(trial > 2*rms) box->nfit = -box->nfit;
       if(trial > 2*rms) box->func = BURN_POSSLOPE;
+      box->sxfit = xs;
+      box->exfit = xs - 1;
       return(-1);
    }
@@ -244,16 +258,23 @@
 /* Reconstruct the fit for comparison with the data col by col */
    for(j=y1; j<=y2; j++) zbuf[j] = exp(ybuf[j]*slope);
-/* And tag on the extras to determine the end of the fit */
+   /* And tag on the extras to determine the end of the fit */
    for(j=yfit; dy*j<dy*(up?y1:y2); j+=dy) {
       if(fitfunc == BURN_EXP) {
-//	 zbuf[j] = exp(Y_SCALE*MAX(0,dy*(j-y0))*slope);
-	 zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
+  	     zbuf[j] = exp(Y_SCALE*dy*(j-y0)*slope);
       } else {
-	 if(dy*(j-y0) <= 0) zbuf[j] = exp(log(Y_SCALE*1.0)*slope);
-	 else zbuf[j] = exp(log(Y_SCALE*dy*(j-y0))*slope);
-      }
-   }
-
-/* Save the fit information */
+	     if(dy*(j-y0) <= 0) zbuf[j] = exp(log(Y_SCALE*1.0)*slope);
+	     else zbuf[j] = exp(log(Y_SCALE*dy*(j-y0))*slope);
+      }
+   }
+
+/* find where the fitted function becomes less than half the rms and cut the trail there*/
+   for(j=y1; j<=y2; j++) {
+      zfun = exp((ybuf[j]*slope)+zero);
+      if(zfun > NEGLIGIBLE_TRAIL*rms) break;
+   }   
+   //For burns, just let it go to the top
+   if(j < y2 && !up) box->eyfit = j;
+
+/* Determine good x-ranges and Save the fit information */
    box->slope = slope;
    box->nfit = 0;
@@ -262,37 +283,36 @@
       zsum = wsum = 0.0;
       for(j=y1; j<=y2; j++) {
-	 if(wbuf[j] > 0 && data[i+j*NX] != NODATA && 
-	    mask[i+j*NX] == MASK_NONE) {
-//	 if(wbuf[j] > 0 && data[i+j*NX] != NODATA) {
-//	    wsum += wbuf[j];
-//	    zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j];
-	    wsum += zbuf[j];
-	    zsum += data[i+j*NX] - bckgnd;
-	 }
+	     if(wbuf[j] > 0 && data[i+j*NX] != NODATA && 
+	     (mask[i+j*NX] == MASK_NONE || mask[i+j*NX] == MASK_SAT_HALO)) {
+//	       if(wbuf[j] > 0 && data[i+j*NX] != NODATA) {
+//	       wsum += wbuf[j];
+//	       zsum += wbuf[j] * (data[i+j*NX] - bckgnd) / zbuf[j];
+	       wsum += zbuf[j];
+	       zsum += data[i+j*NX] - bckgnd;
+	     }
       }
       if(zsum < 0 || wsum <= 0) {
-	 zsum = 0.0;
+	     zsum = 0.0;
       } else {
-	 zsum = zsum / wsum;
-      }
-/* FIXME: what's a really good criterion for negligible fit? */
-/* 100 pixels up fit is zsum or ~zsum/e */
-//      if(zsum > NEGLIGIBLE_TRAIL*rms) {
+	     zsum = zsum / wsum;
+      }
+      /* FIXME: what's a really good criterion for negligible fit? */
+      /* 100 pixels up fit is zsum or ~zsum/e */
+//    if(zsum > NEGLIGIBLE_TRAIL*rms) {
+      /* 100113: but also evaluate at the end of the fit for stubby trails */
       if(zsum > NEGLIGIBLE_TRAIL*rms ||
-/* 100113: but also evaluate at the end of the fit for stubby trails */
-	 zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
-
-/* Ascertain the starting point of where the fit is good */
+	   zbuf[y2]*zsum > NEGLIGIBLE_TRAIL*rms) {
+        /* Ascertain the starting point of where the fit is good */
 //      box->yfit[box->nfit] = yfit;	/* Vanilla, misses center */
 //      box->yfit[box->nfit] = y0;	/* Center, no adjust */
-	 for(j=yfit; dy*j<=dy*(up?y2:y1); j+=dy) {
-	    trial = zsum * zbuf[j];
-/* First time data - fit is closer to sky than is data */
-	    if(data[i+j*NX] > bckgnd + 0.5*trial) break;
-	 }
-	 box->yfit[box->nfit] = j;
-	 box->zero[box->nfit] = zsum;
-	 box->xfit[box->nfit] = i;
-	 box->nfit += 1;
+	     for(j=yfit; dy*j<=dy*(up?y2:y1); j+=dy) {
+	       trial = zsum * zbuf[j];
+/*        First time data - fit is closer to sky than is data */
+	       if(data[i+j*NX] > bckgnd + 0.5*trial) break;
+	     }
+	     box->yfit[box->nfit] = j;
+	     box->zero[box->nfit] = zsum;
+	     box->xfit[box->nfit] = i;
+	     box->nfit += 1;
       }
    }
@@ -306,4 +326,5 @@
       box->fiterr = FIT_ALL_GONE;
    }
+
 //   if(box->nfit <= 0) fprintf(stderr, "WHOA, got nfit = 0\n");
    return(0);