Index: /trunk/psastro/Makefile
===================================================================
--- /trunk/psastro/Makefile	(revision 5360)
+++ /trunk/psastro/Makefile	(revision 5360)
@@ -0,0 +1,115 @@
+default: psphot
+help:
+	@echo "USAGE: make psphot"
+
+CC      =       gcc
+SRC     =       src
+BIN     =       bin
+
+DESTBIN =       $(HOME)/src/bin/$(ARCH)
+
+LPSLIB  :=      $(shell pslib-config --libs)
+IPSLIB  :=      $(shell pslib-config --cflags)
+
+INCS	= 	$(IPSLIB)
+# LIBS	= 	-lpsmodule $(LPSLIB)
+LIBS	= 	$(LPSLIB)
+CFLAGS	=	$(INCS) -std=c99 -Wall -Werror -g
+# CFLAGS	=	$(INCS) -std=c99 -g
+LFLAGS	=	$(LIBS) 
+
+PSPHOT = \
+$(SRC)/psphot.$(ARCH).o            \
+$(SRC)/psphotArguments.$(ARCH).o   \
+$(SRC)/psphotSetup.$(ARCH).o	   \
+$(SRC)/psphotImageStats.$(ARCH).o  \
+$(SRC)/psphotSourceStats.$(ARCH).o \
+$(SRC)/psphotChoosePSF.$(ARCH).o   \
+$(SRC)/psphotApplyPSF.$(ARCH).o	   \
+$(SRC)/psphotFitGalaxies.$(ARCH).o \
+$(SRC)/psphotOutput.$(ARCH).o      \
+$(SRC)/psphotMarkPSF.$(ARCH).o     \
+$(SRC)/psphotSubtractPSF.$(ARCH).o \
+$(SRC)/psphotSortBySN.$(ARCH).o    \
+$(SRC)/psphotDefinePixels.$(ARCH).o \
+$(SRC)/psphotMagnitudes.$(ARCH).o  \
+$(SRC)/psLibUtils.$(ARCH).o	   \
+$(SRC)/psLine.$(ARCH).o		   \
+$(SRC)/psMinimize.$(ARCH).o	   \
+$(SRC)/psPolynomials.$(ARCH).o	   \
+$(SRC)/psEllipse.$(ARCH).o         \
+$(SRC)/psModulesUtils.$(ARCH).o	   \
+$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
+$(SRC)/pmPSF.$(ARCH).o             \
+$(SRC)/pmPSFtry.$(ARCH).o          \
+$(SRC)/psImageData.$(ARCH).o	   \
+$(SRC)/pmModelGroup.$(ARCH).o       \
+$(SRC)/pmObjects_EAM.$(ARCH).o
+
+MODELS = \
+$(SRC)/models/pmModel_GAUSS.c      \
+$(SRC)/models/pmModel_PGAUSS.c     \
+$(SRC)/models/pmModel_QGAUSS.c     \
+$(SRC)/models/pmModel_SGAUSS.c     \
+$(SRC)/models/pmModel_TGAUSS.c      
+
+$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS)
+
+# deprecated
+
+FITSOURCE = \
+$(SRC)/fitsource.$(ARCH).o \
+$(SRC)/onesource.$(ARCH).o \
+$(SRC)/fs_args.$(ARCH).o \
+$(SRC)/psphot-utils.$(ARCH).o \
+$(SRC)/psPolynomials.$(ARCH).o \
+$(SRC)/psUtils.$(ARCH).o \
+$(SRC)/setup.$(ARCH).o \
+$(SRC)/LocalSky.$(ARCH).o
+
+psphot: $(BIN)/psphot.$(ARCH)
+$(BIN)/psphot.$(ARCH) : $(PSPHOT)
+$(PSPHOT): $(SRC)/psphot.h
+
+fitsource: $(BIN)/fitsource.$(ARCH)
+$(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
+$(FITSOURCE): $(SRC)/psphot.h
+
+INSTALL = psphot
+
+# 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: /trunk/psastro/src/pmAstrom.c
===================================================================
--- /trunk/psastro/src/pmAstrom.c	(revision 5360)
+++ /trunk/psastro/src/pmAstrom.c	(revision 5360)
@@ -0,0 +1,277 @@
+# include "pmAstrom.h"
+
+// input: st1, st2 are psArray of psPlane? pmAstromObj
+
+static double maxOffset;  // maximum allowed offset between lists, in raw pixels
+static double Scale;      // grid pixel scale
+static double Offset;     // deltas to pixels
+
+bool p_pmAstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
+
+    if (fabs(dP) > maxOffset) return false;
+    if (fabs(dQ) > maxOffset) return false;
+
+    *pOut = dP / Scale + Offset;
+    *qOut = dQ / Scale + Offset;
+    return true;
+}
+
+/* Illustration of the grid bins
+   dP        bin
+   -35:-25 -> 0     bin = dP / Scale + Offset
+   -25:-15 -> 1     Scale = 10
+   -15:-05 -> 2     Offset = 3.5
+   -05:+05 -> 3     nPix = 3 (maxOffset = 35 = (nPix + 0.5)*dPix
+   +05:+15 -> 4     dPix = 10
+   +15:+25 -> 5
+   +25:+35 -> 6
+
+   maxOffsetRequest = 30
+   nPix = (int) (maxOffset / dPix + 0.5);
+   maxOffset = (nPix + 0.5)*Scale;
+*/
+
+// apply the measured FPA offset and rotation (stat) to the fpa astrom structures
+psFPA *pmAstromApplyGridMatch (psFPA *fpa, pmAstromGridMatchStat stat) {
+
+    // stat.angle modifies fpa.toTangentPlane (effective angle)
+    // stat.dP, stat.dQ modifies fpa.projection (Ro, Do)
+
+    return (fpa);
+}
+
+pmAstromGridMatchStat pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
+
+    pmAstromGridMatchStat minStat, newStat;
+
+    // find center of the st2 field (focal-plane coords)
+    pMin = 1e10; pMax = -1e10;
+    qMin = 1e10; qMax = -1e10;
+    for (int i = 0; i < st2->n; i++) {
+	ob2 = (pmAstromObj *)st2->data[i];
+	pMin = PS_MIN (ob2->P, pMin);
+	pMax = PS_MAX (ob2->P, pMax);
+	qMin = PS_MIN (ob2->Q, qMin);
+	qMax = PS_MAX (ob2->Q, qMax);
+    }
+    pCenter = 0.5*(pMin + pMax);
+    qCenter = 0.5*(qMin + qMax);
+
+    double minAngle = psMetadataLookupF32 (&status, config, "MIN_ANGLE");
+    double maxAngle = psMetadataLookupF32 (&status, config, "MAX_ANGLE");
+    double delAngle = psMetadataLookupF32 (&status, config, "DEL_ANGLE");
+
+    minStat.minMetric = 1e10;
+    for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
+	st2r = pmAstromRotateObj (st2, angle, pCenter, qCenter);
+	newStat = pmAstromGridMatchAngle (st1, st2r, config);
+	newStat.angle = angle;
+	if (newStat.minMetric < minStat.minMetric) {
+	    minStat = newStat;
+	}
+    }
+    return (minStat);
+}
+
+// rotate the focal-plane coordinates about the center coordinate
+// angle specified in radians
+psArray *pmAstromRotateObj (psArray *old, double angle, double pCenter, double qCenter) {
+
+    pmAstromObj *ob;
+
+    psArray *new = psArrayAlloc (old->n);
+
+    cs = cos(angle);
+    sn = sin(angle);
+    
+    for (int i = 0; i < old->n; i++) {
+
+	oldObj = (pmAstromObj *)old->data[i];
+	newObj = pmAstromObjCopy (oldObj);
+
+	P = oldObj->P - pCenter;
+	Q = oldObj->Q - qCenter;
+	
+	newObj->P = P*cs + Q*sn;
+	newObj->Q = Q*cs - P*sn;
+	
+	new->data[i] = newObj;
+    }
+    return (new);
+}
+
+// match the two lists using the binned delta-delta max
+pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
+
+    // input lists must be projected to common focal plane
+    pmAstromGridMatchStats matchStats;
+    
+    double gridOffset = psMetadataLookupF32 (&status, config, "GRID_OFFSET");
+    double gridScale  = psMetadataLookupF32 (&status, config, "GRID_SCALE");
+
+    // set the static scaling factors
+    nPixHalf = (int)(gridOffset / gridScale + 0.5);
+    nPix     = 2*nPixHalf + 1;
+
+    maxOffpix = gridScale * (nPix + 0.5);
+    Offset    = maxOffpix / gridScale; 
+    Scale     = gridScale;
+
+    // ** can we assume the allocated image is init-ed?
+    gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
+    gridDP = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    gridDQ = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32);
+    psImageInit (gridNP);
+    psImageInit (gridDP);
+    psImageInit (gridDQ);
+    psImageInit (gridD2);
+
+    // short-cut names for grid images
+    NP = gridNP->data.F32;
+    DP = gridDP->data.F32;
+    DQ = gridDQ->data.F32;
+    D2 = gridD2->data.F32;
+
+    // accumulate grids for focal plane (L,M) matches
+    for (int i = 0; i < st1->n; i++) {
+	ob1 = (pmAstromObj *)st1->data[i];
+	for (int j = 0; j < st2->n; j++) {
+	    ob2 = (pmAstromObj *)st2->data[i];
+	    dP = ob1->P - ob2->P;
+	    if (fabs(dP) > maxOffset) continue;
+
+	    dQ = ob1->Q - ob2->Q;
+	    if (fabs(dQ) > maxOffset) continue;
+
+	    // find bin coordinates for this delta-delta
+	    p_pmAstromGridBin (&iP, &iQ, dP, dQ);
+		
+	    // accumulate dP2, dQ2? 
+	    NP[iQ][iP] ++;
+	    DP[iQ][iP] += dP;
+	    DQ[iQ][iP] += dQ;
+	    D2[iQ][iP] += PS_SQR(dP) + PS_SQR(dQ);
+	}
+    }
+
+    // find the max pixel
+    stats = psStatsAlloc (PS_STAT_MAX);
+    stats = psImageStats (stats, gridNP);
+
+    // only check bins with at least 1/2 of max bin
+    minNpt = 0.5*stats->max;
+
+    // find the 'best' bin
+    minMetric = minVar = 1e10;
+    minP = minQ = -1;
+    for (int j = 0; j < gridNP->nY; j++) {
+	for (int i = 0; i < gridNP->nX; i++) {
+
+	    if (NP[j][i] < minNpts) continue;
+
+	    // this metric emphasize a narrow peak with lots of sources over one with few.
+	    var = fabs((D2[j][i]/NP[j][i]) - PS_SQR(DP[j][i]/NP[j][i]) - PS_SQR(DQ[j][i]/NP[j][i]));
+	    metric = var / PS_SQR(PS_SQR(NP[j][i]));
+	    
+	    if (metric < minMetric) {
+		minMetric = metric;
+		minVar    = var;
+		minP      = i;
+		minQ      = j;
+	    }
+	}
+    }
+
+    // convert the bin to delta-delta
+    matchStats.dP        = DP[minQ][minP] / NP[minQ][minP];
+    matchStats.dQ        = DQ[minQ][minP] / NP[minQ][minP];
+    matchStats.minMetric = minMetric;
+    matchStats.minVar    = minVar;
+    matchStats.nMatch    = NP[minQ][minP];
+
+    return (matchStats);
+}
+
+psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) {
+
+    // assume the lists are sorted, or sort in place?
+
+    // sort st1 by P
+    // sort st2 by P
+
+    int i = 0;
+    int j = 0;
+
+    double RADIUS = psMetadataLookupR32 (&status, myHeader, "ASTROM_MATCH_RADIUS");
+    double RADIUS_SQR = PS_SQR (RADIUS);
+
+    psArray *matches = psArrayAlloc (100);
+    matches->n = 0;
+
+    while ((i < st1->n) && (j < st2->n)) {
+
+	dP = ((pmAstromObj *)st1->data[i])->P - ((pmAstromObj *)st2->data[i])->P;
+	if (dP < -RADIUS) {
+	    i++;
+	    continue;
+	}
+	if (dP > +RADIUS) {
+	    j++;
+	    continue;
+	}
+
+	jStart = j;
+	while ((dX > -RADIUS) && (j < st2->n)) {
+	    
+	    dP = ((pmAstromObj *)st1->data[i])->P - ((pmAstromObj *)st2->data[i])->P;
+	    dQ = ((pmAstromObj *)st1->data[i])->Q - ((pmAstromObj *)st2->data[i])->Q;
+	    dR = dP*dP + dQ*dQ;
+
+	    if (dR > RADIUS_SQR) {
+		j++;
+		continue;
+	    }
+	    
+	    // got a match; add to output list
+	    pmAstromMatch *match = pmAstromMatchAlloc (i, j);
+	    psArrayAdd (matches, match, 100);
+
+	    j++;
+	}
+	j = jStart;
+	i++;
+    }
+    return (match);
+}
+
+pmAstromMatchedListFit (pmFPA *fpa, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+
+    psVector *X = psVectorAlloc (match->n);
+    psVector *Y = psVectorAlloc (match->n);
+    psVector *x = psVectorAlloc (match->n);
+    psVector *y = psVectorAlloc (match->n);
+  
+    // take the matched stars, first fit 
+    for (i = 0; i < match->n; i++) {
+
+	pair = match->data[i];
+	ob1 = st1->data[pair->i1];
+	ob2 = st2->data[pair->i2];
+
+	X->data.F64[i] = ob1->P;
+	Y->data.F64[i] = ob1->Q;
+
+	x->data.F64[i] = ob2->P;
+	y->data.F64[i] = ob2->Q;
+
+	// use one or the other...
+	psPolynomial2D *xt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
+	psPolynomial2D *xt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
+
+	psPolynomial2D *yt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
+	psPolynomial2D *yt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
+    }
+
+    // apply fitted polynomials to fpa
+}
Index: /trunk/psastro/src/pmAstrom.h
===================================================================
--- /trunk/psastro/src/pmAstrom.h	(revision 5360)
+++ /trunk/psastro/src/pmAstrom.h	(revision 5360)
@@ -0,0 +1,35 @@
+# include <stdio.h>
+# include <strings.h>  // for strcasecmp
+# include <unistd.h>   // for unlink
+# include <pslib.h>
+
+psArray *pmAstromGridMatch (psArray *s1, psArray *s2, pmAstromGridMatchOpt *opt);
+
+typedef struct {
+    double X, Y;
+    double P, Q;
+    double L, M;
+    double R, D;
+    double Mag, dMag;
+} pmAstromObj;
+
+typedef struct {
+    int i1;
+    int i2;
+} pmAstromMatch;
+
+typedef struct {
+    double minAngle;
+    double maxAngle;
+    double delAngle;
+} pmAstromOpt;
+
+typedef struct {
+    double angle;
+    double dP;
+    double dQ;
+    double minMetric;
+    double minVar;
+    int    nMatch;
+} pmAstromGridMatchStats;
+
Index: /trunk/psastro/src/psastro.c
===================================================================
--- /trunk/psastro/src/psastro.c	(revision 5360)
+++ /trunk/psastro/src/psastro.c	(revision 5360)
@@ -0,0 +1,51 @@
+# include "psastro.h"
+
+int main (int argc, char **argv) {
+
+    psMetadata *header = NULL;
+    pmAstromGridMatchStat stat;
+    pmAstromMatch *match;
+
+    // load configuration information
+    psMetadata *config = psastroArguments (&argc, argv);
+
+    // define optional parameters (keep in psMetadata?)
+    pmAstromOpt *options = psastroDefineOptions (config);
+
+    // load the input data (cmp file)
+    psArray *rawstars = pmSourcesReadCMP (&header, argv[1]);
+
+    psFPA *fpa = pmBuildFPA (header, config);
+
+    // limit fit to bright stars only 
+    psArray *subset = pspsastroSelectBrightStars (config, rawstars);
+
+    // use the header & config info to project rawstars on the focal plane
+    psastroProjectRawstars (subset, header, config);
+
+    // load the corresponding reference data (DVO command)
+    psArray *refstars = psastroLoadReference (argv[3]);
+
+    // use the header & config info to project refstars on the focal plane
+    psastroProjectRefstars (refstars, header, config);
+
+    // find initial offset / rotation
+    stat = pmAstromGridMatch (subset, refstars, options);
+
+    pmAstromModifyFPA (fpa, stat);
+
+    // use fit result to re-project rawstars
+    psastroProjectRawstars (fpa, subset);
+    psastroProjectRefstars (fpa, refstars);
+    
+    // use small radius to match stars
+    match = pmAstromRadiusMatch (rawstars, refstars, options);
+
+    // fit astrometric terms
+    output = pmAstromMatchedListFit (fpa, subset, refstars, match, options);
+
+    // write out data (cmp file)
+    psastroWriteCMP (fpa, header, rawstars, argv[2]);
+
+    exit (0);
+}
Index: /trunk/psastro/src/psastro.h
===================================================================
--- /trunk/psastro/src/psastro.h	(revision 5360)
+++ /trunk/psastro/src/psastro.h	(revision 5360)
@@ -0,0 +1,55 @@
+# include <stdio.h>
+# include <strings.h>  // for strcasecmp
+# include <unistd.h>   // for unlink
+# include <pslib.h>
+# include "psLibUtils.h"
+# include "pmObjects_EAM.h"
+# include "psModulesUtils.h"
+# include "pmPSF.h"
+# include "pmPSFtry.h"
+# include "pmModelGroup.h"
+
+typedef struct {
+    psImage *image;
+    psImage *mask;
+    psImage *weight;
+    psMetadata *header;
+} eamReadout;
+
+// top-level psphot functions
+psMetadata  *psphotArguments (int *argc, char **argv);
+eamReadout  *psphotSetup (psMetadata *config);
+psStats     *psphotImageStats (eamReadout *imdata, psMetadata *config);
+psArray     *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks);
+pmPSF       *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
+bool         psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool         psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
+void         psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+
+bool         psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE);
+bool         psphotSubtractPSF (pmSource *source);
+int 	     psphotSortBySN (const void **a, const void **b);
+int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+bool 	     psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
+
+psArray     *pmPeaksSigmaLimit (eamReadout *imdata, psMetadata *config, psStats *sky);
+pmModel     *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
+pmModel     *pmSourceSelectModel (pmSource *source);
+
+// eamReadout functions
+eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header);
+
+// output functions
+bool 	     pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	     pmSourcesWriteOBJ  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	     pmSourcesWriteCMP  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	     pmSourcesWriteCMF  (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	     pmSourcesWriteSX   (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+
+bool 	     pmPeaksWriteText (psArray *sources, char *filename);
+bool 	     pmMomentsWriteText (psArray *sources, char *filename);
+bool 	     pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf);
+bool 	     pmModelWriteFLTs (psArray *sources, char *filename);
+bool 	     pmModelWriteNULLs (psArray *sources, char *filename);
+
+int  	     pmSourcesDophotType (pmSource *source);
Index: /trunk/psastro/src/psastroArguments.c
===================================================================
--- /trunk/psastro/src/psastroArguments.c	(revision 5360)
+++ /trunk/psastro/src/psastroArguments.c	(revision 5360)
@@ -0,0 +1,86 @@
+# include "psphot.h"
+
+static void usage (void) ;
+
+psMetadata *psphotArguments (int *argc, char **argv) {
+
+  int N, Nfail;
+  int mode = PS_DATA_STRING | PS_META_REPLACE;
+
+  // basic pslib options
+  fprintf (stderr, "starting... %s\n", psLibVersion());
+  psLogSetFormat ("M");
+  psLogArguments (argc, argv);
+  psTraceArguments (argc, argv);
+
+  // optional mask image - add to config
+  char *mask = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-mask"))) {
+    psArgumentRemove (N, argc, argv);
+    mask = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  // optional weight image - add to config
+  char *weight = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-weight"))) {
+    psArgumentRemove (N, argc, argv);
+    weight = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  // optional output residual image - add to config
+  char *resid = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-resid"))) {
+    psArgumentRemove (N, argc, argv);
+    resid = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  // optional output residual image - add to config
+  char *photcode = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
+    psArgumentRemove (N, argc, argv);
+    photcode = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  if (*argc != 4) usage ();
+
+  // load config information
+  psMetadata *config = psMetadataAlloc ();
+  psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
+  config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
+
+  // identify input image & optional weight & mask images
+  // command-line entries override config-file entries
+  psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
+  psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
+
+  if (mask != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
+  }
+  if (weight != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
+  }
+  if (resid != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
+  }
+  if (photcode != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
+  } else {
+    psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
+  }    
+  return (config);
+}
+
+static void usage (void) {
+
+    fprintf (stderr, "USAGE: psphot (image) (output) (config)\n");
+    fprintf (stderr, "options: \n");
+    fprintf (stderr, "  -mask  (filename)\n");
+    fprintf (stderr, "  -weight (filename)\n");
+    fprintf (stderr, "  -resid (filename)\n");
+    fprintf (stderr, "  -photcode (photcode)\n");
+    exit (2);
+}
Index: /trunk/psastro/src/psastroUtils.c
===================================================================
--- /trunk/psastro/src/psastroUtils.c	(revision 5360)
+++ /trunk/psastro/src/psastroUtils.c	(revision 5360)
@@ -0,0 +1,130 @@
+# include "psastro.h"
+
+// crude function to load CMP data (header + ascii)
+// XXX EAM : assumes fixed line, expects NSTARS entries
+psArrray *pmSourcesReadCMP (psMetadata **header, char *filename) {
+
+    psMetadata *myHeader;
+    char line[80];
+
+    psFits *fits = psFitsAlloc (filename);
+    myHeader = psFitsReadHeader (fits);
+    psFree (fits);
+
+    // how many lines in the header?
+    // XXX EAM : is this calculation robust?
+    nLines = header->list->n;
+    nBytes = nLines * 80;
+    if (nBytes % 2880) {
+	nBlock = 1 + (int)(nBytes / 2880);
+	nBytes = nBlock * 2880;
+    }
+
+    // re-open, seek to end of header
+    FILE *f = fopen (filename, "r");
+    if (f == NULL) {
+	psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
+	return NULL;
+    }
+    fseek (f, nBytes, SEEK_SET);
+
+    // prepare array to store data
+    nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
+    psArray *stars = psArrayAlloc (nStars);
+    stars->n = 0;
+
+    // we have fixed bytes / line : use that info
+    for (i = 0; i < nStars; i++) {
+	fread (line, 1, 67, f);
+	sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
+	
+	pmAstromObj *obj = pmAstromObjAlloc (X, Y, Mag, dMag);
+	psArrayAdd (stars, obj, 100);
+    }
+    fclose (f);
+
+    *header = myHeader;
+    return (stars);
+}
+
+// sort by mag (descending)
+int psastroSortByMag (const void **a, const void **b)
+{
+    pmSource *A = *(pmSource **)a;
+    pmSource *B = *(pmSource **)b;
+
+    psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
+    psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
+    if (isnan (fA)) fA = 0;
+    if (isnan (fB)) fB = 0;
+
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (-1);
+    if (diff < FLT_EPSILON) return (+1);
+    return (0);
+}
+
+psArray *psastroSelectBrightStars (psMetadata *config, psArray *stars) {
+
+    stars = psArraySort (stars, psastroSortByMag);
+
+    // add exclusions for objects on some basis?
+    int MAX_NSTARS = pmMetadataLookupS32 (&status, config, "MAX_NSTARS");
+
+    nSubset = PS_MIN (MAX_NSTARS, stars->n);
+
+    psArray *subset = psArrayAlloc (nSubset);
+
+    for (int i = 0; i < nSubset; i++) {
+	subset->data[i] = stars->data[i];
+    }
+
+    return (subset);
+}
+
+// fake this by loading from a text file
+psArray *psastroLoadReference (char *filename) {
+
+    psMetadata *myHeader;
+    char line[80];
+
+    psFits *fits = psFitsAlloc (filename);
+    myHeader = psFitsReadHeader (fits);
+    psFree (fits);
+
+    // how many lines in the header?
+    // XXX EAM : is this calculation robust?
+    nLines = header->list->n;
+    nBytes = nLines * 80;
+    if (nBytes % 2880) {
+	nBlock = 1 + (int)(nBytes / 2880);
+	nBytes = nBlock * 2880;
+    }
+
+    // re-open, seek to end of header
+    FILE *f = fopen (filename, "r");
+    if (f == NULL) {
+	psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
+	return NULL;
+    }
+    fseek (f, nBytes, SEEK_SET);
+
+    // prepare array to store data
+    nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
+    psArray *stars = psArrayAlloc (nStars);
+    stars->n = 0;
+
+    // we have fixed bytes / line : use that info
+    for (i = 0; i < nStars; i++) {
+	fread (line, 1, 67, f);
+	sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
+	
+	pmAstromObj *obj = pmAstromObjAlloc (X, Y, Mag, dMag);
+	psArrayAdd (stars, obj, 100);
+    }
+    fclose (f);
+
+    *header = myHeader;
+    return (stars);
+}
+