Index: /branches/ccl_branches/pstest/Makefile
===================================================================
--- /branches/ccl_branches/pstest/Makefile	(revision 41121)
+++ /branches/ccl_branches/pstest/Makefile	(revision 41121)
@@ -0,0 +1,59 @@
+default: tests
+help:
+	@echo "USAGE: make pstest"
+
+CC      =       gcc
+SRC     =       src
+BIN     =       bin
+
+DESTBIN =       $(HOME)/src/bin/$(ARCH)
+
+LPSLIB  :=      $(shell pslib-config --libs)
+IPSLIB  :=      $(shell pslib-config --cflags)
+
+INCS	= 	$(IPSLIB)
+LIBS	= 	$(LPSLIB)
+CFLAGS	=	$(INCS) -std=c99 -Wall -Werror -g
+LFLAGS	=	$(LIBS) 
+
+TESTS = \
+$(BIN)/tst_psPolynomial.$(ARCH)
+
+tests: $(TESTS)
+
+# dependancy rules for binary code #########################
+.PRECIOUS: %.$(ARCH).o
+.PRECIOUS: $(BIN)/%.$(ARCH)
+
+%.$(ARCH).o : %.c
+	$(CC) $(CFLAGS) -c $< -o $@
+
+$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
+	$(CC) $^ -o $@ $(LFLAGS)
+	@echo "done with $@"
+
+$(DESTBIN)/%: $(BIN)/%.$(ARCH)
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
+	rm -f $(DESTBIN)/$*
+	cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$*
+
+$(INSTALL): % : $(BIN)/%.$(ARCH)
+
+%.clean :
+	rm -f $(BIN)/$*.$(ARCH)
+	rm -f $(SRC)/*.$(ARCH).o
+
+%.install:
+	make $(DESTBIN)/$*
+
+# utilities #################################################
+
+install:
+	for i in $(INSTALL); do make $$i.install; done
+
+clean:	
+	rm -f $(BIN)/*.$(ARCH)
+	rm -f `find . -name "*.o"`
+	rm -f `find . -name "*~"`
+	rm -f `find . -name "#*"`
Index: /branches/ccl_branches/pstest/src/tst_psPolynomial.c
===================================================================
--- /branches/ccl_branches/pstest/src/tst_psPolynomial.c	(revision 41121)
+++ /branches/ccl_branches/pstest/src/tst_psPolynomial.c	(revision 41121)
@@ -0,0 +1,526 @@
+#include <stdio.h>
+#include <math.h>
+#include "pslib_strict.h"
+#include "psTest.h"
+
+static psS32 testPolynomial1DEval_01(void);
+static psS32 testPolynomial1DEvalVector_01(void);
+static psS32 testPolynomial1DEvalVector_02(void);
+static psS32 testVectorFitPolynomial1D_01(void);
+static psS32 testVectorFitPolynomial1D_02(void);
+static psS32 testVectorClipFitPolynomial1D_01(void);
+
+/* list of test functions in test suite */
+testDescription tests[] = {
+                              {testPolynomial1DEval_01,1,"psPolynomial1DEval_01",0,false},
+                              {testPolynomial1DEvalVector_01,2,"psPolynomial1DEvalVector_01",0,false},
+                              {testPolynomial1DEvalVector_02,2,"psPolynomial1DEvalVector_02",0,false},
+                              {testVectorFitPolynomial1D_01,2,"psVectorFitPolynomial1D_01",0,false},
+                              {testVectorFitPolynomial1D_02,2,"psVectorFitPolynomial1D_02",0,false},
+                              {testVectorClipFitPolynomial1D_01,2,"psVectorFitClipPolynomial1D_01",0,false},
+                              {NULL}
+                          };
+
+/* basic main function: just need to define the test suite name */
+int main(int argc, char* argv[])
+{
+    //psTraceSetLevel("VectorFitPolynomial1DOrd", 6);
+    psLogSetFormat("HLNM");
+    psLogSetLevel(PS_LOG_INFO);
+    return !runTestSuite(stderr,"psPolynomialExamples", tests, argc, argv);
+}
+
+# define X0 +5.0
+# define X1 +3.0
+# define X2 -1.0
+# define X3 -2.0
+
+
+static psPolynomial1D *makePolynomial1D(int Nord, psPolynomialType type, int ordMask) {
+
+    // define the polynomial
+    psPolynomial1D *poly = psPolynomial1DAlloc(Nord, type);
+    switch (Nord) {
+      case 3:
+	poly->coeff[3] = (ordMask==3) ?  0:X3;
+      case 2:
+	poly->coeff[2] = (ordMask==2) ?  0:X2;
+      case 1:
+	poly->coeff[1] = (ordMask==1) ?  0:X1;
+      case 0:
+	poly->coeff[0] = (ordMask==0) ?  0:X0;
+	break;
+      default:
+	break;
+    }
+    return (poly);
+}
+
+static void fillVectors(psVector *x, psVector *y, int Nord, int Npts, psPolynomialType type) {
+
+    psF32 X;//, P1, P2, P3;
+    double c[] = {(double)X0,(double)X1,(double)X2,(double)X3};
+    //double sv,d,dd;
+    double sum, Tn, Tn_1, Tn_2;
+
+    for (int i = 0; i < Npts; i++) {
+	X = (float)i/Npts - 0.5;
+	x->data.F32[i] = X;
+	switch (type) {
+	  case PS_POLYNOMIAL_ORD:
+	    switch (Nord) {
+	      case 0:
+		y->data.F32[i] = c[0];
+		break;
+	      case 1:
+		y->data.F32[i] = c[0] + c[1]*X;
+		break;
+	      case 2:
+		y->data.F32[i] = c[0] + c[1]*X + c[2]*X*X;
+		break;
+	      case 3:
+		y->data.F32[i] = c[0] + c[1]*X + c[2]*X*X + c[3]*X*X*X;
+		break;
+	      default:
+		y->data.F32[i] = 0;
+	    }
+	    break;
+	  case PS_POLYNOMIAL_CHEB:
+	    switch (Nord) {
+	      case 0:
+		y->data.F32[i] = c[0];
+		break;
+	      case 1:
+		y->data.F32[i] =c[0] + c[1]*X;
+		break;
+	      default:
+		//Clenshaw's formula approach (NR 5.5)
+		//d=0.0;
+		//dd=0.0;
+		//for (int j=Nord; j>=1; j--) {
+		//    sv = d;
+		//    d = (double)2.0*(double)X*d - dd + c[j];
+		//    dd = sv;
+		//}
+		//y->data.F32[i] = (double)X*d-dd+(double)0.5*c[0];
+
+		//Straightforward approach
+		Tn_1=X; //T1
+		Tn_2=1; //T0
+		sum=0.;
+		for (int j=2; j<=Nord; j++) {
+		    Tn = 2.0*X*Tn_1 - Tn_2;
+		    sum = sum + c[j]*Tn;
+		    Tn_2 = Tn_1;
+		    Tn_1 = Tn;
+		}
+		y->data.F32[i] = sum + c[1]*X + c[0] - 0.5*c[0];
+	    }
+	    break;
+	}
+    }
+    return;
+}
+
+# define NPTS 10
+
+int testPolynomial1DEval_01(void) {
+    
+    //# define NORD 3
+    int NORD;
+    int maxNORD = 3;
+    int retval = 0;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types
+	for (int j = 1; j <= maxNORD; j++) {  //loop over the polynomial orders
+	    int nGood = 0, nBad = 0; 
+	    NORD = j; 
+	    //printf("Testing NORD=%d\n", NORD);
+	    psPolynomial1D *poly = makePolynomial1D(NORD, type[k], -1); 
+	    
+	    //create x vector, known y vector: 
+	    psVector *x = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    psVector *y = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    fillVectors (x, y, NORD, NPTS, type[k]); 
+	    
+	    psF32 xV, yV; 
+	    
+	    for (int i = 0; i < NPTS; i++) { 
+		xV = x->data.F32[i]; 
+		yV = psPolynomial1DEval (poly, xV); 
+		//Change some values around to make it fail!
+		//if (i>90 && NORD==2) {yV = yV + 0.1;}
+		if (yV == y->data.F32[i]) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    printf("%.3f   -->  %.10f  & %.10f\tdel=%e\n",xV,yV,y->data.F32[i],y->data.F32[i]-yV); 
+		} 
+	    } 
+	    psFree (x); 
+	    psFree (y); 
+	    psFree (poly); 
+	    
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+	}
+    }
+    return retval;
+}
+
+
+int testPolynomial1DEvalVector_01(void) {
+    
+    //# define NORD 3
+    int NORD;
+    int maxNORD = 3;
+    int retval = 0;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+    
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types
+	for (int j = 1; j <= maxNORD; j++) {  //loop over the polynomial orders
+	    int nGood = 0, nBad = 0; 
+	    NORD = j; 
+	    //printf("Testing NORD=%d\n", NORD);
+	    psPolynomial1D *poly = makePolynomial1D(NORD, type[k], -1); 
+	    
+	    // create x vector, known y vector:
+	    psVector *x = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    psVector *y = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    fillVectors (x, y, NORD, NPTS, type[k]); 
+	    
+	    psVector *xV = x; 
+	    psVector *yV = psPolynomial1DEvalVector(poly, xV); 
+	    
+	    for (int i = 0; i < NPTS; i++) { 
+		//Change some values around to make it fail!
+		//if (i>90 && NORD==2) {yV->data.F32[i] = yV->data.F32[i] + 0.1;}
+		if (yV->data.F32[i] == y->data.F32[i]) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    printf("%.3f   -->  %.3f  & %.3f\tdel=%e\n",xV->data.F32[i],yV->data.F32[i],y->data.F32[i],y->data.F32[i]-yV->data.F32[i]); 
+		} 
+	    } 
+	    psFree (x); 
+	    psFree (y); 
+	    psFree (yV); 
+	    psFree (poly); 
+	    
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+	}
+    }
+    return retval;
+}
+int testPolynomial1DEvalVector_02(void) {
+    //  check that psPolynomial1DEvalVector obeys the psPolynomial1D.mask 
+    int NORD;
+    int maxNORD = 3;
+    int retval = 0;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+    
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types
+	for (int imask=0; imask<=maxNORD; imask++) {  // mask one order at a time
+	    int j= maxNORD;
+	    NORD = j;
+	    int nGood = 0, nBad = 0; 
+	    psPolynomial1D *poly       = makePolynomial1D(NORD, type[k], -1);
+	    psPolynomial1D *polynomask = makePolynomial1D(NORD, type[k], imask);
+	    poly->mask[imask]=1;
+
+	    // create x vector, known y vector:
+	    psVector *x = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    psVector *y = psVectorAlloc (NPTS, PS_TYPE_F32); 
+	    fillVectors (x, y, NORD, NPTS, type[k]); 
+	    
+	    psVector *xV = x; 
+	    psVector *yV = psPolynomial1DEvalVector(poly, xV); 
+	    psVector *yVnm = psPolynomial1DEvalVector(polynomask, xV); 
+	    
+	    for (int i = 0; i < NPTS; i++) { 
+		if (yV->data.F32[i] == yVnm->data.F32[i]) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    printf("x:%.3f  y:%.3f , %.3f , %.3f\tdel=%e\n",xV->data.F32[i],y->data.F32[i],yV->data.F32[i],yVnm->data.F32[i],yVnm->data.F32[i]-yV->data.F32[i]); 
+		} 
+	    } 
+	    psFree (x); 
+	    psFree (y); 
+	    psFree (yV); 
+	    psFree (yVnm); 
+	    psFree (poly); 
+	    psFree (polynomask); 
+	    
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+	}
+    }
+    return retval;
+}
+
+
+
+int testVectorFitPolynomial1D_01(void) { 
+    int retval=0;
+    int NORD, maxNORD=3;
+    float tolerance = 1.0E-5;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+    
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types   XXX: not using Chebyshevs right now!!!
+	for (int j = 1; j <= maxNORD; j++) {  //loop over the polynomial orders 
+	    int nGood = 0, nBad = 0; 
+	    NORD = j;
+	    printf("Testing NORD=%d\n", NORD);
+	    psPolynomial1D *poly = makePolynomial1D(NORD, type[k], -1);
+	    psPolynomial1D *polycopy = makePolynomial1D(NORD, type[k], -1);
+	    psVector *f = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *fErr = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *x = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    
+	    fillVectors(x, f, NORD, NPTS, type[k]);
+	    psTime *now = psTimeGetNow(PS_TIME_UTC);
+	    const psRandom *r = psRandomAllocSpecific(PS_RANDOM_TAUS, now->sec);
+	    for (int i = 0; i < NPTS; i++) {
+		double errgen = psRandomGaussian(r);
+		double amp = 1.0E-3;
+		fErr->data.F32[i] = amp*errgen;
+		//fErr->data.F32[i] = 0.;
+	    }
+	    psFree(now);
+	    psFree(r);
+
+	    //f->data.F32[NPTS-1] = X0 + 10.*X1*x->data.F32[NPTS-1];
+	    //fErr->data.F32[NPTS-1] = 0.0000001;
+
+
+	    psVectorFitPolynomial1D(poly, NULL, 0, f, fErr, x);
+	    for (int i = 0; i <= NORD; i++) { 
+		if (poly->coeff[i] <= polycopy->coeff[i]+tolerance 
+		      &&
+		    poly->coeff[i] >= polycopy->coeff[i]-tolerance) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    printf("polycopy->coeff[%d]=%.20f    poly->coeff[%d]=%.20f\n",i,polycopy->coeff[i],i,poly->coeff[i]); 
+		}
+		if (polycopy->coeffErr[i] != poly->coeffErr[i]) { 
+		    printf("polycopy->coeffErr[%d]=%.20f    poly->coeffErr[%d]=%.20f\n",i,polycopy->coeffErr[i],i,poly->coeffErr[i]); 
+		}		
+	    } 
+
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+
+	    printf("NORD: %d   retval: %d\n",NORD,retval);
+	    psFree(poly);
+	    psFree(polycopy);
+	    psFree(f);
+	    psFree(fErr);
+	    psFree(x);	    
+	}
+    }
+    return retval;
+}
+
+
+int testVectorFitPolynomial1D_02(void) { 
+    // this time test that the masking in psPolynomial1D works correctly
+    int retval=0;
+    int NORD, maxNORD=3;
+    float tolerance = 1.0E-5;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+    
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types   XXX: not using Chebyshevs right now!!!
+	for (int imask=0; imask<=maxNORD; imask++) {  // mask one order at a time
+	    int j= maxNORD;
+	    NORD = j;
+	    int nGood = 0, nBad = 0; 
+	    psPolynomial1D *poly       = makePolynomial1D(NORD, type[k], -1);
+	    psPolynomial1D *polynomask = makePolynomial1D(NORD, type[k], imask);
+	    poly->mask[imask]=1;
+	    //poly and polynomask look the same now.
+
+	    psVector *mask = psVectorAlloc(NPTS, PS_TYPE_U8);
+	    psVector *f = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *fErr = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *x = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    
+	    fillVectors(x, f, NORD, NPTS, type[k]);
+	    
+	    // add random errors to the vector
+	    psTime *now = psTimeGetNow(PS_TIME_UTC);
+	    const psRandom *r = psRandomAllocSpecific(PS_RANDOM_TAUS, now->sec);
+	    for (int i = 0; i < NPTS; i++) {
+		double errgen = psRandomGaussian(r);
+		double amp = 1.0E-3;
+		fErr->data.F32[i] = amp*errgen;
+		//fErr->data.F32[i] = 0.;
+	    }
+	    psFree(now);
+	    psFree(r);
+	    
+	    
+	    //mask->data.U8[imask] = 0xff;
+	    //for (int ii=0; ii<=mask->n; ii++) {
+	    //printf("  %d",mask->data.U8[ii]);
+	    //}
+	    //printf("\n");
+
+	    poly->mask[imask] = 1;
+	    for (int ii=0; ii<=poly->nX; ii++) {
+		printf("%x  ", poly->mask[ii]);
+	    }
+	    printf("\n");
+
+	    psVectorFitPolynomial1D(poly, NULL, 0xff, f, fErr, x);
+	    psVectorFitPolynomial1D(polynomask, NULL, 0xff, f, fErr, x);
+
+	    psVector *xx = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *yy = psPolynomial1DEvalVector(poly, xx);
+	    psVector *del;
+	    psBinaryOp(del, f, "-", yy);
+	    psPolynomial1D *polydel = makePolynomial1D(NORD, type[k], -1);
+	    psVectorFitPolynomial1D(polydel, NULL,0xff, f, fErr, x);
+	    
+
+
+	    //mask->data.U8[imask] = 0;
+	    poly->mask[imask] = 0;
+	    for (int i = 0; i <= NORD; i++) { 
+		printf("polynomask->coeff[%d]=%.8f    poly->coeff[%d]=%.8f    polydel->coeff[%d]=%.8f\n",i,polynomask->coeff[i],i,poly->coeff[i],i,polydel->coeff[i]); 
+		if (poly->coeff[i] <= polynomask->coeff[i]+tolerance 
+		    &&
+		    poly->coeff[i] >= polynomask->coeff[i]-tolerance) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    //printf("polynomask->coeff[%d]=%.8f    poly->coeff[%d]=%.8f\n",i,polynomask->coeff[i],i,poly->coeff[i]); 
+		}
+		if (polynomask->coeffErr[i] != poly->coeffErr[i]) { 
+		    //printf("polynomask->coeffErr[%d]=%.8f    poly->coeffErr[%d]=%.8f\n",i,polynomask->coeffErr[i],i,poly->coeffErr[i]); 
+		}		
+	    }
+	    
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+	    
+	    printf("Mask Ord: %d   retval: %d\n", imask, retval);
+	    psFree(poly);
+	    psFree(polynomask);
+	    psFree(mask);
+	    psFree(f);
+	    psFree(fErr);
+	    psFree(x);	    
+	    psFree(xx);
+	    psFree(yy);
+	    psFree(del);
+	    psFree(polydel);
+
+	}
+    }
+    return retval;
+}
+
+
+
+
+
+int testVectorClipFitPolynomial1D_01(void) { 
+    int retval=0;
+    int NORD, maxNORD=3;
+    double Nsigma = 3.;
+    int Niter= 3.;
+    float tolerance = 1.0E-5;
+    psPolynomialType type[] = {PS_POLYNOMIAL_ORD, PS_POLYNOMIAL_CHEB};
+    
+    for (int k = 0; k<1; k++) {  //loop over the 2 polynomial types  XXX: psVectorClipFirPolynomial1D does not
+	                          //                                  work with Chebyshev now
+	for (int j = 1; j <= maxNORD; j++) {  //loop over the polynomial orders
+	    int nGood = 0, nBad = 0; 
+	    NORD = j;
+	    printf("Testing NORD=%d\n", NORD);
+	    psPolynomial1D *poly = makePolynomial1D(NORD, type[k], -1);
+	    psPolynomial1D *polycopy = makePolynomial1D(NORD, type[k], -1);
+	    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	    stats->clipSigma = Nsigma;
+	    stats->clipIter  = Niter;
+	    psVector *f = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *fErr = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    psVector *x = psVectorAlloc(NPTS, PS_TYPE_F32);
+	    
+	    fillVectors(x, f, NORD, NPTS, type[k]);
+	    psTime *now = psTimeGetNow(PS_TIME_UTC);
+	    const psRandom *r = psRandomAllocSpecific(PS_RANDOM_TAUS, now->sec);
+	    for (int i = 0; i < NPTS; i++) {
+		double errgen = psRandomGaussian(r);
+		double amp = 1.0E-3;
+		fErr->data.F32[i] = amp*errgen;
+		//fErr->data.F32[i] = 0.;
+	    }
+	    psFree(now);
+	    psFree(r);
+
+	    //f->data.F32[NPTS-1] = X0 + 10.*X1*x->data.F32[NPTS-1];
+	    //fErr->data.F32[NPTS-1] = 0.0000001;
+
+
+	    psVectorClipFitPolynomial1D(poly, stats, NULL, 2, f, fErr, x);
+	    for (int i = 0; i <= NORD; i++) { 
+		if (poly->coeff[i] <= polycopy->coeff[i]+tolerance 
+		      &&
+		    poly->coeff[i] >= polycopy->coeff[i]-tolerance) { 
+		    nGood++; 
+		} else { 
+		    nBad++; 
+		    printf("polycopy->coeff[%d]=%.20f    poly->coeff[%d]=%.20f\n",i,polycopy->coeff[i],i,poly->coeff[i]); 
+		}
+		if (polycopy->coeffErr[i] != poly->coeffErr[i]) { 
+		    printf("polycopy->coeffErr[%d]=%.20f    poly->coeffErr[%d]=%.20f\n",i,polycopy->coeffErr[i],i,poly->coeffErr[i]); 
+		}		
+	    } 
+
+	    if (nBad > 0) { 
+		psError(PS_ERR_UNKNOWN,true,"%d incorrect polynomial evaluations for NORD=%d\n", nBad, NORD); 
+		retval = retval + pow(2,NORD); 
+	    } 
+
+	    printf("NORD: %d   retval: %d\n",NORD,retval);
+	    psFree(poly);
+	    psFree(polycopy);
+	    psFree(stats);
+	    psFree(f);
+	    psFree(fErr);
+	    psFree(x);	    
+	}
+    }
+    return retval;
+}
+
+
+
+
+
+
+
+//test psStats
+
+//test psGaussian
+
+//test psMinimization (and psMinimizeLMChi2 and psMinimizeChi2Powell?)
+
+//test psRandomUniform / psRandomGausian / psRandomPoisson
+
+//test psSpline?
Index: /branches/ccl_branches/simtest/Makefile.am
===================================================================
--- /branches/ccl_branches/simtest/Makefile.am	(revision 41121)
+++ /branches/ccl_branches/simtest/Makefile.am	(revision 41121)
@@ -0,0 +1,15 @@
+
+install_files = \
+	sim.pro \
+	images.pro
+
+installdir = $(datadir)/pantasks/modules
+
+install_DATA = $(install_files)
+
+install-data-hook:
+	chmod 0755 $(installdir) 
+
+EXTRA_DIST = $(install_files)
+
+ACLOCAL_AMFLAGS = -I m4
Index: /branches/ccl_branches/simtest/autogen.sh
===================================================================
--- /branches/ccl_branches/simtest/autogen.sh	(revision 41121)
+++ /branches/ccl_branches/simtest/autogen.sh	(revision 41121)
@@ -0,0 +1,103 @@
+#!/bin/sh
+# Run this to generate all the initial makefiles, etc.
+
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+ORIGDIR=`pwd`
+cd $srcdir
+
+PROJECT=simtest
+TEST_TYPE=-f
+# change this to be a unique filename in the top level dir
+FILE=autogen.sh
+
+DIE=0
+
+LIBTOOLIZE=libtoolize
+ACLOCAL="aclocal $ACLOCAL_FLAGS"
+AUTOHEADER=autoheader
+AUTOMAKE=automake
+AUTOCONF=autoconf
+
+#($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $LIBTOOLIZE installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/"
+#        DIE=1
+#}
+
+($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $ACLOCAL installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+#($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $AUTOHEADER installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+#        DIE=1
+#}
+
+($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOMAKE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOCONF installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+test $TEST_TYPE $FILE || {
+        echo "You must run this script in the top-level $PROJECT directory"
+        exit 1
+}
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+#$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed"
+$ACLOCAL || echo "$ACLOCAL failed"
+#$AUTOHEADER || echo "$AUTOHEADER failed"
+$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed"
+$AUTOCONF || echo "$AUTOCONF failed"
+
+cd $ORIGDIR
+
+run_configure=true
+for arg in $*; do
+    case $arg in
+        --no-configure)
+            run_configure=false
+            ;;
+        *)
+            ;;
+    esac
+done
+
+if $run_configure; then
+    $srcdir/configure --enable-maintainer-mode "$@"
+    echo
+    echo "Now type 'make' to compile $PROJECT."
+else
+    echo
+    echo "Now run 'configure' and 'make' to compile $PROJECT."
+fi
Index: /branches/ccl_branches/simtest/configure.ac
===================================================================
--- /branches/ccl_branches/simtest/configure.ac	(revision 41121)
+++ /branches/ccl_branches/simtest/configure.ac	(revision 41121)
@@ -0,0 +1,17 @@
+AC_PREREQ(2.59)
+
+AC_INIT([simtest], [0.9.0], [ipp-support@ifa.hawaii.edu])
+AC_CONFIG_SRCDIR([sim.pro])
+
+AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2])
+AM_MAINTAINER_MODE
+
+AC_PROG_INSTALL
+
+MODULE_DIR="${datadir}/pantasks/modules"
+AC_SUBST(MODULE_DIR)
+
+AC_CONFIG_FILES([
+  Makefile
+])
+AC_OUTPUT
Index: /branches/ccl_branches/simtest/images.pro
===================================================================
--- /branches/ccl_branches/simtest/images.pro	(revision 41121)
+++ /branches/ccl_branches/simtest/images.pro	(revision 41121)
@@ -0,0 +1,210 @@
+
+# set global parameters for the simulations
+macro mkimages.init
+ if ($0 != 7)
+  echo "USAGE: mkimage.empty (Nx) (Ny) (RDNOISE) (SKY) (Nstars) (psfSigma)"
+  break
+ end
+
+ $mkiNx = $1
+ $mkiNy = $2
+ $mkiRD = $3
+ $mkiSKY = $4
+ $mkiPSFsigma = $6
+
+ gaussdev mki_noisev {$mkiNx*$mkiNy} 0.0 1.0
+ dimenup mki_noisev mki_noisei $mkiNx $mkiNy
+ mcreate blank $mkiNx $mkiNy
+end
+
+# generate a fake image with bias and sky, no stars
+macro mkimage.empty
+ if ($0 != 2)
+  echo "USAGE: mkimage.empty (buffer)"
+  break
+ end
+
+ set $1 = blank + $mkiRD*mki_noisei + sqrt($mkiSKY)*mki_noisei + $mkiSKY
+end
+
+# set global parameters for the simulations
+macro mkimages.stars
+ if ($0 != 7)
+  echo "USAGE: mkimage.empty (Nx) (Ny) (RDNOISE) (SKY) (Nstars) (psfSigma)"
+  break
+ end
+
+ $mkiNx = $1
+ $mkiNy = $2
+ $mkiRD = $3
+ $mkiSKY = $4
+ $mkiPSFsigma = $6
+
+ gaussdev mki_noisev {$mkiNx*$mkiNy} 0.0 1.0
+ dimenup mki_noisev mki_noisei $mkiNx $mkiNy
+ mcreate blank $mkiNx $mkiNy
+
+ # mkstars $5 0.6
+ mkstars $5 0.4
+ # mkstars $5 0.3
+end
+
+# generate a fake image with bias, sky, stars, gaussian PFS
+macro mkimage.stars.gauss
+ if ($0 != 2)
+  echo "USAGE: mkimage.stars (buffer)"
+  break
+ end
+
+ set _tmp.stars = mki_stars_image
+
+ mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkgauss
+
+ # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
+ set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
+end
+
+# generate a fake image with bias, sky, stars, gaussian PFS
+macro mkimage.stars.power
+ if ($0 != 2)
+  echo "USAGE: mkimage.stars (buffer)"
+  break
+ end
+
+ set _tmp.stars = mki_stars_image
+
+ mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkpower
+
+ # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
+ set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
+end
+
+# generate a fake image with bias, sky, stars, gaussian PFS
+macro mkimage.stars.wing
+ if ($0 != 2)
+  echo "USAGE: mkimage.stars (buffer)"
+  break
+ end
+
+ set _tmp.stars = mki_stars_image
+
+ mkimage.smooth _tmp.stars _tmp.smooth $mkiPSFsigma mkwing
+
+ # noise sigma = sqrt (_tmp.smooth + $mkiRD^2 + $mkiSKY)
+ set $1 = _tmp.smooth + $mkiSKY + sqrt(_tmp.smooth + $mkiSKY + $mkiRD^2)*mki_noisei
+end
+
+# generate the vectors of star coordinates, fluxes and magnitudes
+macro mkstars
+ if ($0 != 3)
+   echo "USAGE: mkstars (Nstars) (alpha)"
+   break
+ end
+
+ local sigma alpha Nstars
+ $Nstars = $1
+ $alpha = $2
+ # alpha is dlogN/dmag = 0.6 for N ~ S^-3/2
+
+ # $sigma is flux of star with S/N = 1 (sigma_sky * 4pi sigma_psf^2)
+ $sigma = sqrt($mkiRD^2 + $mkiSKY) * 4*3.1416*$mkiPSFsigma^2
+ 
+ # faintest magnitudes:
+ # include stars down to 1 sigma
+ $M1 = -2.5*log($sigma)
+
+ # brightest magnitude in sim
+ $M0 = $M1 - log($alpha*$Nstars/0.43)/$alpha
+
+ create _tmp 0 $Nstars
+ set mki_stars_x = $mkiNx * rnd(_tmp)
+ set mki_stars_y = $mkiNy * rnd(_tmp)
+
+ # random variable between 0 and Nstars
+ set _tmpvar = $Nstars * rnd(_tmp)
+
+ # random mags drawn from random variable
+ set mki_stars_mag = $M0 + log(_tmpvar*$alpha/0.43)/$alpha
+
+ set mki_stars_flux = ten(-0.4*mki_stars_mag)
+ sort mki_stars_mag mki_stars_flux mki_stars_x mki_stars_y
+
+ # insert the stars
+ mcreate mki_stars_image $mkiNx $mkiNy
+ for i 0 $Nstars
+   zap mki_stars_image mki_stars_x[$i] mki_stars_y[$i] 1 1 -v mki_stars_flux[$i]
+ end
+end
+
+####
+macro mkimage.smooth
+ if ($0 != 5)
+   echo "USAGE: smooth (in) (out) (sigma) (func)"
+   break
+ end
+ fft2d $1 0 to sm_I_r sm_I_i
+ set sm_g = zero($1)
+ keyword sm_g NAXIS1 SMnx
+ keyword sm_g NAXIS2 SMny
+ $4 sm_g $3 -c 0 0
+ $4 sm_g $3 -c $SMnx 0
+ $4 sm_g $3 -c 0 $SMny
+ $4 sm_g $3 -c $SMnx $SMny
+ stats -q sm_g - - - -
+ set sm_g = sm_g/$TOTAL
+ fft2d sm_g 0 to sm_G_r sm_G_i
+ set sm_r = sm_G_r*sm_I_r - sm_G_i*sm_I_i
+ set sm_i = sm_G_r*sm_I_i + sm_G_i*sm_I_r
+ fft2d -inverse sm_r sm_i to sm_I_r sm_I_i
+
+ set $2 = sm_I_r*$SMnx*$SMny*$SMnx*$SMny
+
+ if (1)
+  delete sm_I_r 
+  delete sm_I_i
+  delete sm_G_r 
+  delete sm_G_i
+  delete sm_r
+  delete sm_i
+  delete sm_g
+ end
+end
+
+macro mkpower
+ if ($0 != 6)
+  echo "USAGE: mkpower (buffer) (sigma) -c X Y"
+  break
+ end
+
+ delete -q tmp
+
+ keyword $1 NAXIS1 tmpnx
+ keyword $1 NAXIS2 tmpny
+ mcreate tmp $tmpnx $tmpny
+ set tmpx = xramp(tmp) - $4
+ set tmpy = yramp(tmp) - $5
+
+# $sigma_y = $2 * $AxisRatio
+# $sigma_xy = dsin(2*$Theta)*(
+
+ set tmpz = tmpx^2 / (2 * $2^2) + tmpy^2 / (2 * $2^2)
+ set $1 = $1 + 1 / (1 + tmpz + tmpz^2.5)
+end
+
+macro mkwing
+ if ($0 != 6)
+  echo "USAGE: mkpower (buffer) (sigma) -c X Y"
+  break
+ end
+
+ delete -q tmp
+
+ keyword $1 NAXIS1 tmpnx
+ keyword $1 NAXIS2 tmpny
+ mcreate tmp $tmpnx $tmpny
+ set tmpx = xramp(tmp) - $4
+ set tmpy = yramp(tmp) - $5
+
+ set tmpz = tmpx^2 / (2 * $2^2) + tmpy^2 / (2 * $2^2)
+ set $1 = $1 + 1 / (1 + tmpz + tmpz^2)
+end
Index: /branches/ccl_branches/simtest/sim.pro
===================================================================
--- /branches/ccl_branches/simtest/sim.pro	(revision 41121)
+++ /branches/ccl_branches/simtest/sim.pro	(revision 41121)
@@ -0,0 +1,262 @@
+### Generate simulated test data for testing the Image Processing Pipeline
+
+module images.pro
+
+$flatlevel = 1000
+$biaslevel = 500
+$skyrate  = 200
+$darkrate = 100
+$shutter = 0.5
+
+macro cleanup
+  exec rm -f *.fit
+  exec rm -f *.jpg
+  exec rm -f *.fits
+  exec rm -rf logs
+end
+
+$simdir = `ipp_datapath.pl path://SIMTEST`
+
+macro init
+  exec pxadmin -delete -dbname simtest
+  exec pxadmin -create -dbname simtest
+  exec mkdir -p $simdir/raw
+  exec mkdir -p $simdir/proc
+  mkbiases
+  mkdarks
+  mkshutter
+  mkflats
+  mkobjects
+  exec ipp_serial_inject.pl $simdir/raw/*.fits --workdir path://SIMTEST/proc --dbname simtest
+end
+
+### create a bias detrun:
+# exec dettool -pretend -definebyquery -inst SIMTEST -filter NONE -det_type bias -select_exp_type bias -workdir path://SIMTEST/mkbias.0 -dbname simtest
+
+### If something goes wrong:
+# exec dettool -updatedetrun -det_id 1 -state stop -dbname simtest
+
+### create a dark detrun:
+# exec dettool -pretend -definebyquery -inst SIMTEST -filter NONE -det_type dark -select_exp_type dark -workdir path://SIMTEST/mkdark.0 -dbname simtest
+
+### create a shutter detrun:
+# exec dettool -pretend -definebyquery -inst SIMTEST -filter NONE -det_type shutter -select_exp_type shutter -workdir path://SIMTEST/mkshutter.0 -dbname simtest
+
+### create a flat detrun:
+# exec dettool -pretend -definebyquery -inst SIMTEST -filter B -det_type flat -select_exp_type flat -select_filter B -workdir path://SIMTEST/mkflat.0 -dbname simtest 
+# exec dettool -pretend -definebyquery -inst SIMTEST -filter V -det_type flat -select_exp_type flat -select_filter V -workdir path://SIMTEST/mkflat.1 -dbname simtest 
+
+# generate a set of fake biases
+macro mkbiases
+  echo "generating 10 bias images"
+  for i 0 10
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set bias = blank + $mkiRD*mki_noisei + $biaslevel
+    keyword bias AIRMASS  -wf {1 + 0.1*$i}
+    keyword bias FILTER   -w NONE
+    keyword bias POSANGLE -wf 0.0
+    keyword bias RA       -wf 0.0
+    keyword bias DEC      -wf 0.0
+    keyword bias EXPTIME  -wf 0.0
+    keyword bias DARKTIME -wf 0.0
+    keyword bias DATE-OBS -w "2006/10/01"
+    sprintf time "00:00:%02d" $i
+    keyword bias UTC-OBS  -w $time
+    keyword bias CCDSUM   -w "1 1"
+    keyword bias TELESCOP -w simscope
+    keyword bias INSTRUME -w SIMTEST
+    keyword bias OBSTYPE  -w BIAS
+    keyword bias OBJECT   -w BIAS
+    sprintf name "$simdir/raw/bias.%02d.fits" $i
+    keyword bias GAIN     -wf 1.0
+    keyword bias RDNOISE  -wf $mkiRD
+    wd bias $name
+  end
+end
+
+macro mkdarks
+  echo "generating 20 dark images"
+  for i 0 20
+    $exptime = 3 * ($i + 1)
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set dark = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime)*mki_noisei + $darkrate*$exptime + $biaslevel
+    keyword dark AIRMASS  -wf {1 + 0.1*$i}
+    keyword dark FILTER   -w  NONE
+    keyword dark POSANGLE -wf 0.0
+    keyword dark RA       -wf 0.0
+    keyword dark DEC      -wf 0.0
+    keyword dark EXPTIME  -wf $exptime
+    keyword dark DARKTIME -wf $exptime
+    keyword dark DATE-OBS -w "2006/10/01"
+    sprintf time "00:01:%02d" $i
+    keyword dark UTC-OBS  -w $time
+    keyword dark CCDSUM   -w "1 1"
+    keyword dark TELESCOP -w simscope
+    keyword dark INSTRUME -w SIMTEST
+    keyword dark OBSTYPE  -w DARK
+    keyword dark OBJECT   -w DARK
+    keyword dark GAIN     -wf 1.0
+    keyword dark RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/dark.%02d.fits" $i
+    wd dark $name
+  end  
+end
+
+macro mkshutter
+  echo "generating 10 shutter images"
+  for i 0 10
+    $exptime = 3 * ($i + 1)
+    $opentime = $exptime - $shutter
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set shutter = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime + $flatlevel*$opentime)*mki_noisei + $flatlevel*$opentime + $darkrate*$exptime + $biaslevel
+    keyword shutter AIRMASS  -wf 1.0
+    keyword shutter FILTER   -w  NONE
+    keyword shutter POSANGLE -wf 0.0
+    keyword shutter RA       -wf 0.0
+    keyword shutter DEC      -wf 0.0
+    keyword shutter EXPTIME  -wf $exptime
+    keyword shutter DARKTIME -wf $exptime
+    keyword shutter DATE-OBS -w "2006/10/01"
+    sprintf time "00:01:%02d" $i
+    keyword shutter UTC-OBS  -w $time
+    keyword shutter CCDSUM   -w "1 1"
+    keyword shutter TELESCOP -w simscope
+    keyword shutter INSTRUME -w SIMTEST
+    keyword shutter OBSTYPE  -w SHUTTER
+    keyword shutter OBJECT   -w SHUTTER
+    keyword shutter GAIN     -wf 1.0
+    keyword shutter RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/shutter.%02d.fits" $i
+    wd shutter $name
+  end
+end
+
+macro mkflats
+  echo "generating 20 flat images in 2 filters"
+  $exptime = 10.0
+  $opentime = $exptime - $shutter
+  for i 0 10
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set flat = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime + $flatlevel*$opentime)*mki_noisei + $darkrate*$exptime + $flatlevel*$opentime + $biaslevel
+    keyword flat AIRMASS  -wf {1 + 0.1*$i}
+    keyword flat FILTER   -w  V
+    keyword flat POSANGLE -wf 0.0
+    keyword flat RA       -wf 0.0
+    keyword flat DEC      -wf 0.0
+    keyword flat EXPTIME  -wf $exptime
+    keyword flat DARKTIME -wf $exptime
+    keyword flat DATE-OBS -w "2006/10/01"
+    sprintf time "00:01:%02d" $i
+    keyword flat UTC-OBS  -w $time
+    keyword flat CCDSUM   -w "1 1"
+    keyword flat TELESCOP -w simscope
+    keyword flat INSTRUME -w SIMTEST
+    keyword flat OBSTYPE  -w FLAT
+    keyword flat OBJECT   -w FLAT
+    keyword flat GAIN     -wf 1.0
+    keyword flat RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/flat.%02d.fits" $i
+    wd flat $name
+  end  
+
+  for i 10 20
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set flat = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime + $flatlevel*$opentime)*mki_noisei + $darkrate*$exptime + $flatlevel*$opentime + $biaslevel
+    keyword flat AIRMASS  -wf {1 + 0.1*$i}
+    keyword flat FILTER   -w  B
+    keyword flat POSANGLE -wf 0.0
+    keyword flat RA       -wf 0.0
+    keyword flat DEC      -wf 0.0
+    keyword flat EXPTIME  -wf $exptime
+    keyword flat DARKTIME -wf $exptime
+    keyword flat DATE-OBS -w "2006/10/01"
+    sprintf time "00:02:%02d" $i
+    keyword flat UTC-OBS  -w $time
+    keyword flat CCDSUM   -w "1 1"
+    keyword flat TELESCOP -w simscope
+    keyword flat INSTRUME -w SIMTEST
+    keyword flat OBSTYPE  -w FLAT
+    keyword flat OBJECT   -w FLAT
+    keyword flat GAIN     -wf 1.0
+    keyword flat RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/flat.%02d.fits" $i
+    wd flat $name
+  end  
+end
+
+macro mkobjects
+  echo "generating 10 object images in 2 filters"
+  for i 0 5
+    $exptime = 3 * ($i + 1)
+    $opentime = $exptime - $shutter
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set obj = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime + $skyrate*$opentime)*mki_noisei + $darkrate*$exptime + $skyrate*$opentime + $biaslevel
+    keyword obj AIRMASS  -wf {1 + 0.1*$i}
+    keyword obj FILTER   -w  V
+    keyword obj POSANGLE -wf 0.0
+    keyword obj RA       -wf 0.0
+    keyword obj DEC      -wf 0.0
+    keyword obj EXPTIME  -wf $exptime
+    keyword obj DARKTIME -wf $exptime
+    keyword obj DATE-OBS -w "2006/10/01"
+    sprintf time "00:01:%02d" $i
+    keyword obj UTC-OBS  -w $time
+    keyword obj CCDSUM   -w "1 1"
+    keyword obj TELESCOP -w simscope
+    keyword obj INSTRUME -w SIMTEST
+    keyword obj OBSTYPE  -w OBJECT
+    keyword obj OBJECT   -w OBJECT
+    keyword obj GAIN     -wf 1.0
+    keyword obj RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/obj.%02d.fits" $i
+    wd obj $name
+  end  
+
+  for i 5 10
+    $exptime = 3 * ($i - 4)
+    $opentime = $exptime - $shutter
+    mkimages.init 1024 1024 5 1000 0 2.0
+    set obj = blank + sqrt($mkiRD*$mkiRD + $darkrate*$exptime + $skyrate*$opentime)*mki_noisei + $darkrate*$exptime + $skyrate*$opentime + $biaslevel
+    keyword obj AIRMASS  -wf {1 + 0.1*$i}
+    keyword obj FILTER   -w  B
+    keyword obj POSANGLE -wf 0.0
+    keyword obj RA       -wf 0.0
+    keyword obj DEC      -wf 0.0
+    keyword obj EXPTIME  -wf $exptime
+    keyword obj DARKTIME -wf $exptime
+    keyword obj DATE-OBS -w "2006/10/01"
+    sprintf time "00:02:%02d" $i
+    keyword obj UTC-OBS  -w $time
+    keyword obj CCDSUM   -w "1 1"
+    keyword obj TELESCOP -w simscope
+    keyword obj INSTRUME -w SIMTEST
+    keyword obj OBSTYPE  -w OBJECT
+    keyword obj OBJECT   -w OBJECT
+    keyword obj GAIN     -wf 1.0
+    keyword obj RDNOISE  -wf $mkiRD
+    sprintf name "$simdir/raw/obj.%02d.fits" $i
+    wd obj $name
+  end  
+end
+
+macro junk
+-definebyquery Optional arguments, with default values:
+    -det_type          ()                         define the type of detrend run (required)
+    -mode              (master)                   define the mode of this detrend run
+    -exp_type          ()                         search for exp_type
+    -inst              ()                         search for camera
+    -telescope         ()                         search for telescope
+    -filter            ()                         search for filter
+    -uri               ()                         search for uri
+    -set_exp_type      ()                         define exposure type
+    -set_filter        ()                         define filter
+    -set_airmass       (nan)            define airmass
+    -set_exp_time      (nan)            define exposure time
+    -set_ccd_temp      (nan)            define ccd tempature
+    -set_posang        (nan)            define rotator position angle
+    -set_object        ()                         define exposure object
+    -set_registered    (2006-12-13T22:40:18.5)    time detrend run was registered
+    -set_use_begin     ()                         start of detrend run applicable peroid
+    -set_use_end       ()                         end of detrend run applicable peroid
+    -pretend           (FALSE)                    print the exposures that would be included in the detrend run and exit
+end