Index: unk/stac/src/Makefile.median
===================================================================
--- /trunk/stac/src/Makefile.median	(revision 3669)
+++ 	(revision )
@@ -1,35 +1,0 @@
-SHELL = /bin/sh
-CC = gcc
-CFLAGS += -g -pg -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE -DCRFLUX -DTESTING
-PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lcfitsio -lm -lxml2 -lmysqlclient 
-LDFLAGS += $(PSLIB)
-
-TARGET = median
-
-.PHONY: tags clean empty test profile optimise
-
-.c.o:
-		$(CC) -c $(CFLAGS) $(OPTFLAGS) $<
-
-median:		median.o
-		$(CC) $(CFLAGS) -o $@ median.o $(LDFLAGS) $(OPTFLAGS)
-
-clean:
-		-$(RM) *.o gmon.* profile.txt
-
-empty:		clean
-		-$(RM) $(TARGET) TAGS
-
-test:		median
-		./median test_3.fits.shift.1
-
-# Run profiling.
-profile:	CFLAGS += -pg
-profile:	empty $(TARGET)
-		$(RM) gmon.*
-		for ((i = 0; i < 20; i++)) \
-		do \
-			$(MAKE) -f Makefile.median test; \
-			mv -f gmon.out gmon.$$i; \
-		done
-		gprof -p -q -l -c $(TARGET) gmon.* > profile.txt
Index: unk/stac/src/fitTrans.c
===================================================================
--- /trunk/stac/src/fitTrans.c	(revision 3669)
+++ 	(revision )
@@ -1,183 +1,0 @@
-#include <stdio.h>
-#include "pslib.h"
-#include "fitTrans.h"
-
-int main(int argc, char **argv)
-{
-    // Set trace levels
-    psTraceSetLevel(".", 0);
-    psTraceSetLevel("fitTrans.config", 10);
-    psTraceSetLevel("fitTrans.read", 10);
-    psTraceSetLevel("stac.checkMemory", 10);
-
-    // Set logging level
-    psLogSetLevel(9);
-
-    // Parse command line
-    fitTransConfig *config = fitTransParseConfig(argc, argv);
-
-    printf("Order: %d\n", config->order);
-
-    // Read the input files
-    psVector *xReference = NULL;
-    psVector *yReference = NULL;
-    psVector *xInput = NULL;
-    psVector *yInput = NULL;
-    (void)fitTransRead(&xReference, &yReference, config->reference);
-    (void)fitTransRead(&xInput, &yInput, config->input);
-    int nStars = xReference->n;		// Number of stars
-    if (nStars != xInput->n) {
-	fprintf(stderr, "Number of stars in reference (%d) and input lists (%d) do not match.\n",
-		xReference->n, xInput->n);
-	exit(EXIT_FAILURE);
-    }
-
-    psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
-    psPlane *outCoord = psAlloc(sizeof(psPlane)); // Output Coordinates
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics
-    psVector *mask = psVectorAlloc(nStars, PS_TYPE_U8);	// Mask for stars
-    for (int s = 0; s < nStars; s++) {
-	mask->data.U8[s] = 0;
-    }
-
-#if 0
-    for (int i = 0; i < nStars; i++) {
-	xReference->data.F32[i] = (xReference->data.F32[i] - 512.0)/512.0;
-	yReference->data.F32[i] = (yReference->data.F32[i] - 512.0)/512.0;
-	xInput->data.F32[i] = (xInput->data.F32[i] - 512.0)/512.0;
-	yInput->data.F32[i] = (yInput->data.F32[i] - 512.0)/512.0;
-    }
-#endif
-
-    // Solve for the transformation
-    //psPlaneTransform *transform = fitTransSolve(NULL, xReference, yReference, xInput, yInput, NULL, config);
-    psPlaneTransform *transform = fitTransSolve(NULL, xReference, yReference, xInput, yInput, mask, config);
-
-    // Iterate to do rejection
-    bool keepGoing = true;
-    for (int i = 0; i < config->nIter && keepGoing; i++) {
-
-	psVector *xDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in x
-	psVector *yDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in y
-	for (int s = 0; s < nStars; s++) {
-	    inCoord->x = xInput->data.F32[s];
-	    inCoord->y = yInput->data.F32[s];
-	    (void)psPlaneTransformApply(outCoord, transform, inCoord);
-	    xDiff->data.F32[s] = xReference->data.F32[s] - outCoord->x;
-	    yDiff->data.F32[s] = yReference->data.F32[s] - outCoord->y;
-	    //printf("%f %f\n", xDiff->data.F32[s], yDiff->data.F32[s]);
-	}
-
-	(void)psVectorStats(stats, xDiff, mask, 0);
-	float xSD = stats->sampleStdev;	// Standard deviation in x
-	float xMean = stats->sampleMean;
-	(void)psVectorStats(stats, yDiff, mask, 0);
-	float ySD = stats->sampleStdev;	// Standard deviation in y
-	float yMean = stats->sampleMean;
-	printf("Mean is %f, %f\n", xMean, yMean);
-	printf("Standard deviation is %f, %f\n", xSD, ySD);
-
-	// Do the clipping
-	int nRej = 0;			// Number rejected
-	for (int s = 0; s < nStars; s++) {
-	    if (!mask->data.U8[s] && ((fabsf(xDiff->data.F32[s]) > config->nSigma * xSD) ||
-				      (fabsf(yDiff->data.F32[s]) > config->nSigma * ySD))) {
-		mask->data.U8[s] = 1;
-		nRej++;
-		printf("     Clipping star %d\n", s);
-	    }
-	}
-	printf("%d stars clipped on iteration %d.\n", nRej, i + 1);
-	psFree(xDiff);
-	psFree(yDiff);
-
-	if (nRej == 0) {
-	    keepGoing = false;		// Nothing more to reject
-	} else {
-	    fitTransSolve(transform, xReference, yReference, xInput, yInput, mask, config);
-	}
-    }
-
-    psVector *xDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in x
-    psVector *yDiff = psVectorAlloc(xReference->n, PS_TYPE_F32); // Difference vector in y
-    for (int s = 0; s < nStars; s++) {
-	inCoord->x = xInput->data.F32[s];
-	inCoord->y = yInput->data.F32[s];
-	(void)psPlaneTransformApply(outCoord, transform, inCoord);
-	xDiff->data.F32[s] = xReference->data.F32[s] - outCoord->x;
-	yDiff->data.F32[s] =yReference->data.F32[s] - outCoord->y;
-    }
-    
-    (void)psVectorStats(stats, xDiff, mask, 0);
-    float xSD = stats->sampleStdev;	// Standard deviation in x
-    float xMean = stats->sampleMean;	// Mean of x
-    (void)psVectorStats(stats, yDiff, mask, 0);
-    float ySD = stats->sampleStdev;	// Standard deviation in y
-    float yMean = stats->sampleMean;	// Mean of y
-    printf("Final mean is %f, %f\n", xMean, yMean);
-    printf("Final standard deviation is %f, %f\n", xSD, ySD);
-    psFree(xDiff);
-    psFree(yDiff);
-
-    psFree(inCoord);
-    psFree(outCoord);
-    psFree(mask);
-    psFree(stats);
-
-    int order = config->order + 1;	// Polynomial order
-    int nCoeff = order * (order + 1) / 2; // Number of coefficients
-    printf("x' = ");
-    for (int i = 0; i < order; i++) {
-	for (int j = 0; j < order - i; j++) {
-	    if (i == 0 && j == 0) {
-		printf ("%f\n", transform->x->coeff[0][0]);
-	    } else {
-		printf ("     + %f x^%d y^%d\n", transform->x->coeff[i][j], i, j);
-	    }
-	}
-    }
-    printf("\n");
-    printf("y' = ");
-    for (int i = 0; i < order; i++) {
-	for (int j = 0; j < order - i; j++) {
-	    if (i == 0 && j == 0) {
-		 printf ("%f\n", transform->y->coeff[0][0]);
-	    } else {
-		printf ("     + %f x^%d y^%d\n", transform->y->coeff[i][j], i, j);
-	    }
-	}
-    }
-
-    if (config->mapFile != NULL) {
-	FILE *fp = fopen(config->mapFile, "w");	// File pointer
-	if (fp == NULL) {
-	    fprintf(stderr, "Unable to open %s for writing.\n", fp);
-	    exit(EXIT_FAILURE);
-	} else {
-	    fprintf(fp, "%d\n", config->order);
-	    for (int i = 0; i < order; i++) {
-		for (int j = 0; j < order - i; j++) {
-		    fprintf(fp, "%.10f ", transform->x->coeff[i][j]);
-		}
-	    }
-	    fprintf(fp, "\n");
-	    for (int i = 0; i < order; i++) {
-		for (int j = 0; j < order - i; j++) {
-		    fprintf(fp, "%.10f ", transform->y->coeff[i][j]);
-		}
-	    }
-	    fprintf(fp, "\n");
-	    fclose(fp);
-	}
-    }
-
-    // Clean up
-    psFree(xReference);
-    psFree(yReference);
-    psFree(xInput);
-    psFree(yInput);
-    psFree(transform);
-    psFree(config);
-
-    stacCheckMemory();
-}
Index: unk/stac/src/fitTrans.h
===================================================================
--- /trunk/stac/src/fitTrans.h	(revision 3669)
+++ 	(revision )
@@ -1,77 +1,0 @@
-#ifndef FIT_TRANS_H
-#define FIT_TRANS_H
-
-
-// Temporary implementation
-psPlaneTransform *psPlaneTransformAlloc(psS32 n1, psS32 n2);
-
-/************************************************************************************************************/
-
-// Configuration for the program
-typedef struct {
-    int verbose;			// Verbosity level
-    int order;				// Polynomial order of the fit
-    float nSigma;			// Number of standard deviations at which to clip
-    int nIter;				// Number of rejection iterations
-    const char *reference;		// Name of reference star list
-    const char *input;			// Name of input star list
-    const char *mapFile;		// Name of output map file
-} fitTransConfig;
-
-
-/************************************************************************************************************/
-
-// fitTransConfig.c
-
-// The help message
-void help(const char *name
-    );
-
-// Parse the command line arguments
-fitTransConfig *fitTransParseConfig(int argc, // Number of arguments
-				    char **argv	// Arguments
-    );
-
-/************************************************************************************************************/
-
-// fitTransRead.c
-
-// Read an input file
-bool fitTransRead(psVector **x,		// x positions
-		  psVector **y,		// y positions
-		  const char *filename	// Filename from which to read
-    );
-
-/************************************************************************************************************/
-
-// stacCheckMemory.c
-// This stuff was written for STAC, but it should work anywhere, so I'll just plug it in here.
-
-// Check memory
-void stacCheckMemory(void);
-
-// Print out the problem
-void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block.
-		       const char *file, ///< the file in which the problem originated
-		       psS32 lineno	///< the line number in which the problem originated
-    );
-
-// Print out a memblock when it's allocated --- this function used as a callback
-psMemoryId stacMemPrint(const psMemBlock *ptr);
-
-/************************************************************************************************************/
-
-// fitTransSolve.c
-
-// Solve for the transformation
-psPlaneTransform *fitTransSolve(psPlaneTransform *trans,// Output transformation, or NULL
-				psVector *xReference, psVector *yReference, // Reference coordinates
-				psVector *xInput, psVector *yInput, // Input coordinates
-				psVector *mask,	// Mask for coordinates
-				fitTransConfig *config // Configuration
-    );
-
-/************************************************************************************************************/
-
-
-#endif
Index: unk/stac/src/fitTransConfig.c
===================================================================
--- /trunk/stac/src/fitTransConfig.c	(revision 3669)
+++ 	(revision )
@@ -1,92 +1,0 @@
-#include <stdio.h>
-#include <unistd.h>
-#include <getopt.h>
-#include "pslib.h"
-#include "fitTrans.h"
-
-
-void help(const char *name		// Program name
-    )
-{
-    fprintf (stderr, "fitTrans: Fit image transformation given two matched star lists.\n"
-	     "Usage: %s [-h] [-v] [-o ORDER] [-s SIGMA] [-n NITER] [-m MAPFILE] REFERENCE INPUT\n"
-	     "where\n"
-	     "\t-h           Help (this info)\n"
-	     "\t-v           Increase verbosity level\n"
-	     "\t-o ORDER     Order for polynomial fit: 1, 2, 3 (1)\n"
-	     "\t-s SIGMA     Number of standard deviations at which to reject (3)\n"
-	     "\t-n NITER     Number of rejection iterations (1)\n"
-	     "\t-m MAPFILE   Name of file to which to write the transformation map\n"
-	     "\tREFERENCE    Reference list: index, x, y, magnitude\n"
-	     "\tINPUT        Input list: index, x, y, magnitude\n",
-	     name
-	);
-}
-
-
-fitTransConfig *fitTransParseConfig(int argc, // Number of arguments
-				    char **argv	// Arguments
-    )
-{
-    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 variable
-
-    // Initialise the configuration
-    fitTransConfig *config = psAlloc(sizeof(fitTransConfig));
-    config->verbose = 0;		// Verbosity level
-    config->order = 1;			// Order for fit
-    config->nSigma = 3.0;		// Number of sigma at which to clip
-    config->nIter = 1;			// Number of rejection iterations
-    config->mapFile = NULL;		// Name of map file
-
-    /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvo:s:n:m:")) != -1) {
-        switch (opt) {
-	  case 'h':
-	    help(programName);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
-            config->verbose++;
-            break;
-	  case 'o':
-	    if (sscanf(optarg, "%d", &config->order) != 1) {
-		help(programName);
-	    }
-	    break;
-	  case 's':
-	    if (sscanf(optarg, "%f", &config->nSigma) != 1) {
-		help(programName);
-	    }
-	    break;
-	  case 'n':
-	    if (sscanf(optarg, "%d", &config->nIter) != 1) {
-		help(programName);
-	    }
-	    break;
-	  case 'm':
-	    config->mapFile = optarg;
-	    break;
-	  default:
-	    help(programName);
-	}
-    }
-
-    argc -= optind;
-    argv += optind;
-    if (argc != 2) {
-	help(programName);
-        exit(EXIT_FAILURE);
-    }
-
-    config->reference = argv[0];	// Reference list
-    config->input = argv[1];		// Input list
-
-    psTrace("fitTrans.config", 1, "Reference file is %s\n", config->reference);
-    psTrace("fitTrans.config", 1, "Input file is %s\n", config->input);
-    
-    return config;
-}
Index: unk/stac/src/fitTransRead.c
===================================================================
--- /trunk/stac/src/fitTransRead.c	(revision 3669)
+++ 	(revision )
@@ -1,41 +1,0 @@
-#include <stdio.h>
-#include "pslib.h"
-#include "fitTrans.h"
-
-#define BUFFER_SIZE 100
-
-bool fitTransRead(psVector **x,		// x positions
-		  psVector **y,		// y positions
-		  const char *filename	// Filename from which to read
-    )
-{
-    FILE *fp;				// File pointer
-
-    fp = fopen(filename,"r");
-    if (fp == NULL) {
-	psError("fitTrans.read", "Unable to open file %s\n", filename);
-	return false;
-    }
-
-    int index = 0;			// Index of star
-    float magnitude = 0.0;		// Magnitude of star
-    int num = 0;			// Number of stars read
-
-    *x = psVectorAlloc(BUFFER_SIZE, PS_TYPE_F32);
-    *y = psVectorAlloc(BUFFER_SIZE, PS_TYPE_F32);
-
-    while (fscanf(fp, "%d %f %f %f", &index, &((*x)->data.F32[num]), &((*y)->data.F32[num]), &magnitude) != EOF) {
-	num++;
-	if (num % BUFFER_SIZE) {
-	    (void)psVectorRealloc(*x, num+BUFFER_SIZE);
-	    (void)psVectorRealloc(*y, num+BUFFER_SIZE);
-	}
-    }
-
-    (*x)->n = num;
-    (*y)->n = num;
-
-    psTrace("fitTrans.read", 1, "Read %d entries from %s\n", num, filename);
-
-    return true;
-}
Index: unk/stac/src/fitTransSolve.c
===================================================================
--- /trunk/stac/src/fitTransSolve.c	(revision 3669)
+++ 	(revision )
@@ -1,127 +1,0 @@
-#include <stdio.h>
-#include <assert.h>
-#include "pslib.h"
-#include "fitTrans.h"
-
-
-psPlaneTransform *fitTransSolve(psPlaneTransform *trans,// Output coefficients, or NULL
-				psVector *xReference, psVector *yReference, // Reference coordinates
-				psVector *xInput, psVector *yInput, // Input coordinates
-				psVector *mask,	// Mask for coordinates
-				fitTransConfig *config // Configuration
-    )
-{
-    int order = config->order + 1;	// Order of polynomials
-    int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
-
-    assert(!trans || ((trans->x->nX == order) && (trans->x->nY == order) && (trans->y->nX == order) &&
-		      (trans->y->nY == order)));
-    assert(xReference);
-    assert(yReference);
-    assert(xReference->n == yReference->n);
-    assert(xInput);
-    assert(yInput);
-    assert(xInput->n == yInput->n);
-    assert(xReference->n == xInput->n);
-    assert(!mask || xReference->n == mask->n);
-
-    // Check inputs
-    if (trans == NULL) {
-	trans = psPlaneTransformAlloc(order, order);
-	// Initialise transformation
-	for (int i = 0; i < order; i++) {
-	    for (int j = 0; j < order; j++) {
-		trans->x->coeff[i][j] = 0.0;
-		trans->y->coeff[i][j] = 0.0;
-		if (i + j < order) {
-		    trans->x->mask[i][j] = 0;
-		    trans->y->mask[i][j] = 0;
-		} else {
-		    trans->x->mask[i][j] = 1;
-		    trans->y->mask[i][j] = 1;
-		}		    
-	    }
-	}
-    }
-
-    // Initialise calculations
-    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
-    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
-    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
-    for (int i = 0; i < nCoeff; i++) {
-	for (int j = 0; j < nCoeff; j++) {
-	    matrix->data.F64[i][j] = 0.0;
-	}
-	xVector->data.F64[i] = 0.0;
-	yVector->data.F64[i] = 0.0;
-    }
-
-    // Create fake polynomial to use in evaluation
-    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
-    for (int i = 0; i < order; i++) {
-	for (int j = 0; j < order; j++) {
-	    fakePoly->coeff[i][j] = 1.0;// Set all coeffecients to 1
-	    fakePoly->mask[i][j] = 1;	// Mask all coefficients; unmask to evaluate
-	}
-    }
-
-    // Iterate over stars that haven't been masked
-    int numStars = 0;			// Number of stars used
-    for (int k = 0; k < xReference->n; k++) {
-	if (!mask || ! mask->data.U8[k]) {
-	    numStars++;
-
-	    // Iterate over the polynomial coefficients
-	    for (int i = 0, ijIndex = 0; i < order; i++) {
-		for (int j = 0; j < order - i; j++, ijIndex++) {
-
-		    fakePoly->mask[i][j] = 0;
-		    double ijPoly = psDPolynomial2DEval((double)xInput->data.F32[k], (double)yInput->data.F32[k], fakePoly);
-		    fakePoly->mask[i][j] = 1;
-		    
-		    for (int m = 0, mnIndex = 0; m < order; m++) {
-			for (int n = 0; n < order - m; n++, mnIndex++) {
-			    
-			    fakePoly->mask[m][n] = 0;
-			    double mnPoly = psDPolynomial2DEval((double)xInput->data.F32[k], (double)yInput->data.F32[k], fakePoly);
-			    fakePoly->mask[m][n] = 1;
-
-			    matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
-			}
-		    }
-		    
-		    xVector->data.F64[ijIndex] += ijPoly * (double)xReference->data.F32[k];
-		    yVector->data.F64[ijIndex] += ijPoly * (double)yReference->data.F32[k];
-		}
-	    } // Iterating over coefficients
-	}
-    } // Iterating over stars
-
-    printf("Solving based on %d stars...\n", numStars);
-
-    // Solution via LU Decomposition
-    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64);	// Permutation vector for LU Decomposition
-    psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix);	// LU decomposed matrix
-    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
-    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
-    psFree(permutation);
-    psFree(luMatrix);
-
-    // Stuff coefficients into transformation
-    for (int i = 0, ijIndex = 0; i < order; i++) {
-	for (int j = 0; j < order - i; j++, ijIndex++) {
-	    trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
-	    trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
-	}
-    }
-
-    // Clean up
-    psFree(xSolution);
-    psFree(ySolution);
-    psFree(matrix);
-    psFree(xVector);
-    psFree(yVector);
-    psFree(fakePoly);
-
-    return trans;
-}
Index: unk/stac/src/median.c
===================================================================
--- /trunk/stac/src/median.c	(revision 3669)
+++ 	(revision )
@@ -1,59 +1,0 @@
-#include <stdio.h>
-#include "pslib.h"
-
-double getTime(void)
-/* Gets the current time.  Got this from Nick Kaiser's fetchpix.c */
-{
-    struct timeval tv;
-    gettimeofday(&tv, NULL);
-    return(tv.tv_sec + 1.e-6 * tv.tv_usec);
-}
-
-int main(int argc, char *argv[])
-{
-#if 0
-    psTraceSetLevel(".", 10);
-    psTraceSetLevel(".psLib", 10);
-    psTraceSetLevel(".psLib.dataManip", 10);
-    psTraceSetLevel(".psLib.dataManip.psFunctions", 10);
-    psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0);
-    psTraceSetLevel(".psLib.dataManip.psFunctions.vectorBinDisect", 0);
-    psTraceSetLevel(".pslib.stats.robust", 10);
-#endif
-
-    if (argc != 2) {
-	exit(EXIT_FAILURE);
-    }
-
-    const char *filename = argv[1];
-
-    psFits *imageFile = psFitsAlloc(filename);
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
-    psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
-    if (image == NULL) {
-	psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", filename);
-	exit(EXIT_FAILURE);
-    }
-    psTrace("median", 4, "Image %s is %dx%d\n", filename, image->numCols, image->numRows);
-    psFree(imageRegion);
-    // Convert to 32-bit floating point, in necessary
-    if (image->type.type != PS_TYPE_F32) {
-	psTrace("stac.read.images", 3, "Converting %s to floating point in memory....", filename);
-	psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
-	psFree(image);
-	image = temp;
-    }
-    psFree(imageFile);
-
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-    double startTime = getTime();
-    stats = psImageStats(stats, image, NULL, 0);
-    double endTime = getTime();
-    printf("Median: %f\n", stats->robustMedian);
-    printf("Time: %f\n", endTime - startTime);
-
-    psFree(stats);
-    psFree(image);
-
-    exit(EXIT_SUCCESS);
-}
Index: unk/stac/src/photScale.c
===================================================================
--- /trunk/stac/src/photScale.c	(revision 3669)
+++ 	(revision )
@@ -1,341 +1,0 @@
-#include <stdio.h>
-#include <getopt.h>
-#include <math.h>
-#include "pslib.h"
-#include "psMinimize.h"
-
-#define BUFFER 100			// Buffer size for star list
-#define MAXITER 100			// Maximum iterations for fitting
-#define TOLERANCE 1e-2			// Fitting tolerance
-#define ZP 25				// Zero point magnitude
-
-#define SQUARE(x) ((x)*(x))
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-
-
-psVector *gaussian2d_deriv(psImage *deriv, // Derivatives wrt the parameters
-			   const psVector *params,// Parameters with which to evaluate
-			   const psArray *position // Position at which to evaluate
-    )
-{
-    float bg = params->data.F32[0];
-    float norm = params->data.F32[1];
-    float x0 = params->data.F32[2];
-    float y0 = params->data.F32[3];
-    float sigmaU = params->data.F32[4];
-    float sigmaV = params->data.F32[5];
-    float theta = params->data.F32[6];
-
-    psVector *gauss = psVectorAlloc(position->n, PS_TYPE_F32);
-    for (int i = 0; i < position->n; i++) {
-	psVector *pos = position->data[i];
-
-	float x = pos->data.F32[0];
-	float y = pos->data.F32[1];
-
-	float xSub = x - x0;
-	float ySub = y - y0;
-
-	float cosTheta = cos(theta);
-	float sinTheta = sin(theta);
-
-	float invSigmaU2 = 1.0 / SQUARE(sigmaU);
-	float invSigmaV2 = 1.0 / SQUARE(sigmaV);
-
-	float u = - xSub * cosTheta + ySub * sinTheta;
-	float v = xSub * cosTheta + ySub * sinTheta;
-
-	float uOnSigmaU2 = u * invSigmaU2;
-	float vOnSigmaV2 = v * invSigmaV2;
-
-	float exponent = SQUARE(u) * invSigmaU2 + SQUARE(v) * invSigmaV2;
-	float bell = expf(- 0.5 * exponent);
-	gauss->data.F32[i] = norm * bell + bg;
-
-	float deriv_bg = 1.0;
-	float deriv_norm = bell;
-	float deriv_x0 = norm * bell * (- uOnSigmaU2 + vOnSigmaV2) * cosTheta;
-	float deriv_y0 = norm * bell * (uOnSigmaU2 + vOnSigmaV2) * sinTheta;
-	float deriv_sigmaU = norm * bell * uOnSigmaU2 * u / sigmaU;
-	float deriv_sigmaV = norm * bell * vOnSigmaV2 * v / sigmaV;
-	float deriv_theta = - norm * bell * (uOnSigmaU2 * (xSub * sinTheta + ySub * cosTheta) +
-					     vOnSigmaV2 * (-xSub * sinTheta + ySub * cosTheta));
-
-	deriv->data.F32[i][0] = deriv_bg;
-	deriv->data.F32[i][1] = deriv_norm;
-	deriv->data.F32[i][2] = deriv_x0;
-	deriv->data.F32[i][3] = deriv_y0;
-	deriv->data.F32[i][4] = deriv_sigmaU;
-	deriv->data.F32[i][5] = deriv_sigmaV;
-	deriv->data.F32[i][6] = deriv_theta;
-    }
-
-    return gauss;
-}    
-
-psVector *gaussian2d_noderiv(const psVector *params,// Parameters with which to evaluate
-			     const psArray *position // Position at which to evaluate
-			     )
-{
-    float bg = params->data.F32[0];
-    float norm = params->data.F32[1];
-    float x0 = params->data.F32[2];
-    float y0 = params->data.F32[3];
-    float sigmaU = params->data.F32[4];
-    float sigmaV = params->data.F32[5];
-    float theta = params->data.F32[6];
-
-    psVector *gauss = psVectorAlloc(position->n, PS_TYPE_F32);
-    for (int i = 0; i < position->n; i++) {
-	psVector *pos = position->data[i];
-
-	float x = pos->data.F32[0];
-	float y = pos->data.F32[1];
-
-	float xSub = x - x0;
-	float ySub = y - y0;
-
-	float cosTheta = cos(theta);
-	float sinTheta = sin(theta);
-
-	float invSigmaU2 = 1.0 / SQUARE(sigmaU);
-	float invSigmaV2 = 1.0 / SQUARE(sigmaV);
-
-	float u = - xSub * cosTheta + ySub * sinTheta;
-	float v = xSub * cosTheta + ySub * sinTheta;
-
-	float uOnSigmaU2 = u * invSigmaU2;
-	float vOnSigmaV2 = v * invSigmaV2;
-
-	float exponent = SQUARE(u) * invSigmaU2 + SQUARE(v) * invSigmaV2;
-	float bell = expf(- 0.5 * exponent);
-	gauss->data.F32[i] = norm * bell + bg;
-
-    }
-
-    return gauss;
-}
-
-
-
-void help(char *programName)
-{
-    fprintf(stderr, "photScale: calculates the photometric scaling between two images\n"
-	    "Usage: %s [-r REJECT] [-i NUM_ITER] REFERENCE INPUT STARS\n"
-	    "where\n"
-	    "\t-r REJECT       Rejection limit, standard deviations (2.5)\n"
-	    "\t-i NUM_ITER     Number of rejection iterations (3)\n"
-	    "\tREFERENCE       Reference image\n"
-	    "\tINPUT           Input image\n"
-	    "\tSTARS           File with stellar positions, x y on each line\n",
-	    programName);
-}
-
-
-int main(int argc, char *argv[])
-{
-    int verbose = 0;			// Verbosity level
-    int numIter = 3;			// Number of iterations
-    float reject = 2.5;			// Rejection limit (sigma)
-    int size = 5;			// Size of stamps
-
-    psTraceSetLevel(".", 10);
-    psTraceSetLevel("read", 10);
-    psTraceSetLevel("phot", 10);
-    psTraceSetLevel("gauss", 9);
-    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
-
-
-    /* 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, "hvr:i:s:")) != -1) {
-        switch (opt) {
-	  case 'h':
-	    help(argv[-1]);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
-	    verbose++;
-	    break;
-	  case 'r':
-	    if (sscanf(optarg, "%f", &reject) != 1) {
-		help(argv[-1]);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'i':
-	    if (sscanf(optarg, "%d", &numIter) != 1) {
-		help(argv[-1]);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 's':
-	    if (sscanf(optarg, "%d", &size) != 1) {
-		help(argv[-1]);
-		exit(EXIT_FAILURE);
-	    }
-	  default:
-	    help(argv[-1]);
-	    exit(EXIT_FAILURE);
-	}
-    }
-
-    /* Get the remaining, mandatory values */
-    argc -= optind;
-    argv += optind;
-    if (argc != 3) {
-        help(argv[-1]);
-        exit(EXIT_FAILURE);
-    }
-
-    const char *reference = argv[0];
-    const char *input = argv[1];
-    const char *stars = argv[2];
-
-    // Read images
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
-    psFits *referenceFile = psFitsAlloc(reference); // FITS file
-    psImage *refImage = psFitsReadImage(NULL, referenceFile, *imageRegion, 0); // Reference image
-    if (refImage == NULL) {
-	psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", reference);
-	exit(EXIT_FAILURE);
-    }
-    if (refImage->type.type != PS_TYPE_F32) {
-	// Convert to F32
-	psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32);
-	psFree(refImage);
-	refImage = temp;
-    }
-    psFree(referenceFile);
-    psImageClipNaN(refImage, 0.0);
-
-#if 0
-    psFits *inputFile = psFitsAlloc(input); // FITS file
-    psImage *inImage = psFitsReadImage(NULL, inputFile, *imageRegion, 0); // Input image
-    if (inImage == NULL) {
-	psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", input);
-	exit(EXIT_FAILURE);
-    }
-    if (inImage->type.type != PS_TYPE_F32) {
-	// Convert to F32
-	psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32);
-	psFree(inImage);
-	refImage = temp;
-    }
-    psFree(inputFile);
-#endif
-
-    // List of star positions
-    psArray *starList = psArrayAlloc(BUFFER);
-    FILE *starFile = fopen(stars, "r");
-    if (starFile == NULL) {
-	fprintf(stderr, "Unable to open %s\n", stars);
-	exit(EXIT_FAILURE);
-    }
-    float x, y;
-    int num = 0;			// Number of stars
-    while (fscanf(starFile, "%f %f", &x, &y) == 2) {
-	psPlane *coords = psPlaneAlloc();
-	coords->x = x;
-	coords->y = y;
-	starList->data[num] = coords;
-	num++;
-	if (num % BUFFER) {
-	    starList = psArrayRealloc(starList, num + BUFFER);
-	}
-    }
-    starList->n = num;
-    psTrace("read", 1, "%d stars read.\n", num);
-
-    psStats *medianStats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psImageStats(medianStats, refImage, NULL, 0);
-    float refMedian = medianStats->sampleMedian;
-
-    psTrace("phot", 3, "Image median: %f\n", refMedian);
-
-    psVector *refFlux = psVectorAlloc(num, PS_TYPE_F32);
-    psVector *refBG = psVectorAlloc(num, PS_TYPE_F32);
-
-    psVector *inFlux = psVectorAlloc(num, PS_TYPE_F32);
-    psVector *inBG = psVectorAlloc(num, PS_TYPE_F32);
-    psMinimization *min = psMinimizationAlloc(MAXITER, TOLERANCE);
-    for (int i = 0; i < num; i++) {
-	psPlane *coords = starList->data[i];
-
-	// Iterate over postage stamp
-	int xCen = (int)(coords->x + 0.5);
-	int yCen = (int)(coords->y + 0.5);
-	int xMin = MAX(0, xCen - size);
-	int xMax = MIN(refImage->numCols - 1, xCen + size);
-	int yMin = MAX(0, yCen - size);
-	int yMax = MIN(refImage->numRows - 1, yCen + size);
-
-	// Tweak position if necessary
-	if (xCen >= refImage->numCols) {
-	    xCen = refImage->numCols - 1;
-	}
-	if (yCen >= refImage->numRows) {
-	    yCen = refImage->numRows - 1;
-	}
-	if (xCen < 0) {
-	    xCen = 0;
-	}
-	if (yCen < 0) {
-	    yCen = 0;
-	}
-
-	psTrace("phot", 3, "Postage stamp: %d,%d-->%d,%d\n", xMin, yMin, xMax, yMax);
-
-	int numPix = (xMax - xMin + 1) * (yMax - yMin + 1);
-	psArray *pixels = psArrayAlloc(numPix);
-	psVector *fluxes = psVectorAlloc(numPix, PS_TYPE_F32);
-	psVector *errors = psVectorAlloc(numPix, PS_TYPE_F32);
-
-	int index = 0;
-	for (int y = yMin; y <= yMax; y++) {
-	    for (int x = xMin; x <= xMax; x++) {
-		psVector *pixel = psVectorAlloc(2, PS_TYPE_F32);
-		pixel->data.F32[0] = x;
-		pixel->data.F32[1] = y;
-		pixels->data[index] = pixel;
-		fluxes->data.F32[index] = refImage->data.F32[y][x];
-		errors->data.F32[index] = sqrt(refImage->data.F32[y][x]);
-		index++;
-	    }
-	}
-
-	psVector *params = psVectorAlloc(7, PS_TYPE_F32);
-	params->data.F32[0] = refMedian;// Background
-	params->data.F32[1] = refImage->data.F32[yCen][xCen]; // Normalisation
-	params->data.F32[2] = coords->x; // x position
-	params->data.F32[3] = coords->y; // y position
-	params->data.F32[4] = 1.0;	// width in x
-	params->data.F32[5] = 1.0;	// width in y
-	params->data.F32[6] = M_PI/4;	// width in xy
-
-	bool fitOK = psMinimizeLMChi2(min, NULL, params, NULL, pixels, fluxes, NULL, gaussian2d_deriv);
-	//bool fitOK = psMinimizeChi2Powell(min, params, NULL, pixels, fluxes, NULL, gaussian2d_noderiv);
-
-	printf("Final result: %d %d %f --> %f %f %f %f %f %f %f\n", fitOK, min->iter, min->value,
-	       params->data.F32[0], params->data.F32[1], params->data.F32[2], params->data.F32[3],
-	       params->data.F32[4], params->data.F32[5], params->data.F32[6]);
-
-	refFlux->data.F32[i] = params->data.F32[1] * 2.0 * M_PI * params->data.F32[4] * params->data.F32[5];
-	refBG->data.F32[i] = params->data.F32[0];
-	
-	printf("%f %f --> %f %f %f %f\n", coords->x, coords->y, params->data.F32[2], params->data.F32[3],
-	       refFlux->data.F32[i], refBG->data.F32[i]);
-
-    }
-
-    for (int i = 0; i < num; i++) {
-	psPlane *coords = starList->data[i];
-	printf("%f %f --> %f %f\n", coords->x, coords->y, -2.5*log10f(refFlux->data.F32[i]) + ZP,
-	       refBG->data.F32[i]);
-    }
-
-}