Changeset 5993

Timestamp:

Jan 15, 2006, 8:30:16 AM (20 years ago)

Author:

eugene

Message:

API cleanup, removed old test files, some re-organization, changed FLT/GAL to EXT

Location:

trunk/psphot

Files:

• Makefile (modified) (1 diff)
• src/pmSourceFitFixed.c (deleted)
• src/pmSourceFitSet.c (modified) (1 diff)
• src/psBicube.c (deleted)
• src/psLibUtils.c (deleted)
• src/psModulesUtils.c (modified) (1 diff)
• src/psPolynomialUtils.c (added)
• src/psphot.c (modified) (1 diff)
• src/psphot.h (modified) (3 diffs)
• src/psphotApResid.c (modified) (3 diffs)
• src/psphotApplyPSF.c (modified) (1 diff)
• src/psphotArguments.c (modified) (2 diffs)
• src/psphotBasicDeblend.c (modified) (3 diffs)
• src/psphotBlendFit.c (modified) (2 diffs)
• src/psphotChoosePSF.c (modified) (4 diffs)
• src/psphotDefinePixels.c (modified) (2 diffs)
• src/psphotEnsemblePSF.c (modified) (7 diffs)
• src/psphotEvalFLT.c (modified) (3 diffs)
• src/psphotEvalPSF.c (modified) (2 diffs)
• src/psphotFitGalaxies.c (modified) (5 diffs)
• src/psphotFixedPSF.c (deleted)
• src/psphotFullFit.c (modified) (2 diffs)
• src/psphotMagnitudes.c (modified) (1 diff)
• src/psphotMarkPSF.c (deleted)
• src/psphotModelTest.c (modified) (3 diffs)
• src/psphotOutput.c (modified) (16 diffs)
• src/psphotRadiusChecks.c (modified) (1 diff)
• src/psphotReapplyPSF.c (deleted)
• src/psphotReplaceUnfit.c (modified) (1 diff)
• src/psphotRoughClass.c (modified) (1 diff)
• src/psphotSourceFits.c (modified) (7 diffs)
• src/psphotSparseMatrix.c (deleted)
• src/psphotSubtractPSF.c (deleted)
• src/psphotTest.c (modified) (2 diffs)
• src/psphotTestArguments.c (modified) (2 diffs)

trunk/psphot/Makefile

@@ -24 +24 @@
 $(SRC)/psphotArguments.$(ARCH).o    \
 $(SRC)/psphotSetup.$(ARCH).o        \
+$(SRC)/psphotModelTest.$(ARCH).o    \
 $(SRC)/psphotImageStats.$(ARCH).o   \
+$(SRC)/psphotImageBackground.$(ARCH).o \
+$(SRC)/psphotFindPeaks.$(ARCH).o    \
 $(SRC)/psphotSourceStats.$(ARCH).o  \
+$(SRC)/psphotRoughClass.$(ARCH).o   \
+$(SRC)/psphotBasicDeblend.$(ARCH).o \
 $(SRC)/psphotChoosePSF.$(ARCH).o    \
+$(SRC)/psphotEnsemblePSF.$(ARCH).o  \
+$(SRC)/psphotFullFit.$(ARCH).o      \
+$(SRC)/psphotBlendFit.$(ARCH).o     \
 $(SRC)/psphotApplyPSF.$(ARCH).o     \
-$(SRC)/psphotFixedPSF.$(ARCH).o     \
-$(SRC)/psphotEnsemblePSF.$(ARCH).o  \
-$(SRC)/psphotReapplyPSF.$(ARCH).o   \
 $(SRC)/psphotFitGalaxies.$(ARCH).o  \
-$(SRC)/psphotOutput.$(ARCH).o       \
-$(SRC)/psphotMarkPSF.$(ARCH).o      \
+$(SRC)/psphotReplaceUnfit.$(ARCH).o \
+$(SRC)/psphotApResid.$(ARCH).o      \
+$(SRC)/psphotOutput.$(ARCH).o       
+
+SUPPORT = \
+$(SRC)/psphotModelGroupInit.$(ARCH).o  \
+$(SRC)/psphotGrowthCurve.$(ARCH).o  \
+$(SRC)/psphotEvalPSF.$(ARCH).o      \
+$(SRC)/psphotEvalFLT.$(ARCH).o      \
+$(SRC)/psphotSourceFits.$(ARCH).o   \
 $(SRC)/psphotSortBySN.$(ARCH).o     \
 $(SRC)/psphotDefinePixels.$(ARCH).o \
 $(SRC)/psphotMagnitudes.$(ARCH).o   \
-$(SRC)/psphotRadiusChecks.$(ARCH).o \
-$(SRC)/psphotReplaceUnfit.$(ARCH).o \
-$(SRC)/psphotEvalPSF.$(ARCH).o      \
-$(SRC)/psphotEvalFLT.$(ARCH).o      \
-$(SRC)/psphotFullFit.$(ARCH).o      \
+$(SRC)/psphotRadiusChecks.$(ARCH).o 
+
+PSMODULES = \
+$(SRC)/psLine.$(ARCH).o             \
+$(SRC)/psSparse.$(ARCH).o           \
+$(SRC)/psPolynomialUtils.$(ARCH).o  \
+$(SRC)/psImageData.$(ARCH).o        \
+$(SRC)/psModulesUtils.$(ARCH).o     \
 $(SRC)/pmSourceContour.$(ARCH).o    \
-$(SRC)/psLine.$(ARCH).o             \
-$(SRC)/psModulesUtils.$(ARCH).o     \
-$(SRC)/pmPeaksSigmaLimit.$(ARCH).o  \
-$(SRC)/pmSourceFitFixed.$(ARCH).o   \
-$(SRC)/psSparse.$(ARCH).o           \
-$(SRC)/psImageData.$(ARCH).o        \
-$(SRC)/psphotModelTest.$(ARCH).o    \
-$(SRC)/psphotImageBackground.$(ARCH).o \
-$(SRC)/psphotBasicDeblend.$(ARCH).o    \
-$(SRC)/psphotModelGroupInit.$(ARCH).o  \
-$(SRC)/pmModelFitSet.$(ARCH).o  \
-$(SRC)/pmSourceFitSet.$(ARCH).o        \
-$(SRC)/psphotSourceFits.$(ARCH).o        \
-$(SRC)/psphotBlendFit.$(ARCH).o        \
-$(SRC)/psphotApResid.$(ARCH).o \
-$(SRC)/psphotAssessPSF.$(ARCH).o \
-$(SRC)/psBicube.$(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)
-
-MODELTEST = \
-$(SRC)/modelTestFitSource.$(ARCH).o \
-$(SRC)/modelTestArguments.$(ARCH).o \
-$(SRC)/psphotModelGroupInit.$(ARCH).o  \
-$(SRC)/psModulesUtils.$(ARCH).o     \
-$(SRC)/psImageData.$(ARCH).o        \
-$(SRC)/psphotSetup.$(ARCH).o        \
-$(SRC)/psphotDefinePixels.$(ARCH).o \
-$(SRC)/psphotModelTest.$(ARCH).o
+$(SRC)/pmModelFitSet.$(ARCH).o      \
+$(SRC)/pmSourceFitSet.$(ARCH).o     
 
 TEST = \
-$(SRC)/psphotTest.$(ARCH).o \
-$(SRC)/psphotTestArguments.$(ARCH).o \
-$(SRC)/psphotApResid.$(ARCH).o
+$(SRC)/psphotTest.$(ARCH).o     \
+$(SRC)/psBicube.$(ARCH).o           \
+$(SRC)/psphotTestArguments.$(ARCH).o
 
 psphot: $(BIN)/psphot.$(ARCH)
-$(BIN)/psphot.$(ARCH) : $(PSPHOT)
+$(BIN)/psphot.$(ARCH) : $(PSPHOT) $(PSMODULES) $(SUPPORT)
 $(PSPHOT) : $(SRC)/psphot.h
 
-modeltest.install: psphotModelTest.install
-modeltest: $(BIN)/psphotModelTest.$(ARCH)
-$(BIN)/psphotModelTest.$(ARCH) : $(MODELTEST)
-$(MODELTEST): $(SRC)/psphot.h
+test.install: psphotTest.install
+test: $(BIN)/psphotTest.$(ARCH)
+$(DESTBIN)/test : $(DESTBIN)/psphotTest
 
-test-install: psphotTest.install
-test: $(BIN)/psphotTest.$(ARCH)
 $(BIN)/psphotTest.$(ARCH) : $(TEST)
 $(TEST): $(SRC)/psphot.h

trunk/psphot/src/pmSourceFitSet.c

@@ -91 +91 @@
     }
 
-    // PSF model only fits first 4 parameters, FLT model fits all
+    // PSF model only fits first 4 parameters, EXT model fits all
     int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1;
 

trunk/psphot/src/psModulesUtils.c

@@ -123 +123 @@
     return (new);
 }
+
+float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) {
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    float flux;
+
+    psImage *Pi = Mi->pixels;
+    psImage *Pj = Mj->pixels;
+
+    psImage *Ti = Mi->mask;
+    psImage *Tj = Mj->mask;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+    
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+    
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.U8[yi][xi]) continue;
+            if (Tj->data.U8[yj][xj]) continue;
+            flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
+        }
+    }
+    return (flux);
+}

trunk/psphot/src/psphot.c

@@ -67 +67 @@
 
       case 3:
-        // fit extended objects with galaxy models
         psphotApplyPSF (imdata, config, sources, psf, sky);
         break;
 
       case 4:
-        // fit extended objects with galaxy models
         psphotApplyPSF (imdata, config, sources, psf, sky);
-        psphotFitGalaxies (imdata, config, sources, sky);
+        psphotFitExtended (imdata, config, sources, sky);
         break;
     }

trunk/psphot/src/psphot.h

@@ -24 +24 @@
 psMetadata     *psphotArguments (int *argc, char **argv);
 eamReadout     *psphotSetup (psMetadata *config);
+bool            psphotModelTest (eamReadout *imdata, psMetadata *config);
 psStats        *psphotImageStats (eamReadout *imdata, psMetadata *config);
+psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky);
 psArray        *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky);
 psArray        *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks);
 bool            psphotRoughClass (psArray *sources, psMetadata *config);
+bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
 pmPSF          *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
-bool            psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool            psphotFixedPSF (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);
+// optional object analysis steps
+bool            psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool            psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool            psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool            psphotReplaceUnfit (psArray *sources);
+bool            psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool            psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
+bool            psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf);
+
+// basic support functions
+pmModel        *pmModelCopy (pmModel *model);
+pmModel        *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
+float           pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
+psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
+void            psphotModelGroupInit (void);
 int             psphotSortBySN (const void **a, const void **b);
 int             psphotSortByY (const void **a, const void **b);
-int             psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+bool            psphotGrowthCurve (eamReadout *imdata, pmPSF *psf);
+void            psphotTestArguments (int *argc, char **argv);
+
+// functions to set the correct source pixels
+bool            psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
+bool            psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model);
+bool            psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type);
+bool            psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model);
 bool            psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
 bool            psphotRedefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius);
-void            psphotModelGroupInit (void);
-
-bool            pmSourceFitFixed (pmSource *source, pmModel *model);
-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
@@ -57 +69 @@
 bool            pmSourcesWriteSX   (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
 
+bool            pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf);
+bool            pmModelWriteEXTs (psArray *sources, char *filename);
+bool            pmModelWriteNULLs (psArray *sources, char *filename);
 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);
 
+bool            psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image);
+bool            psphotDumpMoments (psMetadata *config, psArray *sources);
+int             psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+bool            psphotUpdateHeader (psMetadata *header, psMetadata *config);
 int             pmSourcesDophotType (pmSource *source);
+bool            psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key);
 
-// psphotModelTest functions
-bool            psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-float           psphotCrossProduct (pmSource *Mi, pmSource *Mj);
-psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky);
-bool            psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+// PSF / DBL / EXT evaluation functions
+bool            psphotEvalPSF (pmSource *source, pmModel *model);
+bool            psphotEvalDBL (pmSource *source, pmModel *model);
+bool            psphotEvalEXT (pmSource *source, pmModel *model);
 
-pmModel        *pmModelCopy (pmModel *model);
-psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
-bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
+// functions to support simultaneous multi-source fitting
+bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
+bool            pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
+psF32           pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
+bool            pmModelFitSetInit (pmModelType type);
 
-bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool psphotInitLimitsPSF (psMetadata *config);
-bool psphotEvalPSF (pmSource *source, pmModel *model);
-bool psphotEvalDBL (pmSource *source, pmModel *model);
-bool psphotEvalFLT (pmSource *source, pmModel *model);
-bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
-bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model);
-bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type);
-bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model);
-bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image);
-bool psphotReplaceUnfit (psArray *sources);
-bool psphotDumpMoments (psMetadata *config, psArray *sources);
-bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf);
-bool psphotAssessPSF (eamReadout *imdata, psMetadata *config, pmPSF *psf);
+//  functions to support the source fitting process
+bool            psphotInitLimitsPSF (psMetadata *config);
+bool            psphotInitLimitsEXT (psMetadata *config, psStats *sky);
+bool            psphotFitBlend (eamReadout *imdata, pmSource *source);
+bool            psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources);
+bool            psphotFitPSF (eamReadout *imdata, pmSource *source);
+pmModel        *psphotFitEXT (eamReadout *imdata, pmSource *source);
+psArray        *psphotFitDBL (eamReadout *imdata, pmSource *source);
 
-bool psphotWritePSF (pmPSF *psf, char *filename);
-pmPSF *psphotReadPSF (char *filename);
+// eamReadout functions
+eamReadout     *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header);
 
-bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...);
-bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...);
-bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...);
-psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder);
-psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder);
-psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder);
-bool psphotModelTest (eamReadout *imdata, psMetadata *config);
+// bicubic interpolation
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
+psPlane         psImageBicubeMin (psPolynomial2D *poly);
 
-bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
-bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
-psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
-bool pmModelFitSetInit (pmModelType type);
-
-bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-
-bool psphotInitLimitsFLT (psMetadata *config, psStats *sky);
-
-bool psphotFitPSF (eamReadout *imdata, pmSource *source);
-bool psphotFitBlend (eamReadout *imdata, pmSource *source);
-bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources);
-
-pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source);
-psArray *psphotFitDBL (eamReadout *imdata, pmSource *source);
-
-psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
-psPlane psImageBicubeMin (psPolynomial2D *poly);
-void psphotTestArguments (int *argc, char **argv);
-
-bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf);
-
+// optional mode for clip fit?
 psPolynomial4D *psVectorChiClipFitPolynomial4D(
     psPolynomial4D *poly,
@@ -134 +121 @@
     const psVector *z,
     const psVector *t);
-
-// XXX deprecated
-// bool psphotApResidFullFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag);
-// bool psphotApResidPolyFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag);

trunk/psphot/src/psphotApResid.c

@@ -1 +1 @@
 # include "psphot.h"
-
-psPolynomial4D *psVectorChiClipFitPolynomial4D(
-    psPolynomial4D *poly,
-    psStats *stats,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x,
-    const psVector *y,
-    const psVector *z,
-    const psVector *t)
-{
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(fErr, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
-
-    // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
-    if (isfinite(stats->max)) {
-        maxClipSigma = +fabs(stats->max);
-    } else {
-        maxClipSigma = +fabs(stats->clipSigma);
-    }
-    if (isfinite(stats->min)) {
-        minClipSigma = -fabs(stats->min);
-    } else {
-        minClipSigma = -fabs(stats->clipSigma);
-    }
-    psVector *fit   = NULL;
-    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
-
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-
-    for (int N = 0; N < stats->clipIter; N++) {
-        int Nkeep = 0;
-
-        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
-        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
-        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
-
-        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
-        psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
-
-        // set mask if pts are not valid
-        // we are masking out any point which is out of range
-        // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
-            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
-                continue;
-            }
-            float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev);
-            if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) {
-                if (mask != NULL) {
-                    mask->data.U8[i] |= 0x01;
-                }
-                continue;
-            }
-            if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) {
-                if (mask != NULL) {
-                    mask->data.U8[i] |= 0x01;
-                }
-                continue;
-            }
-            Nkeep ++;
-        }
-
-        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
-                 Nkeep, x->n);
-
-        stats->clippedNvalues = Nkeep;
-        psFree (fit);
-    }
-    // Free local temporary variables
-    psFree (resid);
-
-    if (poly == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-    return(poly);
-}
 
 // measure the aperture residual statistics
@@ -347 +239 @@
     psf->nApResid = residStats->clippedNvalues;
 
-    /*
-      (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
-      (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
-      (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
-    */
+    // save results for later output
+    psMetadataAdd (config, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
+    psMetadataAdd (config, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
+    psMetadataAdd (config, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
+    psMetadataAdd (config, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
+    psMetadataAdd (config, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss);
+    psMetadataAdd (config, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
 
     psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
@@ -370 +264 @@
 }
 
+/*
+  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+*/
+

trunk/psphot/src/psphotApplyPSF.c

@@ -25 +25 @@
 
         // use the source moments, etc to guess basic model parameters
-        pmModel *modelFLT = pmSourceModelGuess (source, psf->type); 
+        pmModel *modelEXT = pmSourceModelGuess (source, psf->type); 
 
         // set PSF parameters for this model
-        pmModel *model = pmModelFromPSF (modelFLT, psf);
-        psFree (modelFLT);
+        pmModel *model = pmModelFromPSF (modelEXT, psf);
+        psFree (modelEXT);
 
         source->modelPSF = model;

trunk/psphot/src/psphotArguments.c

@@ -123 +123 @@
     psMetadataAdd (config, PS_LIST_TAIL, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
 
-    config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
+    // config file values override defaults set here
+    config = psMetadataConfigParse (config, &Nfail, argv[3], TRUE);
 
     psMetadataItem *item = NULL;
@@ -129 +130 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
         psMetadataAddItem (config, item, PS_LIST_TAIL, PS_META_REPLACE);
-        // psFree (item);
+        psFree (item);
     }
-    // psFree (iter);
+    psFree (iter);
 
     // identify input image & optional weight & mask images

trunk/psphot/src/psphotBasicDeblend.c

@@ -36 +36 @@
     psArray *overlap = psArrayAlloc (100);
 
-    // XXX make sure this results in decreasing, not increasing, SN
+    // examine sources in decreasing SN order
     for (int i = sources->n - 1; i >= 0; i--) {
         N = index->data.U32[i];
@@ -73 +73 @@
         // generate source contour (1/4 peak counts)
         if (overlap->n > 0) {
-            // XXX EAM : make the 1/4 user-defined.
-            // XXX keep threshold from dropping too low (N*sky.sigma)
+            // set the threshold based on user inputs
             threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; 
             threshold = PS_MAX (threshold, minThreshold);
 
-            // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now
+            // generate a basic contour (set of x,y coordinates at-or-below flux level)
             psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold);
             if (contour == NULL) continue;
@@ -97 +96 @@
                     if (xv->data.F32[j+1] < testSource->peak->x) break;
 
-                    # if (0)
-                    int xp0 = source->moments->x - source->pixels->col0;
-                    int xp1 = source->peak->x - source->pixels->col0;
-                    int xp2 = testSource->moments->x - source->pixels->col0;
-                    int xp3 = testSource->peak->x - source->pixels->col0;
-
-                    int yp0 = source->moments->y - source->pixels->row0;
-                    int yp1 = source->peak->y - source->pixels->row0;
-                    int yp2 = testSource->moments->y - testSource->pixels->row0;
-                    int yp3 = testSource->peak->y - testSource->pixels->row0;
-                    
-                    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n",
-                             source->peak->x, source->peak->y,
-                             source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1],
-                             testSource->peak->x, testSource->peak->y,
-                             testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
-                       );
-                    # endif
-                    
                     testSource->mode |= PM_SOURCE_BLEND;
 

trunk/psphot/src/psphotBlendFit.c

@@ -1 +1 @@
 # include "psphot.h"
 
-bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
-{ 
+// XXX I don't like this name
+bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
 
     psTimerStart ("psphot");
@@ -10 +10 @@
     
     psphotInitLimitsPSF (config);
-    psphotInitLimitsFLT (config, sky);
+    psphotInitLimitsEXT (config, sky);
     psphotInitRadiusPSF (config, sky, psf->type);
 

trunk/psphot/src/psphotChoosePSF.c

@@ -2 +2 @@
 
 // try PSF models and select best option
+pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) { 
 
-pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) 
-{ 
     bool            status;
     char           *modelName;
@@ -55 +54 @@
     }
     psFree (iter);
-    // psFree (list); XXX EAM - is list freed with iter?
+    // psFree (list); XXX is list freed with iter?
     psFree (stars);
 
@@ -78 +77 @@
     for (int i = 0; i < try->sources->n; i++) {
         pmSource *source = try->sources->data[i];
-        if (try->mask->data.U8[i] & PSFTRY_MASK_FLT_FAIL) {
+        if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) {
             source->mode &= ~PM_SOURCE_PSFSTAR;
         }
@@ -86 +85 @@
     pmPSF *psf = psMemCopy(try->psf);
     psFree (models);
+
+    psMetadataAdd (config, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "number of stars used to make PSF", psf->nPSFstars);
 
     modelName = pmModelGetType (psf->type);

trunk/psphot/src/psphotDefinePixels.c

@@ -49 +49 @@
     if (extend) {
         // Grab a new subimage
-        // psFree (mySource->pixels);
-        // psFree (mySource->weight);
-        // psFree (mySource->mask);
+        psFree (mySource->pixels);
+        psFree (mySource->weight);
+        psFree (mySource->mask);
 
-        // fprintf (stderr, "extend %f,%f\n", x, y);
         mySource->pixels = psImageSubset(imdata->image,  newRegion);
         mySource->weight = psImageSubset(imdata->weight, newRegion);
@@ -62 +61 @@
     return extend;
 }
-
-
-//**** this function was ALWAYS restricting to model->radius
-//     should have left it larger (PSF_FIT_RAD) for the faint sources

trunk/psphot/src/psphotEnsemblePSF.c

@@ -1 +1 @@
 # include "psphot.h"
-
-psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
-
-    int ix = x - image->col0;
-    int iy = y - image->row0;
-
-    psF32 *Fm = &image->data.F32[iy - 1][ix];
-    psF32 *Fo = &image->data.F32[iy + 0][ix];
-    psF32 *Fp = &image->data.F32[iy + 1][ix];
-
-    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
-    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
-    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
-    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
-    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
-
-    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
-    poly->mask[2][2] = 1;
-    poly->mask[1][2] = 1;
-    poly->mask[2][1] = 1;
-
-    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
-
-    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
-    poly->coeff[0][1] = (Fyp - Fym)/6.0;
-    
-    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
-    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
-    
-    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
-    
-    return (poly);
-}
-
-psPlane psImageBicubeMin (psPolynomial2D *poly) {
-
-    psPlane min;
-
-    min.xErr = min.yErr = 0;
-
-    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
-
-    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
-    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
-    return (min);
-}
 
 bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
@@ -107 +61 @@
         }
 
-        // XXX EAM : add option to use FLT or PSF form
+        // XXX EAM : add option to use EXT or PSF form
         // use the source moments, etc to guess basic model parameters
-        pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); 
+        pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 
         if (inSource->mode &  PM_SOURCE_SATSTAR) {
-            modelFLT->params->data.F32[2] = inSource->moments->x;
-            modelFLT->params->data.F32[3] = inSource->moments->y;
+            modelEXT->params->data.F32[2] = inSource->moments->x;
+            modelEXT->params->data.F32[3] = inSource->moments->y;
         } else {
             // peak-up on peak (for non-sat objects)
@@ -123 +77 @@
 
             psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
-                     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
+                     modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy);
             
             // if min point is too deviant, keep the old value
             if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
-                modelFLT->params->data.F32[2] = min.x + ix;
-                modelFLT->params->data.F32[3] = min.y + iy;
+                modelEXT->params->data.F32[2] = min.x + ix;
+                modelEXT->params->data.F32[3] = min.y + iy;
             }
             psFree (bicube);
@@ -134 +88 @@
 
         // set PSF parameters for this model
-        pmModel *model = pmModelFromPSF (modelFLT, psf);
-        psFree (modelFLT);
+        pmModel *model = pmModelFromPSF (modelEXT, psf);
+        psFree (modelEXT);
 
         // save the original coords
@@ -185 +139 @@
 
         // diagonal element (auto-cross-product)
-        f = psphotCrossProduct (Mi, Mi);
+        f = pmSourceCrossProduct (Mi, Mi);
         psSparseMatrixElement (sparse, i, i, f);
 
         // find the image x model value
-        f = psphotCrossProduct (Fi, Mi);
+        f = pmSourceCrossProduct (Fi, Mi);
         psSparseVectorElement (sparse, i, f);
 
@@ -203 +157 @@
             
             // got an overlap; calculate cross-product and add to output array
-            f = psphotCrossProduct (Mi, Mj);
+            f = pmSourceCrossProduct (Mi, Mj);
             psSparseMatrixElement (sparse, j, i, f); 
         }
@@ -235 +189 @@
     return true;
 }
-
-
-float psphotCrossProduct (pmSource *Mi, pmSource *Mj) {
-
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-    float flux;
-
-    psImage *Pi = Mi->pixels;
-    psImage *Pj = Mj->pixels;
-
-    psImage *Ti = Mi->mask;
-    psImage *Tj = Mj->mask;
-
-    Xs = PS_MAX (Pi->col0, Pj->col0);
-    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
-    
-    Ys = PS_MAX (Pi->row0, Pj->row0);
-    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
-    
-    xIs = Xs - Pi->col0;
-    xJs = Xs - Pj->col0;
-    yIs = Ys - Pi->row0;
-    yJs = Ys - Pj->row0;
-
-    xIe = Xe - Pi->col0;
-    yIe = Ye - Pi->row0;
-
-    flux = 0;
-    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
-        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.U8[yi][xi]) continue;
-            if (Tj->data.U8[yj][xj]) continue;
-            flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
-        }
-    }
-    return (flux);
-}

trunk/psphot/src/psphotEvalFLT.c

@@ -1 +1 @@
 # include "psphot.h"
 
-bool psphotEvalFLT (pmSource *source, pmModel *model)
+bool psphotEvalEXT (pmSource *source, pmModel *model)
 { 
     int keep;
 
-    // do we actually have a valid FLT model?
+    // do we actually have a valid EXT model?
     if (model == NULL) {
         source->mode &= ~PM_SOURCE_FITTED;
@@ -22 +22 @@
       case PM_MODEL_OFFIMAGE:
       default:
-        psLogMsg ("psphot", 5, "GAL fail fit\n");
+        psLogMsg ("psphot", 5, "EXT fail fit\n");
         source->mode |= PM_SOURCE_FAIL;
         return false;
     }
 
-    // unless we prove otherwise, this object is a galaxy.
-    // XXX EAM : do we need to save the old type?
-    source->type = PM_SOURCE_GALAXY;
+    // unless we prove otherwise, this object is extended
+    source->type = PM_SOURCE_EXTENDED;
 
-    // the following source->mode information pertains to modelPSF:
-    source->mode |= PM_SOURCE_FLTMODEL;
+    // the following source->mode information pertains to modelEXT:
+    source->mode |= PM_SOURCE_EXTMODEL;
 
     // if the object has a fitted peak below 0, the fit did not converge cleanly
@@ -44 +43 @@
 
     // poor-quality fit; only keep if nothing else works...
-    psLogMsg ("psphot", 5, "GAL poor fit\n");
+    psLogMsg ("psphot", 5, "EXT poor fit\n");
     source->mode |= PM_SOURCE_POOR;
     return false;

trunk/psphot/src/psphotEvalPSF.c

@@ -19 +19 @@
 
 // saturated stars should fall outside (larger), but have peaks above SATURATION
-// galaxies should be larger, cosmic rays smaller
+// extended sources should be larger, cosmic rays smaller
 // we also reject objects with S/N too low or ChiSquare to high
 
@@ -138 +138 @@
     }
 
-    // object appears to be large, suspected galaxy
+    // object appears to be large, suspected extended source
     if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) {
-        source->type = PM_SOURCE_GALAXY;
+        source->type = PM_SOURCE_EXTENDED;
         return false;
     }

trunk/psphot/src/psphotFitGalaxies.c

@@ -1 +1 @@
 # include "psphot.h"
 
-bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
+bool psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats)
 { 
     bool  status;
@@ -12 +12 @@
     psTimerStart ("psphot");
 
-    float GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
-    float GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
-    char       *modelName  = psMetadataLookupPtr (&status, config, "GAL_MODEL");
+    float EXT_MIN_SN       = psMetadataLookupF32 (&status, config, "EXT_MIN_SN");
+    float EXT_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS");
+    char       *modelName  = psMetadataLookupPtr (&status, config, "EXT_MODEL");
     pmModelType modelType  = pmModelSetType (modelName);
 
-    psphotInitRadiusFLT (config, skyStats, modelType);
+    psphotInitRadiusEXT (config, skyStats, modelType);
 
     for (int i = 0; i < sources->n; i++) {
@@ -23 +23 @@
 
         // only fit suspected extended sources
-        if (source->type != PM_SOURCE_GALAXY) continue;
+        if (source->type != PM_SOURCE_EXTENDED) continue;
         if (source->mode  & PM_SOURCE_BLEND) continue; 
 
         // skip possible extended sources below fit threshold
-        if (source->moments->SN < GAL_MIN_SN) continue;
+        if (source->moments->SN < EXT_MIN_SN) continue;
 
-        // recalculate the source moments using the larger galaxy moments radius
-        status = pmSourceMoments (source, GAL_MOMENTS_RAD);
+        // recalculate the source moments using the larger extended-source moments radius
+        status = pmSourceMoments (source, EXT_MOMENTS_RAD);
         if (!status) {
           source->mode |= PM_SOURCE_FAIL;  // better choice?
@@ -38 +38 @@
         // use the source moments, etc to guess basic model parameters
         pmModel  *model  = pmSourceModelGuess (source, modelType); 
-        source->modelFLT = model;
+        source->modelEXT = model;
         x = model->params->data.F32[2];
         y = model->params->data.F32[3];
 
         // sets the model radius (via source->model) and adjusts pixels as needed
-        psphotCheckRadiusFLT (imdata, source, model);
+        psphotCheckRadiusEXT (imdata, source, model);
 
-        // fit FLT (not PSF) model (set/unset the pixel mask)
+        // fit EXT (not PSF) model (set/unset the pixel mask)
         psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
         pmSourceFitModel (source, model, false);
@@ -61 +61 @@
         }
     }
-    psLogMsg ("psphot", 3, "fit GAL models: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
-    psLogMsg ("psphot", 4, "subtracted %d GAL models\n", Nsub);
+    psLogMsg ("psphot", 3, "fit EXT models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 4, "subtracted %d EXT objects\n", Nsub);
     return (true);
 }

trunk/psphot/src/psphotFullFit.c

@@ -14 +14 @@
     sources = psArraySort (sources, psphotSortBySN);
     
-    // galaxy model parameters
-    float GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
-    float GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+    // extended source model parameters
+    float EXT_MIN_SN       = psMetadataLookupF32 (&status, config, "EXT_MIN_SN");
+    float EXT_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS");
 
-    // 'galaxy' model descriptions
-    char         *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL");
-    pmModelType   modelTypeFLT = pmModelSetType (modelNameFLT);
+    // extended source model descriptions
+    char         *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL");
+    pmModelType   modelTypeEXT = pmModelSetType (modelNameEXT);
 
     psphotInitLimitsPSF (config);
     psphotInitRadiusPSF (config, sky, psf->type);
-    psphotInitRadiusFLT (config, sky, modelTypeFLT);
+    psphotInitRadiusEXT (config, sky, modelTypeEXT);
 
     for (int i = 0; i < sources->n; i++) {
@@ -82 +82 @@
 
         // skip the source if we don't think it is extended
-        if (source->type != PM_SOURCE_GALAXY) goto subLINEAR;
-        if (source->moments->SN < GAL_MIN_SN) goto subLINEAR;
+        if (source->type != PM_SOURCE_EXTENDED) goto subLINEAR;
+        if (source->moments->SN < EXT_MIN_SN) goto subLINEAR;
 
         // add the double-PSF fit mode
         // how do we compare the results?  chisq? 
         
-        // recalculate the source moments using the larger galaxy moments radius
-        if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) goto subLINEAR;
+        // recalculate the source moments using the larger extended-source moments radius
+        if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) goto subLINEAR;
 
         // use the source moments, etc to guess basic model parameters
-        source->modelFLT = pmSourceModelGuess (source, modelTypeFLT); 
-        pmModel *modelFLT = source->modelFLT;
+        source->modelEXT = pmSourceModelGuess (source, modelTypeEXT); 
+        pmModel *modelEXT = source->modelEXT;
 
-        psphotCheckRadiusFLT (imdata, source, modelFLT);
+        psphotCheckRadiusEXT (imdata, source, modelEXT);
 
-        x = modelFLT->params->data.F32[2];
-        y = modelFLT->params->data.F32[3];
+        x = modelEXT->params->data.F32[2];
+        y = modelEXT->params->data.F32[3];
 
-        // fit FLT (not PSF) model (set/unset the pixel mask)
-        psImageKeepCircle (source->mask, x, y, modelFLT->radius, "OR", PSPHOT_MASK_MARKED);
-        pmSourceFitModel (source, modelFLT, false);
-        psImageKeepCircle (source->mask, x, y, modelFLT->radius, "AND", ~PSPHOT_MASK_MARKED);
+        // fit EXT (not PSF) model (set/unset the pixel mask)
+        psImageKeepCircle (source->mask, x, y, modelEXT->radius, "OR", PSPHOT_MASK_MARKED);
+        pmSourceFitModel (source, modelEXT, false);
+        psImageKeepCircle (source->mask, x, y, modelEXT->radius, "AND", ~PSPHOT_MASK_MARKED);
 
         // track model evaluations
-        Niter += modelFLT[0].nIter;
+        Niter += modelEXT[0].nIter;
         Nfit++;
 
-        // does the FLT model succeed?
-        if (psphotEvalFLT (source, source->modelFLT)) {
-            pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false);
+        // does the EXT model succeed?
+        if (psphotEvalEXT (source, source->modelEXT)) {
+            pmSourceSubModel (source->pixels, source->mask, source->modelEXT, false, false);
             source->mode |=  PM_SOURCE_SUBTRACTED;
             source->mode &= ~PM_SOURCE_TEMPSUB;

trunk/psphot/src/psphotMagnitudes.c

@@ -21 +21 @@
         break;
 
-      case PM_SOURCE_GALAXY:
-        model = source->modelFLT;
+      case PM_SOURCE_EXTENDED:
+        model = source->modelEXT;
         if (model == NULL) return NULL;
         radius = model->radius;

trunk/psphot/src/psphotModelTest.c

@@ -6 +6 @@
     bool status;
     int modelType;
+    unsigned int Nfail;
     float obsMag, fitMag, value;
     char name[64];
@@ -18 +19 @@
     char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE");
     if (!status) {
-        psfModeWord = psStringCopy ("FLT");
+        psfModeWord = psStringCopy ("EXT");
     }
     bool psfMode = !strcasecmp (psfModeWord, "PSF");
@@ -26 +27 @@
         char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE");
         if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
-        psf = psphotReadPSF (psfFile);
+        psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE);
+        psf = pmPSFfromMD (psfData);
         modelType = psf->type;
     } else {

trunk/psphot/src/psphotOutput.c

@@ -14 +14 @@
 
     if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
-    if (psfFile != NULL) psphotWritePSF (psf, psfFile);
     psphotSamplePSFs (config, psf, imdata->image);
 
+    if (psfFile != NULL) {
+        psMetadata *psfData = pmPSFtoMD (NULL, psf);
+        psMetadataConfigWrite (psfData, psfFile);
+        psFree (psfData);
+    }
     if (outputFile == NULL) {
         psLogMsg ("output", 3, "no data output file selected");
@@ -166 +170 @@
     float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
     if (!status) ZERO_POINT = 25.0;
-    char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE");
 
     // write necessary information to output header
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS",     sources->n);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT",  PS_DATA_F32 | PS_META_REPLACE, "zero point",          ZERO_POINT);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nPSFstars);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code",     PHOTCODE);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
+    psphotUpdateHeader (imdata->header, config);
+    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n);
 
     // create file, write-out header
@@ -227 +218 @@
         psLineAdd (line, "%2d ",   type);
         psLineAdd (line, "%3.1f ", lsky);
-        psLineAdd (line, "%6.3f ", 99.999);    // should be 'Mgal
+        psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal
         psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); 
-        psLineAdd (line, "%6.2f ", PAR[4]);  // should be 'FHWM x'
-        psLineAdd (line, "%6.2f ", PAR[5]);  // should be 'FHWM y'
-        psLineAdd (line, "%5.1f\n", 0);   // should be theta
+        psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x'
+        psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y'
+        psLineAdd (line, "%5.1f\n", 0); // should be theta
         fwrite (line->line, 1, line->Nline, f);
     }
@@ -256 +247 @@
     float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
     float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
-    char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE");
 
     // write necessary information to output header
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS",     sources->n);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT",  PS_DATA_F32 | PS_META_REPLACE, "zero point",          ZERO_POINT);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code",     PHOTCODE);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
-    psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
-
-    // set NAXIS to 0
-    mdi = psMetadataLookup (imdata->header, "NAXIS");
-    mdi->data.S32 = 0;
-    mdi->type = PS_DATA_S32;
+    psphotUpdateHeader (imdata->header, config);
 
     table = psArrayAlloc (sources->n);
@@ -299 +274 @@
         psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", type);
         psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP",  PS_DATA_F32, "", lsky);
-        psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0);    // should be 'Mgal'
+        psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0);
         psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X",  PS_DATA_F32, "", source->apMag + ZERO_POINT); 
-        psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y",  PS_DATA_F32, "", PAR[4]);  // should be 'FHWM x'
-        psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", PAR[5]);  // should be 'FHWM y'
-        psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_TYPE_S32, "", 0);   // should be theta
-        // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY",   PS_DATA_STRING, "", NULL);
+        psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y",  PS_DATA_F32, "", PAR[4]);
+        psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", PAR[5]);
+        psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_DATA_STRING, "", "0");
+        psMetadataAdd (row, PS_LIST_TAIL, "DUMMY",   PS_DATA_STRING, "", "123");
     
         psArrayAdd (table, 100, row);
@@ -310 +285 @@
     }
 
-    // psAbort ("psphotOutput", "FITS table output not finished");
-
-    // psImage *image = psImageAlloc (1, 1, PS_TYPE_F32);
-    // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE");
-    // psFitsWriteImage (fits, imdata->header, image, 0);
-
+    // write the correct number of sources in the header
+    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n);
+    psMetadataAdd (imdata->header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true);
+
+    // set NAXIS to 0 (we don't write out the data array)
+    mdi = psMetadataLookup (imdata->header, "NAXIS");
+    mdi->data.S32 = 0;
+    mdi->type = PS_DATA_S32;
+
+    // create the basic table header
     theader = psMetadataAlloc ();
     psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");
 
-    psImage *image;
-
-    image = psImageAlloc (512, 512, PS_DATA_F32);
-    for (int j = 0; j < 512; j++) {
-        for (int i = 0; i < 512; i++) {
-            image->data.F32[j][i] = i;
-        }
-    }
-
+    // write out the header and table
     psFits *fits = psFitsOpen (filename, "w");
-    // psFitsWriteHeader (imdata->header, fits);
-    // psFitsWriteTable (fits, theader, table);
-
-    psFitsWriteImage (fits, NULL, image, 0);
-    psFitsWriteTable (fits, NULL, table);
+    psFitsWriteHeader (imdata->header, fits);
+    psFitsWriteTable (fits, theader, table);
     psFitsClose (fits);
 
@@ -348 +316 @@
     pmModelWritePSFs (sources, config, name, psf);
 
-    sprintf (name, "%s.flt.dat", filename);
-    pmModelWriteFLTs (sources, name);
+    sprintf (name, "%s.ext.dat", filename);
+    pmModelWriteEXTs (sources, name);
 
     sprintf (name, "%s.nul.dat", filename);
@@ -417 +385 @@
 
 // dump the sources to an output file
-bool pmModelWriteFLTs (psArray *sources, char *filename) {
+bool pmModelWriteEXTs (psArray *sources, char *filename) {
 
     double dPos, dMag;
@@ -427 +395 @@
     f = fopen (filename, "w");
     if (f == NULL) {
-        psLogMsg ("pmModelWriteFLTs", 3, "can't open output file for moments%s\n", filename);
+        psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename);
         return false;
     }
@@ -435 +403 @@
         pmSource *source = (pmSource *) sources->data[i];
 
-        if (source->type != PM_SOURCE_GALAXY) continue;
+        if (source->type != PM_SOURCE_EXTENDED) continue;
         model = pmSourceMagnitudes (source, NULL, 0.0);
         if (model == NULL) continue;
@@ -491 +459 @@
         source = sources->data[i];
 
-        // skip these sources (in PSF or FLT)
+        // skip these sources (in PSF or EXT)
         if (source->type == PM_SOURCE_STAR) continue;
-        if (source->type == PM_SOURCE_GALAXY) continue;
+        if (source->type == PM_SOURCE_EXTENDED) continue;
             
         if (source->moments == NULL) {
@@ -580 +548 @@
         return (1);
 
-      case PM_SOURCE_GALAXY:
+      case PM_SOURCE_EXTENDED:
         return (2);
 
@@ -597 +565 @@
 }
 
-psImage *pmModelPSFatXY (psImage *image, pmModel *modelFLT, pmPSF *psf, int x, int y, int dx, int dy) {
+psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) {
 
     psRegion region = {x - dx, x + dx, y - dy, y + dy};
@@ -603 +571 @@
     psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
     psImageInit (sample, 0);
-    modelFLT->params->data.F32[2] = x;
-    modelFLT->params->data.F32[3] = y;
-    pmModel *modelPSF = pmModelFromPSF (modelFLT, psf);
+    modelEXT->params->data.F32[2] = x;
+    modelEXT->params->data.F32[3] = y;
+    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
     pmSourceAddModel (sample, NULL, modelPSF, false, false);
     psFree (modelPSF);
@@ -625 +593 @@
     if (output[0] == 0) return false;
 
-    pmModel *modelFLT = pmModelAlloc (psf->type);
-    modelFLT->params->data.F32[0] = 0;
-    modelFLT->params->data.F32[1] = 1;
+    pmModel *modelEXT = pmModelAlloc (psf->type);
+    modelEXT->params->data.F32[0] = 0;
+    modelEXT->params->data.F32[1] = 1;
 
     psFits *fits = psFitsOpen (output, "w");
 
-    sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
+    sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25);
     psFitsWriteImage (fits, NULL, sample, 0);
-    sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
+    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
     psFitsWriteImage (fits, NULL, sample, 0);
-    sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, 25, 25, 25);
+    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25);
     psFitsWriteImage (fits, NULL, sample, 0);
-    sample = pmModelPSFatXY (image, modelFLT, psf, 25, image->numRows - 25, 25, 25);
+    sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25);
     psFitsWriteImage (fits, NULL, sample, 0);
-    sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
+    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
     psFitsWriteImage (fits, NULL, sample, 0);
 
     psFitsClose (fits);
     return (TRUE);
-}
-
-psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder) {
-
-    bool status;
-    char keyword[80];
-
-    // get polynomial orders
-    // XXX add status failures tests
-    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
-    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
-
-    psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
-
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny);
-            poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
-            if (!status) poly->mask[nx][ny] = 1;
-
-            sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny);
-            poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
-        }
-    }
-    return (poly);
-}
-
-psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder) {
-
-    bool status;
-    char keyword[80];
-
-    // get polynomial orders
-    // XXX add status failures tests
-    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
-    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
-    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
-
-    psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, PS_POLYNOMIAL_ORD);
-
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            for (int nz = 0; nz < poly->nZ + 1; nz++) {
-                sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
-                poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
-                if (!status) poly->mask[nx][ny][nz] = 1;
-
-                sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
-                poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
-            }
-        }
-    }
-    return (poly);
-}
-
-psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder) {
-
-    bool status;
-    char keyword[80];
-
-    // get polynomial orders
-    // XXX add status failures tests
-    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
-    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
-    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
-    int nTorder = psMetadataLookupS32 (&status, folder, "NORDER_T");
-
-    psPolynomial4D *poly = psPolynomial4DAlloc (nXorder, nYorder, nZorder, nTorder, PS_POLYNOMIAL_ORD);
-
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            for (int nz = 0; nz < poly->nZ + 1; nz++) {
-                for (int nt = 0; nt < poly->nT + 1; nt++) {
-                    sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
-                    poly->coeff[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
-                    if (!status) poly->mask[nx][ny][nz][nt] = 1;
-                    
-                    sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
-                    poly->coeffErr[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
-                }
-            }
-        }
-    }
-    return (poly);
-}
-
-// XXX : these may need F64, or %g format for output
-bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
-
-    int Nbyte;
-    char tmp;
-    char *root;
-    va_list argp;  
-
-    va_start (argp, format);
-    Nbyte = vsnprintf (&tmp, 0, format, argp);
-    va_end (argp);
-
-    if (!Nbyte) return false;
-
-    va_start (argp, format);
-    root = (char *) psAlloc (Nbyte + 1);
-    memset (root, 0, Nbyte + 1);
-    vsnprintf (root, Nbyte + 1, format, argp);
-    va_end (argp);
-
-    psMetadata *folder = psMetadataAlloc ();
-    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 2D polynomial", folder);
-    psFree (root);
-
-    // specify the polynomial orders
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
-
-    // place polynomial entries on folder
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            if (poly->mask[nx][ny]) continue;
-            psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny);
-            psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny);
-        }
-    }
-    return true;
-}
-
-bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...) {
-
-    int Nbyte;
-    char tmp;
-    char *root;
-    va_list argp;  
-
-    va_start (argp, format);
-    Nbyte = vsnprintf (&tmp, 0, format, argp);
-    va_end (argp);
-
-    if (!Nbyte) return false;
-
-    va_start (argp, format);
-    root = (char *) psAlloc (Nbyte + 1);
-    memset (root, 0, Nbyte + 1);
-    vsnprintf (root, Nbyte + 1, format, argp);
-    va_end (argp);
-
-    psMetadata *folder = psMetadataAlloc ();
-    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 3D polynomial", folder);
-    psFree (root);
-
-    // specify the polynomial orders
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
-
-    // place polynomial entries on folder
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            for (int nz = 0; nz < poly->nZ + 1; nz++) {
-                if (poly->mask[nx][ny][nz]) continue;
-                psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz);
-                psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz);
-            }
-        }
-    }
-    return true;
-}
-
-bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...) {
-
-    int Nbyte;
-    char tmp;
-    char *root;
-    va_list argp;  
-
-    va_start (argp, format);
-    Nbyte = vsnprintf (&tmp, 0, format, argp);
-    va_end (argp);
-
-    if (!Nbyte) return false;
-
-    va_start (argp, format);
-    root = (char *) psAlloc (Nbyte + 1);
-    memset (root, 0, Nbyte + 1);
-    vsnprintf (root, Nbyte + 1, format, argp);
-    va_end (argp);
-
-    psMetadata *folder = psMetadataAlloc ();
-    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 4D polynomial", folder);
-    psFree (root);
-
-    // specify the polynomial orders
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
-    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_T", PS_DATA_S32, "number of z orders", poly->nT);
-
-    // place polynomial entries on folder
-    for (int nx = 0; nx < poly->nX + 1; nx++) {
-        for (int ny = 0; ny < poly->nY + 1; ny++) {
-            for (int nz = 0; nz < poly->nZ + 1; nz++) {
-                for (int nt = 0; nt < poly->nT + 1; nt++) {
-                    if (poly->mask[nx][ny][nz][nt]) continue;
-                    psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz][nt], nx, ny, nz, nt);
-                    psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz][nt], nx, ny, nz, nt);
-                }
-            }
-        }
-    }
-    return true;
-}
-
-bool psphotWritePSF (pmPSF *psf, char *filename) {
-
-    psMetadata *psfdata = psMetadataAlloc ();
-
-    char *modelName = pmModelGetType (psf->type);
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
-
-    int nPar = pmModelParameterCount (psf->type)    ;
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
-
-    for (int i = 0; i < nPar - 4; i++) {
-        psPolynomial2D *poly = psf->params->data[i];
-        psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
-    }
-    psPolynomial4DtoMD (psfdata, psf->ApTrend, "APTREND");
-
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
-
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
-    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
-    
-    psMetadataConfigWrite (psfdata, filename);
-    psFree (psfdata);
-
-    return true;
-}
-
-pmPSF *psphotReadPSF (char *filename) {
-
-    bool status;
-    unsigned int Nfail;
-    char keyword[80];
-
-    psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE);
-
-    char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME");
-    pmModelType type = pmModelSetType (modelName);
-
-    pmPSF *psf = pmPSFAlloc (type);
-
-    int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR");
-    if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count");
-
-    for (int i = 0; i < nPar - 4; i++) {
-        sprintf (keyword, "PSF_PAR%02d", i);
-        psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
-        psPolynomial2D *poly = psPolynomial2DfromMD (folder);
-        psFree (psf->params->data[i]);
-        psf->params->data[i] = poly;
-    }
-    sprintf (keyword, "APTREND");
-    psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
-    psPolynomial4D *poly = psPolynomial4DfromMD (folder);
-    psFree (psf->ApTrend);
-    psf->ApTrend = poly;
-
-    psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID");
-    psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID");
-    psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS");
-
-    psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ");
-    psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS");
-
-    psFree (psfdata);
-    return (psf);
 }
 
@@ -937 +628 @@
     return true;
 }
+
+bool psphotUpdateHeader (psMetadata *header, psMetadata *config) {
+
+    // write necessary information to output header
+    psMetadataItemTransfer (header, config, "NSTARS");
+    psMetadataItemTransfer (header, config, "ZERO_PT");
+    psMetadataItemTransfer (header, config, "APMIFIT");
+    psMetadataItemTransfer (header, config, "dAPMIFIT");
+    psMetadataItemTransfer (header, config, "SKYBIAS");
+    psMetadataItemTransfer (header, config, "PHOTCODE");
+    
+    // write necessary information to output header
+    psMetadataItemTransfer (header, config, "NPSFSTAR");
+    psMetadataItemTransfer (header, config, "SKYBIAS");
+    psMetadataItemTransfer (header, config, "SKYSAT");
+    psMetadataItemTransfer (header, config, "APMIFIT");
+    psMetadataItemTransfer (header, config, "DAPMIFIT");
+    psMetadataItemTransfer (header, config, "NAPMIFIT");
+    psMetadataItemTransfer (header, config, "APLOSS");
+
+    // XXX these need to be defined from elsewhere
+    psMetadataAdd (header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
+    psMetadataAdd (header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
+    psMetadataAdd (header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
+    psMetadataAdd (header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
+    psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
+
+    return true;
+}

trunk/psphot/src/psphotRadiusChecks.c

@@ -31 +31 @@
     }
 
-# if (1)
     bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
     return status;
-
-# else
-
-    if (model->radius > PSF_OUTER_RADIUS) {
-        // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
-        psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
-        return true;
-    } else {
-        fprintf (stderr, "skipping %f,%f\n", model->params->data.F32[2], model->params->data.F32[3]);
-    }
-    return false;
-# endif
 }
 
+static float EXT_OUTER_RADIUS;
+static float EXT_FIT_PADDING;
+static float EXT_FLUX_LIMIT;
+static pmModelRadius modelRadiusEXT;
 
-static float GAL_OUTER_RADIUS;
-static float GAL_FIT_PADDING;
-static float GAL_FLUX_LIMIT;
-static pmModelRadius modelRadiusGAL;
-
-bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type) {
+bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type) {
 
     bool status;
 
-    float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA");
-    GAL_OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "GAL_OUTER_RADIUS");
-    GAL_FIT_PADDING      = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING");
-    GAL_FLUX_LIMIT       = GAL_FIT_NSIGMA * sky->sampleStdev;
+    float EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "EXT_FIT_NSIGMA");
+    EXT_OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "EXT_OUTER_RADIUS");
+    EXT_FIT_PADDING      = psMetadataLookupF32 (&status, config, "EXT_FIT_PADDING");
+    EXT_FLUX_LIMIT       = EXT_FIT_NSIGMA * sky->sampleStdev;
 
     // this function specifies the radius at this the model hits the given flux
-    modelRadiusGAL       = pmModelRadius_GetFunction (type);
+    modelRadiusEXT       = pmModelRadius_GetFunction (type);
     return true;
 }
 
-bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model) {
+bool psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model) {
 
     // set the fit radius based on the object flux limit and the model
-    model->radius = modelRadiusGAL (model->params, GAL_FLUX_LIMIT) + GAL_FIT_PADDING;
-    if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
+    model->radius = modelRadiusEXT (model->params, EXT_FLUX_LIMIT) + EXT_FIT_PADDING;
+    if (isnan(model->radius)) psAbort (__func__, "error in radius");
 
-# if (1)
     // redefine the pixels if needed
     bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
     return status;
-
-# else
-
-    // check if we need to redefine the pixels
-    if (model->radius > GAL_OUTER_RADIUS) {
-        // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
-        // do I need to save the old pixel and mask data?
-        psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius);
-        return true;
-    }
-    return false;
-# endif
 }

trunk/psphot/src/psphotReplaceUnfit.c

@@ -20 +20 @@
     return true;
 }
-
-# if (0)
-if (source->mode &  PM_SOURCE_BLEND) continue;
-if (source->type == PM_SOURCE_SATURATED) continue;
-
-if (source->type == PM_SOURCE_DEFECT) {
-    if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
-    continue; 
-}
-
-if (source->type == PM_SOURCE_GALAXY) {
-    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
-    if (source->mode & PM_SOURCE_POOR) goto addPSF;
-    continue;
-}
-
-// need to skip the successful fits
-if (source->type == PM_SOURCE_STAR) {
-    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
-    if (source->mode & PM_SOURCE_POOR) goto addPSF;
-    continue;
-}
-# endif

trunk/psphot/src/psphotRoughClass.c

@@ -9 +9 @@
     psfClump = pmSourcePSFClump (sources, config);
 
-    // group into STAR, COSMIC, GALAXY, SATURATED, etc.
+    // group into STAR, COSMIC, EXTENDED, SATURATED, etc.
     pmSourceRoughClass (sources, config, psfClump);
 

trunk/psphot/src/psphotSourceFits.c

@@ -41 +41 @@
 
 
-static float GAL_MIN_SN;
-static float GAL_MOMENTS_RAD;
-static pmModelType modelTypeFLT;
-
-bool psphotInitLimitsFLT (psMetadata *config, psStats *sky) {
+static float EXT_MIN_SN;
+static float EXT_MOMENTS_RAD;
+static pmModelType modelTypeEXT;
+
+bool psphotInitLimitsEXT (psMetadata *config, psStats *sky) {
 
     bool status;
 
-    // galaxy model parameters
-    GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
-    GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
-
-    // 'galaxy' model descriptions
-    char *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL");
-    modelTypeFLT = pmModelSetType (modelNameFLT);
-    psphotInitRadiusFLT (config, sky, modelTypeFLT);
+    // extended source model parameters
+    EXT_MIN_SN       = psMetadataLookupF32 (&status, config, "EXT_MIN_SN");
+    EXT_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS");
+
+    // extended source model descriptions
+    char *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL");
+    modelTypeEXT = pmModelSetType (modelNameEXT);
+    psphotInitRadiusEXT (config, sky, modelTypeEXT);
 
     return true;
@@ -63 +63 @@
 bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources) { 
 
-    bool okFLT, okDBL;
-    float chiFLT, chiDBL;
+    bool okEXT, okDBL;
+    float chiEXT, chiDBL;
 
     // skip the source if we don't think it is extended
@@ -70 +70 @@
     if (source->type == PM_SOURCE_DEFECT) return false;
     if (source->type == PM_SOURCE_SATURATED) return false;
-    if (source->moments->SN < GAL_MIN_SN) return false;
-
-    // recalculate the source moments using the larger galaxy moments radius
-    if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) return false;
+    if (source->moments->SN < EXT_MIN_SN) return false;
+
+    // recalculate the source moments using the larger extended-source moments radius
+    if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false;
 
     psTrace ("psphot.blend", 5, "trying blob...\n");
@@ -79 +79 @@
     pmSource *tmpSrc = pmSourceAlloc ();
 
-    pmModel *FLT = psphotFitFLT (imdata, source);
-    okFLT = psphotEvalFLT (tmpSrc, FLT);
-    chiFLT = FLT->chisq / FLT->nDOF;
+    pmModel *EXT = psphotFitEXT (imdata, source);
+    okEXT = psphotEvalEXT (tmpSrc, EXT);
+    chiEXT = EXT->chisq / EXT->nDOF;
 
     psArray *DBL = psphotFitDBL (imdata, source);
@@ -89 +89 @@
     chiDBL = ONE->chisq / ONE->nDOF;
 
-    if (okFLT && okDBL) {
-        psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiFLT, chiDBL);
+    if (okEXT && okDBL) {
+        psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
         // XXX EAM : a bogus bias: need to examine this better
-        if (3*chiFLT > chiDBL) goto keepDBL;
-        goto keepFLT;
-    }
-
-    if (okFLT && !okDBL) goto keepFLT;
-    if (!okFLT && okDBL) goto keepDBL;
+        if (3*chiEXT > chiDBL) goto keepDBL;
+        goto keepEXT;
+    }
+
+    if (okEXT && !okDBL) goto keepEXT;
+    if (!okEXT && okDBL) goto keepDBL;
     return false;
 
-keepFLT:
-    // sub FLT
-    pmSourceSubModel (source->pixels, source->mask, FLT, false, false);
-    psTrace ("psphot.blend", 5, "blob as FLT: %f %f\n", FLT->params->data.F32[2], FLT->params->data.F32[3]);
+keepEXT:
+    // sub EXT
+    pmSourceSubModel (source->pixels, source->mask, EXT, false, false);
+    psTrace ("psphot.blend", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[2], EXT->params->data.F32[3]);
 
     // save new model
-    source->modelFLT = FLT;
+    source->modelEXT = EXT;
     source->mode |=  PM_SOURCE_SUBTRACTED;
     source->mode &= ~PM_SOURCE_TEMPSUB;
@@ -168 +168 @@
     y = PSF->params->data.F32[3];
 
-    // fit FLT (not PSF) model (set/unset the pixel mask)
+    // fit EXT (not PSF) model (set/unset the pixel mask)
     psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED);
     pmSourceFitSet (source, modelSet, true);
@@ -176 +176 @@
 }
 
-pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source) { 
+pmModel *psphotFitEXT (eamReadout *imdata, pmSource *source) { 
 
     float x, y;
 
     // use the source moments, etc to guess basic model parameters
-    pmModel *FLT = pmSourceModelGuess (source, modelTypeFLT); 
-
-    psphotCheckRadiusFLT (imdata, source, FLT);
-
-    x = FLT->params->data.F32[2];
-    y = FLT->params->data.F32[3];
-
-    // fit FLT (not PSF) model (set/unset the pixel mask)
-    psImageKeepCircle (source->mask, x, y, FLT->radius, "OR", PSPHOT_MASK_MARKED);
-    pmSourceFitModel (source, FLT, false);
-    psImageKeepCircle (source->mask, x, y, FLT->radius, "AND", ~PSPHOT_MASK_MARKED);
-
-    return (FLT);
+    pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); 
+
+    psphotCheckRadiusEXT (imdata, source, EXT);
+
+    x = EXT->params->data.F32[2];
+    y = EXT->params->data.F32[3];
+
+    // fit EXT (not PSF) model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, EXT->radius, "OR", PSPHOT_MASK_MARKED);
+    pmSourceFitModel (source, EXT, false);
+    psImageKeepCircle (source->mask, x, y, EXT->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+    return (EXT);
 }
 

trunk/psphot/src/psphotTest.c

@@ -1 +1 @@
 # include "psphot.h"
-bool apResidSetPower (float value);
 
 void psExit (int status, char *process, char *format, ...) {
@@ -16 +15 @@
 int main (int argc, char **argv) {
 
-    // read in a file consisting of:
-    // x y radius fitmag apmag dmag
+    psMetadata *row;
+    psArray *table;
 
-    double fit;
-    double ap;
+    psMetadataItem *mdi;
+    psRegion region = {0,0,0,0};        // a region representing the entire array
 
     psphotTestArguments (&argc, argv);
 
-    FILE *f = fopen (argv[1], "r");
-    if (f == NULL) psExit (2, "test", "no file\n");
+    psFits *file = psFitsOpen (argv[1], "r");
+    psMetadata *header = psFitsReadHeader (NULL, file);
+    psImage *image = psFitsReadImage (NULL, file, region, 0);
+    psFitsClose (file);
 
-    float power = atof (argv[2]);
-    apResidSetPower (power);
+    psMetadataConfigWrite (header, argv[2]);
 
-    psVector *x     = psVectorAlloc (100, PS_TYPE_F64);
-    psVector *y     = psVectorAlloc (100, PS_TYPE_F64);
-    psVector *rad   = psVectorAlloc (100, PS_TYPE_F64);
-    psVector *rflux = psVectorAlloc (100, PS_TYPE_F64);
-    psVector *apres = psVectorAlloc (100, PS_TYPE_F64);
-    psVector *dmag  = psVectorAlloc (100, PS_TYPE_F64);
+    // attempt to write image with NAXIS = 0
+    mdi = psMetadataLookup (header, "NAXIS");
+    mdi->data.S32 = 0;
+    mdi->type = PS_DATA_S32;
+    
+    // create a test image
+    // psImage *tmpimage = psImageAlloc (10, 10, PS_DATA_F32);
 
-    int Npt = 0;
-    while (fscanf (f, "%lf %lf %lf %lf %lf %lf", 
-                   &x->data.F64[Npt], 
-                   &y->data.F64[Npt], 
-                   &rad->data.F64[Npt], 
-                   &fit, &ap, 
-                   &dmag->data.F64[Npt]) != EOF) {
+    // create a test table
+    table = psArrayAlloc (10);
+    table->n = 0;
 
-        if (fabs(ap-fit) > 0.1) continue;
-        if (fit > -10) continue;
+    for (int i = 0; i < 10; i++) {
+        row = psMetadataAlloc ();
+        psMetadataAdd (row, PS_LIST_TAIL, "ROW",   PS_DATA_S32,    "", i);
+        psMetadataAdd (row, PS_LIST_TAIL, "FROW",  PS_TYPE_F32,    "", 0.1*i);
+        psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line");
+    
+        table->data[i] = row;
+    }
+    table->n = 10;
+    
+    psMetadata *theader = psMetadataAlloc ();
+    psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");
 
-        apres->data.F64[Npt] = ap-fit;
-        rflux->data.F64[Npt] = pow(10.0, 0.4*fit);
+    psFits *fits = psFitsOpen (argv[3], "w");
+    // psFitsWriteImage (fits, header, tmpimage, 0);
+    psFitsWriteHeader (header, fits);
+    psFitsWriteTable (fits, theader, table);
+    psFitsClose (fits);
 
-        psVectorExtend (x,     100, 1);
-        psVectorExtend (y,     100, 1);
-        psVectorExtend (rad,   100, 1);
-        psVectorExtend (rflux, 100, 1);
-        psVectorExtend (apres, 100, 1);
-        psVectorExtend (dmag,  100, 1);
-        Npt ++;
-    }
-    x->n = y->n = rad->n = rflux->n = apres->n = dmag->n = Npt;
-
-    psphotApResidPolyFit (x, y, rad, rflux, apres, dmag);
-
+    psFree (header);
+    psFree (image);
     exit (0);
 }

trunk/psphot/src/psphotTestArguments.c

@@ -8 +8 @@
   psArgumentVerbosity (argc, argv);
 
-  if (*argc != 3) usage ();
+  if (*argc != 4) usage ();
 
   return;
@@ -15 +15 @@
 static int usage () {
 
-    fprintf (stderr, "USAGE: psphotModelTest (apresid.dat) (power)\n");
+    fprintf (stderr, "USAGE: psphotTest (input.fits) (output.hdr) (out.fits)\n");
     exit (2);
 }