Index: trunk/stac/scripts/testOffsets.pl
===================================================================
--- trunk/stac/scripts/testOffsets.pl	(revision 3680)
+++ trunk/stac/scripts/testOffsets.pl	(revision 3680)
@@ -0,0 +1,54 @@
+#!/usr/bin/perl
+#
+# testOffsets.pl
+# Concerned that different offsets are producing different sensitivity to finding CRs.
+# For example, half pixel offsets in each direction seem to produce lots of false positives,
+# but tiny offsets produce few CRs found.  Need statistics.
+
+my @offsets = ();
+my @numcr = ();
+my $badNum = 0;
+for (my $iter = 0; $iter < 500; $iter++) {
+    system "./testParams.pl";
+    my $bad = 0;
+    for (my $i = 0; $i < 4; $i++) {
+
+	# Number of CRs found
+	my $wc = `wc -l test_$i.fits.cr`;
+	my ($number) = ($wc =~ /^\s*(\d+)/);
+	print "Found $number CRs in image $i\n";
+	if ($number > 2000) {
+	    $bad = 1;
+	}
+
+	# Offset
+	my $mapFile;
+	open $mapFile, "test_$i.fits.map";
+	my @map = <$mapFile>;
+	print @map . "\n";
+	close $mapFile;
+	my ($x) = ($map[1] =~ /^(\S+)/);
+	my ($y) = ($map[2] =~ /^(\S+)/);
+	my $offset = sqrt($x**2 + $y**2);
+	print "Offset for image $i is $offset\n";
+
+	push @offsets, $offset;
+	push @numcr, $number;
+    }
+
+    if ($bad) {
+	$badNum++;
+	system "mkdir bad_$badNum";
+	system "cp -f test_* testout.fits* bad_$badNum";
+    }
+}
+
+my $cr;
+open $cr, "> cr.dat";
+for (my $i = 0; $i <= $#offsets; $i++) {
+    print $cr "$offsets[$i] $numcr[$i]\n";
+}
+close $cr;
+system "sort -k 1 -n cr.dat > cr.dat.sort";
+
+__END__
Index: trunk/stac/scripts/testParams.pl
===================================================================
--- trunk/stac/scripts/testParams.pl	(revision 3680)
+++ trunk/stac/scripts/testParams.pl	(revision 3680)
@@ -0,0 +1,176 @@
+#!/usr/bin/perl
+#
+# testParams.pl
+# Program to run through the parameters for STAC.
+
+
+@REJECT = (3.0); # List of "reject" parameters
+@FRAC = (0.5); # List of "frac" parameters
+@GRAD = (0.6); # List of "grad" parameters
+@SEEING = (2.5); # List of seeing
+$ITERATIONS = 1; # Number of iterations for each parameter combination
+
+$RUN = "rm -f testout.fits* chi2_*.fits test_[0-3].fits.\{err,grad,mask,rejmap,shift*}; ./stac -v testout.fits test_0.fits test_1.fits test_2.fits test_3.fits";		# Command with which to run the STAC program
+$RUNREJECT = "-k";		# Rejection parameter flag
+$RUNFRAC = "-f";		# Frac parameter flag
+$RUNGRAD = "-G";		# Grad parameter flag
+$CREATEDIR = "/home/mithrandir/price/images/"; # Directory in which the image creation routines reside
+$CREATE = "./mkfakeimages.pl -diag"; # Command for the image creation
+$CREATESEEING = "-seeing";	# Switch for seeing in the image creation
+$CREATENOCRS = "-numcrs 0";	# Switch for no CRs in the image creation
+
+$WORKINGDIR = `pwd`;		# Working directory (current directory)
+chomp($WORKINGDIR);
+
+# Iterate over the parameter space
+open OUTPUT, "> testParams.out";
+foreach $seeing (@SEEING) {
+    foreach $reject (@REJECT) {
+	foreach $frac (@FRAC) {
+	    foreach $grad (@GRAD) {
+		@trueRejected = (); # Number of true CRs rejected in each iteration
+		@falseRejected = (); # Number of false CRs rejected in each iteration
+		@dBonB = (); # Relative difference between input and output slope
+		
+		# Do each iteration
+		for ($i=1;$i<=$ITERATIONS;$i++) {
+		    
+		    ### Run the program on the data with CRs
+		    
+		    # Make the data
+		    chdir $CREATEDIR;
+		    print "$CREATE $CREATESEEING $seeing\n";
+		    system "$CREATE $CREATESEEING $seeing";
+		    system "mv -f test_[0-3]\{.fits,.fits.map,.cr\} $WORKINGDIR";
+		    chdir $WORKINGDIR;
+		    
+		    # Run the program and count the number of CRs
+		    print "$RUN $RUNREJECT $reject $RUNFRAC $frac $RUNGRAD $grad |\n";
+		    open STAC, "$RUN $RUNREJECT $reject $RUNFRAC $frac $RUNGRAD $grad |";
+		    $numCRs = 0;	# Number of CRs masked
+		    while (<STAC>) {
+			print;
+			if (/(\d+) pixels masked in image \d+/) {
+			    $numCRs += $1;
+			}
+		    }
+		    close STAC;
+		    printf "---> %d pixels masked in total.\n", $numCRs;
+
+		    my @trueFlux = (); # Flux of pixels that are real CRs
+		    my @falseFlux = ();	# Flux of pixels that were masked, but not CRs
+		    for (my $i = 0; $i < 4; $i++) {
+			# Read list of CRs injected into image
+			open CR, "test_$i.cr";
+			my %originalCRs = ();
+			while (<CR>) {
+			    my ($x,$y,$flux) = split;
+			    $originalCRs{"$x,$y"} = $flux;
+			}
+			close CR;
+
+			# Read list of CRs masked by program
+			open CR, "test_$i.fits.cr";
+			while (<CR>) {
+			    my ($x,$y,$junk,$flux,$frac,$grad) = split;
+			    # See if it was injected
+			    if (exists($originalCRs{"$x,$y"})) {
+				push @trueFlux, $flux;
+			    } else {
+				push @falseFlux, $flux;
+			    }
+			}
+		    }
+
+		    # Number of true and false CRs
+		    push @trueRejected, scalar @trueFlux;
+		    push @falseRejected, scalar @falseFlux;
+
+		    print "True: " . scalar @trueFlux . "   False: " . scalar @falseFlux . "\n";
+
+		    # Get input slope
+		    system "cat test_[0-3].cr | sort -k 3 -n | awk \'{num++; print num,\$3;}\' > test_in.cr";
+		    open FIT, "guessline.pl test_in.cr |";
+		    while (<FIT>) {
+			if (/b = (\S+)/) {
+			    $inSlope = $1;
+			    last;
+			}
+		    }
+		    close FIT;
+
+		    # Get output slope
+		    open CR, "| sort -n | awk \'{num++; print num,\$1;}\' > test_out.cr";
+		    foreach $flux (@trueFlux) {
+			print CR $flux . "\n";
+		    }
+		    close CR;
+		    
+		    open FIT, "guessline.pl test_out.cr |";
+		    while (<FIT>) {
+			if (/b = (\S+)/) {
+			    $outSlope = $1;
+			    last;
+			}
+		    }
+		    close FIT;
+		    print "inSlope: $inSlope\noutSlope: $outSlope\n";
+
+		    push @dBonB, abs($outSlope - $inSlope) / $inSlope;
+		    
+		}
+		
+		# Calculate the mean
+		$trueMean = mean(\@trueRejected);
+		$falseMean = mean(\@falseRejected);
+		$trueMedian = median(\@trueRejected);
+		$falseMedian = median(\@falseRejected);
+		$dBonBMean = mean(\@dBonB);
+		$dBonBMedian = median(\@dBonB);
+		
+		printf "$seeing\t$reject\t$frac\t$grad\t%.1f\t%.1f\t%.1f\t%.1f\t%.3f\t%.3f\n",$trueMean,$trueMedian,$falseMean,$falseMedian,$dBonBMean,$dBonBMedian;
+		printf OUTPUT "$seeing\t$reject\t$frac\t$grad\t%.1f\t%.1f\t%.1f\t%.1f\t%.3f\t%.3f\n",$trueMean,$trueMedian,$falseMean,$falseMedian,$dBonBMean,$dBonBMedian;
+	    }
+	}
+    }
+}
+close OUTPUT;
+
+
+
+
+# Return the mean
+sub mean
+{
+    my ($arrayRef) = @_;
+
+    my @array = @$arrayRef;
+    my $mean = 0;
+    foreach my $value (@array) {
+	$mean += $value;
+    }
+    return $mean / (scalar @array);
+}
+
+# Return the median
+sub median
+{
+    my ($arrayRef) = @_;
+    my @array = @$arrayRef;
+    # Insertion sort
+    my @sort = ();
+    foreach my $value (@array) {
+	for ($i=0;$i<=$#sort;$i++) {
+	    if ($sort[$i] < $value) {
+		last;
+	    }
+	}
+	splice @sort, $i, 0, $value;
+    }
+    # Median is the middle one.
+    if (scalar @sort % 2 == 1) {
+	return $sort[(scalar @sort  - 1) / 2];
+    } else {
+	return ($sort[(scalar @sort)/2] + $sort[(scalar @sort)/2 - 1])/2;
+    }
+}
Index: trunk/stac/src/combine.c
===================================================================
--- trunk/stac/src/combine.c	(revision 3680)
+++ trunk/stac/src/combine.c	(revision 3680)
@@ -0,0 +1,70 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+#include "combineConfig.h"
+
+int main(int argc, char **argv)
+{
+    // Set trace levels
+    psTraceSetLevel(".", 0);
+    psTraceSetLevel("stac.checkMemory", 10);
+    psTraceSetLevel("stac.read", 10);
+    psTraceSetLevel("stac.scales", 7);
+    psTraceSetLevel("stac.rescale", 10);
+    psTraceSetLevel("stac.combine", 10);
+    psTraceSetLevel("combine", 10);
+
+    // Set logging level
+    psLogSetLevel(9);
+
+    // Parse command line
+    combineConfig *config = combineParseConfig(argc, argv); // Configuration
+
+    psArray *images = stacReadImages(config->inNames); // The images
+    psArray *errors = stacErrorImages(images, config->gain, config->readnoise); // Error images
+
+    // Calculate scales between images
+    psVector *scales = NULL;		// Relative scales between images
+    psVector *offsets = NULL;		// Offsets between images
+    (void)stacScales(&scales, &offsets, images, config->starFile, config->starMap, 0.0, 0.0, config->aper);
+
+    // Set the saturation and bad values
+    psVector *saturated = psVectorAlloc(images->n, PS_TYPE_F32); // Saturation limits
+    psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits
+    for (int i = 0; i < images->n; i++) {
+	saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
+	bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
+    }
+
+    // Rescale the images
+    (void)stacRescale(images, errors, NULL, scales, offsets);
+
+    // Combine the images
+    psImage *combined = NULL;		// Combined image
+    psArray *rejected = NULL;		// Array of rejection masks
+    stacCombine(&combined, &rejected, images, errors, config->nReject, NULL, saturated, bad, config->reject);
+
+    psFits *outFile = psFitsAlloc(config->outName);
+    if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
+	psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName);
+    }
+    psTrace("combine", 1, "Combined image written to %s\n", config->outName);
+    psFree(outFile);
+
+    // Clean up
+    psFree(combined);
+    psFree(rejected);
+    psFree(images);
+    psFree(errors);
+    psFree(saturated);
+    psFree(bad);
+    psFree(scales);
+    psFree(offsets);
+    combineConfigFree(config);
+
+#ifdef TESTING
+    // Check memory for leaks, corruption
+    stacCheckMemory();
+#endif
+
+}
Index: trunk/stac/src/combineConfig.c
===================================================================
--- trunk/stac/src/combineConfig.c	(revision 3680)
+++ trunk/stac/src/combineConfig.c	(revision 3680)
@@ -0,0 +1,161 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <string.h>
+#include "pslib.h"
+#include "combineConfig.h"
+
+void help(const char *programName
+    )
+{
+    fprintf (stderr, "shift: shift an image, given the transformation\n"
+	     "Usage: %s [-h] [-v] [-g GAIN] [-r READNOISE] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k SIGMAREJ] [-n NREJECT] OUT IN1 IN2...\n"
+	     "where\n"
+	     "\t-h           Help (this info)\n"
+	     "\t-v           Increase verbosity level\n"
+	     "\t-g GAIN      Gain in e/ADU (1.0)\n"
+	     "\t-r READNOISE Read noise in e (0.0)\n"
+	     "\t-s SAT       Saturation point (65535)\n"
+	     "\t-b BAD       Bad level (0)\n"
+	     "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
+	     "\t-a APER      Aperture radius for photometry (3.0)\n"
+	     "\t-k SIGMAREJ  k-sigma rejection threshold (3.0)\n"
+	     "\t-n NREJECT   Number of rejection iterations (1)\n"
+	     "\tOUT          Output image\n"
+	     "\tIN1 IN2...   Input images (identical size)\n",
+	     programName
+	);
+}
+
+
+combineConfig *combineConfigAlloc(void)
+{
+    combineConfig *config = psAlloc(sizeof(combineConfig)); // Configuration
+
+    // Parameters with default values
+    config->verbose = 0;		// Verbosity level
+    config->gain = 1.0;			// Gain (e/ADU)
+    config->readnoise = 0.0;		// Read noise (e)
+    config->reject = 4.0;		// Rejection threshold (sigma)
+    config->nReject = 1;		// Number of rejection iterations
+    config->saturated = 65535.0;	// Saturation level
+    config->bad = 0.0;			// Bad level
+    config->outName = NULL;		// Output name
+    config->inNames = NULL;		// Input names;
+    config->starFile = NULL;		// Filename of file containing stars
+    config->starMap = NULL;		// Map for stars
+    config->aper = 3.0;			// Aperture for photometry
+
+    return config;
+}
+
+void combineConfigFree(combineConfig *config)
+{
+    psFree(config->inNames);
+    psFree(config);
+}
+
+
+combineConfig *combineParseConfig(int argc, char *argv[])
+{
+    combineConfig *config = combineConfigAlloc(); // Configuration
+
+    const char *programName = argv[0];	// Program name
+
+    /* Variables for getopt */
+    int opt;   /* Option, from getopt */
+    extern char *optarg;   /* Argument accompanying switch */
+    extern int optind;   /* getopt variables */
+
+    /* Parse command-line arguments using getopt */
+    while ((opt = getopt(argc, argv, "hvg:r:s:b:p:a:k:n:")) != -1) {
+        switch (opt) {
+	  case 'h':
+	    help(programName);
+	    exit(EXIT_SUCCESS);
+	  case 'v':
+            config->verbose++;
+            break;
+	  case 'r':
+            if (sscanf(optarg, "%f", &config->readnoise) != 1) {
+		printf("Unable to read readnoise.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 'g':
+	    if (sscanf(optarg, "%f", &config->gain) != 1) {
+ 		printf("Unable to read gain.\n");
+		help(programName);
+		exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 's':
+            if (sscanf(optarg, "%f", &config->saturated) != 1) {
+ 		printf("Unable to read saturation limit.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 'b':
+            if (sscanf(optarg, "%f", &config->bad) != 1) {
+ 		printf("Unable to read bad limit.\n");
+		help(programName);
+                exit(EXIT_FAILURE);
+	    }
+	    break;
+	  case 'p':
+	    if (argc < optind+1) {
+ 		printf("Unable to read photometric files.\n");
+		help(programName);
+		exit(EXIT_FAILURE);
+	    }
+	    config->starFile = argv[optind-1];
+	    config->starMap = argv[optind++];
+	    // Note: incrementing optind, so I can read more than one parameter.
+	    break;
+	  case 'a':
+            if (sscanf(optarg, "%f", &config->aper) != 1) {
+ 		printf("Unable to read aperture.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+	    break;
+	  case 'k':
+            if (sscanf(optarg, "%f", &config->reject) != 1) {
+ 		printf("Unable to read rejection limit.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+	    break;
+	  case 'n':
+            if (sscanf(optarg, "%d", &config->nReject) != 1) {
+  		printf("Unable to read number of rejection iterations.\n");
+		help(programName);
+                exit(EXIT_FAILURE);
+            }
+	    break;
+	  default:
+	    printf("Bad option: %c\n", opt);
+	    help(programName);
+	    exit(EXIT_FAILURE);
+	}
+    }
+
+    /* Get the remaining, mandatory values */
+    argc -= optind;
+    argv += optind;
+    if (argc < 2) {
+        help(programName);
+        exit(EXIT_FAILURE);
+    }
+    config->outName = argv[0];		// Output filename
+    config->inNames = psArrayAlloc(argc-1); // Input filenames
+    for (int i = 1; i < argc; i++) {
+	config->inNames->data[i-1] = psAlloc(strlen(argv[i]));
+	strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i]));
+    }
+
+    return config;
+}
Index: trunk/stac/src/combineConfig.h
===================================================================
--- trunk/stac/src/combineConfig.h	(revision 3680)
+++ trunk/stac/src/combineConfig.h	(revision 3680)
@@ -0,0 +1,28 @@
+#ifndef COMBINE_CONFIG_H
+#define COMBINE_CONFIG_H
+
+typedef struct {
+    int verbose;			// Verbosity level
+    float gain;				// Gain (e/ADU)
+    float readnoise;			// Read noise (e)
+    float reject;			// Rejection threshold (sigma)
+    int nReject;			// Number of rejection iterations
+    float saturated;			// Saturation level
+    float bad;				// Bad level
+    char *outName;			// Output name
+    psArray *inNames;			// Input names;
+    char *starFile;			// Filename of file containing stars
+    char *starMap;			// Map for stars
+    float aper;				// Aperture for photometry
+} combineConfig;
+
+
+combineConfig *combineConfigAlloc(void);// Allocator
+void combineConfigFree(combineConfig *config); // Deallocator
+
+void help(const char *programName);	// Print help
+combineConfig *combineParseConfig(int argc, char *argv[]); // Parse command line
+
+
+
+#endif
Index: trunk/stac/src/stacRead.c
===================================================================
--- trunk/stac/src/stacRead.c	(revision 3673)
+++ trunk/stac/src/stacRead.c	(revision 3680)
@@ -32,4 +32,8 @@
 	    psFree(image);
 	} else {
+	    int numNaN = psImageClipNaN(image, 0.0);
+	    if (numNaN) {
+		psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN);
+	    }
 	    images->data[i] = image;
 	}
Index: trunk/stac/src/stacWrite.c
===================================================================
--- trunk/stac/src/stacWrite.c	(revision 3680)
+++ trunk/stac/src/stacWrite.c	(revision 3680)
@@ -0,0 +1,77 @@
+#include <stdio.h>
+#include <assert.h>
+#include "pslib.h"
+#include "stac.h"
+
+bool stacWriteMap(const char *mapName,	// Filename to write to
+		  psPlaneTransform *map	// Map to write
+    )
+{
+    assert(mapName);
+
+    FILE *mapFile = fopen(mapName, "w");
+    if (!mapFile) {
+	fprintf(stderr, "Unable to open map file: %s\n", mapName);
+	return false;
+    }
+
+    psDPolynomial2D *xMap = map->x;	// x transform
+    psDPolynomial2D *yMap = map->y;	// y transform
+
+    // A crucial limitation of the current system --- the order of each polynomial must be the same
+    assert(xMap->nX == xMap->nY && yMap->nX == yMap->nY && xMap->nX == yMap->nX);
+    int order = xMap->nX - 1;	// The polynomial order
+    fprintf(mapFile, "%d\n", order);
+    
+    // x coefficients
+    for (int k = 0; k < order + 1; k++) {
+	for (int j = 0; j < k + 1; j++) {
+	    int i = k - j;
+	    if (xMap->mask[i][j]) {
+		fprintf(mapFile, "0.0 ");
+	    } else {
+		fprintf(mapFile, "%g ", xMap->coeff[i][j]);
+	    }
+	}
+    }
+    fprintf(mapFile, "\n");
+
+    // y coefficients
+    for (int k = 0; k < order + 1; k++) {
+	for (int j = 0; j < k + 1; j++) {
+	    int i = k - j;
+	    if (yMap->mask[i][j]) {
+		fprintf(mapFile, "0.0 ");
+	    } else {
+		fprintf(mapFile, "%g ", yMap->coeff[i][j]);
+	    }
+	}
+    }
+    fprintf(mapFile, "\n");
+    
+    fclose(mapFile);
+
+    return true;
+}
+
+
+bool stacWriteMaps(const psArray *names, // Filenames of the input images (will add ".map")
+		   const psArray *maps	// Maps to write
+    )
+{
+    assert(names);
+    assert(maps);
+    assert(names->n == maps->n);
+
+    for (int i = 0; i < names->n; i++) {
+ 	char mapName[MAXCHAR];		// Filename of error image
+	sprintf(mapName, "%s.map", names->data[i]);
+	if (!stacWriteMap(mapName, maps->data[i])) {
+	    return false;
+	}
+    }
+
+    return true;
+}
+
+