Index: /trunk/ppShutter/src/.cvsignore
===================================================================
--- /trunk/ppShutter/src/.cvsignore	(revision 8927)
+++ /trunk/ppShutter/src/.cvsignore	(revision 8928)
@@ -2,6 +2,6 @@
 Makefile
 Makefile.in
-ppMerge
 config.h
 config.h.in
 stamp-h1
+ppShutter
Index: /trunk/ppShutter/src/ppShutter.c
===================================================================
--- /trunk/ppShutter/src/ppShutter.c	(revision 8927)
+++ /trunk/ppShutter/src/ppShutter.c	(revision 8928)
@@ -5,5 +5,5 @@
 void usage () {
 
-    fprintf (stderr, "USAGE: ppShutter (output) (input) (input) (input) ...\n");
+    fprintf (stderr, "USAGE: ppShutter (output) (outscale) (input) (input) (input) ...\n");
     fprintf (stderr, " minimum number of input images: 5\n");
     exit (1);
@@ -19,7 +19,7 @@
 
     // minimum number of input images: 5
-    if (argc < 7) usage ();
+    if (argc < 8) usage ();
 
-    int Ninput = argc - 2;
+    int Ninput = argc - 3;
 
     // load the input images (for now, load them all into memory)
@@ -28,5 +28,5 @@
     psArray *images  = psArrayAlloc (Ninput);
     for (int i = 0; i < Ninput; i++) {
-	fits = psFitsOpen (argv[i-2], "r");
+	fits = psFitsOpen (argv[i+3], "r");
 	files->data[i] = fits;
 
@@ -37,4 +37,5 @@
 	images->data[i] = image;
     }
+    images->n = files->n = headers->n = Ninput;
 
     // find the exposures times for each input images
@@ -43,7 +44,10 @@
     psVector *fluxRef = psVectorAlloc (Ninput, PS_TYPE_F32);
     psVector *normCts = psVectorAlloc (Ninput, PS_TYPE_F32);
+    psVector *normErr = psVectorAlloc (Ninput, PS_TYPE_F32);
     for (int i = 0; i < Ninput; i++) {
 	exptime->data.F32[i] = psMetadataLookupF32 (NULL, headers->data[i], "EXPTIME");
     }
+    exptime->n = fluxRef->n = normCts->n = normErr->n = Ninput;
+
     psVector *index = psVectorSortIndex (NULL, exptime);
 
@@ -52,5 +56,5 @@
 
     // statistics used to measure initial parameters
-    psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+    psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
     // define the reference region
@@ -62,4 +66,5 @@
 	myStats = psImageStats(myStats, image, NULL, 0);
 	fluxRef->data.F32[i] = myStats->sampleMean;
+	fprintf (stderr, "ref flux for image %d: %f\n", i, fluxRef->data.F32[i]);
     }	
 
@@ -79,4 +84,5 @@
 	psPlane *coord = samples->data[j];
 
+	fprintf (stderr, "norm cnts for region %d: ", j);
 	for (int i = 0; i < Ninput; i++) {
 	    psImage *image = images->data[index->data.S32[i]];
@@ -84,16 +90,26 @@
 	    region = psRegionForSquare(coord->x, coord->y, 10);
 	    region = psRegionForImage(image, region);
+	    // fprintf (stderr, "%f,%f:%f,%f : ", region.x0, region.y0, region.x1, region.y1);
 
-	    myStats = psImageStats(myStats, image, NULL, 0);
+	    psImage *raster = psImageSubset (image, region);
+	    myStats = psImageStats(myStats, raster, NULL, 0);
 	    normCts->data.F32[i] = myStats->sampleMean / fluxRef->data.F32[i];
+	    normErr->data.F32[i] = myStats->sampleStdev / fluxRef->data.F32[i];
+	    fprintf (stderr, "%f ", normCts->data.F32[i]);
+	    // fprintf (stderr, "%f\n", myStats->sampleMean);
+	    
+	    psFree (raster);
 	}	
+	fprintf (stderr, "\n");
 	
 	pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, normCts);
-	pmShutterCorrPars *params = pmShutterCorrectionFullFit (exptime, normCts, guess);
+	pmShutterCorrPars *params = pmShutterCorrectionFullFit (exptime, normCts, normErr, guess);
 
 	refOffsetSample->data.F32[j] = params->offref;
     }
+    refOffsetSample->n = Nsamples;
     myStats = psVectorStats (myStats, refOffsetSample, NULL, NULL, 0);
     float refOffset = myStats->sampleMean;
+    fprintf (stderr, "refOffset : %f +/- %f\n", myStats->sampleMean, myStats->sampleStdev);
 
     // now run through all pixels in the images and measure the shutter parameters for each pixel
@@ -101,11 +117,14 @@
     psImage *offset = psImageAlloc (refimage->numCols, refimage->numRows, PS_TYPE_F32);
 
+    int dotInterval = refimage->numRows / 32;
     for (int ny = 0; ny < refimage->numRows; ny++) {
+	if (ny % dotInterval == 0) fprintf (stderr, ".");	
 	for (int nx = 0; nx < refimage->numCols; nx++) {
 	    for (int i = 0; i < Ninput; i++) {
 		psImage *image = images->data[index->data.S32[i]];
 		normCts->data.F32[i] = image->data.F32[ny][nx] / fluxRef->data.F32[i];
+		normErr->data.F32[i] = sqrt(image->data.F32[ny][nx]) / fluxRef->data.F32[i];
 	    }
-	    pmShutterCorrPars *params = pmShutterCorrectionLinFit (exptime, normCts, refOffset);
+	    pmShutterCorrPars *params = pmShutterCorrectionLinFit (exptime, normCts, normErr, refOffset);
 	    pattern->data.F32[ny][nx] = params->scale;
 	    offset->data.F32[ny][nx] = params->offset;
Index: /trunk/ppShutter/test/tests.pro
===================================================================
--- /trunk/ppShutter/test/tests.pro	(revision 8928)
+++ /trunk/ppShutter/test/tests.pro	(revision 8928)
@@ -0,0 +1,71 @@
+
+# make a simple shutter offset map : off = Ax((x-Nx/2)/Nx)^2 + Ay((y-Ny/2)/Ny)^2
+macro mkshutter.model
+ if ($0 != 6)
+   echo "USAGE: mkshutter.model (Nx) (Ny) (Ax) (Ay) (Ao)"
+   break
+ end
+
+ $Nx = $1
+ $Ny = $2
+
+ mcreate tmpb $Nx $Ny
+ set tmpx = xramp(tmpb)
+ set tmpy = yramp(tmpb)
+
+ set tmpX = (2*tmpx/$Nx) - 1
+ set tmpY = (2*tmpy/$Ny) - 1
+
+ set model = $3*tmpX^2 + $4*tmpY^2 + $5
+ delete tmpb tmpx tmpy tmpX tmpY
+end
+
+macro mkshutter.image
+ if ($0 != 4)
+   echo "USAGE: mkshutter.image (buffer) (exptime) (flux)"
+   break
+ end
+
+ set $1 = $3*($2 + model)
+ 
+ gaussdev noisev {$Nx*$Ny} 0.0 1.0
+ dimenup noisev noisem $Nx $Ny
+ set $1 = $1 + noisem*sqrt($1)
+ keyword $1 EXPTIME -wf $2
+end
+
+macro mkshutter.testdata
+  mkshutter.model 800 800 0.1 0.2 0.0
+
+  tv model 0 0.2
+  resize 800 800
+  center 400 400
+
+  mkshutter.image image 0.1 10000
+  wd image im.0.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 0.2 10000
+  wd image im.1.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 0.3 10000
+  wd image im.2.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 0.5 10000
+  wd image im.3.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 1.0 10000
+  wd image im.4.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 2.0 10000
+  wd image im.5.fits -bitpix -32 
+  stv image 0.2
+
+  mkshutter.image image 5.0 10000
+  wd image im.6.fits -bitpix -32 
+  stv image 0.2
+end