Index: trunk/ppShutter/src/ppShutter.c
===================================================================
--- trunk/ppShutter/src/ppShutter.c	(revision 8901)
+++ 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;