Index: trunk/psphot/src/psphotOutput.c
===================================================================
--- trunk/psphot/src/psphotOutput.c	(revision 5802)
+++ trunk/psphot/src/psphotOutput.c	(revision 5828)
@@ -5,10 +5,14 @@
 
     bool status;
+
+    psTimerStart ("psphot");
 
     char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE");
     char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE");
     char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE");
+    char *psfFile    = psMetadataLookupPtr (&status, config, "PSF_OUTPUT_FILE");
 
     if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
+    if (psfFile != NULL) psphotWritePSF (psf, psfFile);
 
     if (outputFile == NULL) {
@@ -44,4 +48,5 @@
 	return;
     }
+
     psAbort ("psphot", "unknown output mode %s", outputMode);
 }
@@ -287,5 +292,4 @@
 
     // set NAXIS to 0
-    fprintf (stderr, "setting naxis\n"); 
     mdi = psMetadataLookup (imdata->header, "NAXIS");
     mdi->data.S32 = 0;
@@ -540,5 +544,4 @@
     }
 
-    fprintf (stderr, "writing out moments\n");
     for (i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -604,8 +607,16 @@
 
 
-bool psphotSamplePSFs (pmPSF *psf, psImage *image) {
-
-  // make sample PSFs for 4 corners and the center
+bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image) {
+
+    bool status;
+
+    // make sample PSFs for 4 corners and the center
     psImage *sample;
+
+    // optional dump of all rough source data
+    char *output = psMetadataLookupPtr (&status, config, "PSF_SAMPLE_FILE");
+    if (!status) return false;
+    if (output == NULL) return false;
+    if (output[0] == 0) return false;
 
     pmModel *modelFLT = pmModelAlloc (psf->type);
@@ -613,5 +624,5 @@
     modelFLT->params->data.F32[1] = 1;
 
-    psFits *fits = psFitsOpen ("sample.psf.fits", "w");
+    psFits *fits = psFitsOpen (output, "w");
 
     sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
@@ -630,4 +641,206 @@
 }
 
+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);
+}
+
+// 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 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);
+    }
+    psPolynomial3DtoMD (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);
+    psPolynomial3D *poly = psPolynomial3DfromMD (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);
+}
+
 bool psphotDumpMoments (psMetadata *config, psArray *sources) {