Changeset 21536 for branches/cnb_branch_20090215/psModules

Timestamp:

Feb 19, 2009, 7:59:50 AM (17 years ago)

Author:

beaumont

Message:

Added visualizations to ppSub. Set up a single variable in pmVisual to control when plots are drawn.

Location:

branches/cnb_branch_20090215/psModules/src

Files:

• astrom/pmAstrometryVisual.c (modified) (27 diffs)
• astrom/pmAstrometryVisual.h (modified) (1 diff)
• camera/pmFPACopy.c (modified) (1 diff)
• detrend/pmShutterCorrection.h (modified) (1 diff)
• extras/pmVisual.c (modified) (5 diffs)
• extras/pmVisual.h (modified) (3 diffs)
• imcombine/Makefile.am (modified) (2 diffs)
• imcombine/pmPSFEnvelope.c (modified) (6 diffs)
• imcombine/pmStack.c (modified) (2 diffs)
• imcombine/pmSubtractionIO.c (modified) (3 diffs)
• imcombine/pmSubtractionVisual.c (modified) (10 diffs)
• imcombine/pmSubtractionVisual.h (modified) (1 diff)
• objects/Makefile.am (modified) (1 diff)
• objects/pmModel.c (modified) (5 diffs)
• objects/pmModel.h (modified) (4 diffs)
• objects/pmPeaks.c (modified) (7 diffs)
• objects/pmSource.c (modified) (7 diffs)
• objects/pmSource.h (modified) (5 diffs)
• objects/pmSourceFitModel.c (modified) (2 diffs)
• objects/pmSourceIO.c (modified) (5 diffs)
• objects/pmSourceIO.h (modified) (4 diffs)
• objects/pmSourceIO_CMF_PS1_V1.c (added)
• objects/pmSourceIO_PS1_CAL_0.c (modified) (19 diffs)
• objects/pmSourceIO_PS1_DEV_1.c (modified) (12 diffs)
• objects/pmSourcePhotometry.c (modified) (11 diffs)
• objects/pmSourcePhotometry.h (modified) (2 diffs)
• objects/pmSourceVisual.c (modified) (6 diffs)
• objects/pmSourceVisual.h (modified) (2 diffs)

branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.c

@@ -23 +23 @@
 
 //variables to determine when things are plotted
-static bool isVisual             = false;
 static bool plotGridMatch        = true;
 static bool plotTweak            = true;
@@ -48 +47 @@
 /* Initialization Routines  */
 
-bool pmAstromSetVisual (bool mode) {
-    isVisual = mode;
-    return true;
-}
-
-
 bool pmAstromVisualClose()
 {
@@ -68 +61 @@
 {
     // make sure we want to plot this
-    if (!plotRawStars || !isVisual) return true;
+    if (!plotRawStars || !pmVisualIsVisual()) return true;
 
     //set up plot region
     if (!pmVisualInitWindow (&kapa, "psastro:plots")){
-        isVisual = false;
         return false;
     }
@@ -222 +214 @@
 
     // pause and wait for user input:
-    pmVisualAskUser(&plotRawStars, &isVisual);
+    pmVisualAskUser(&plotRawStars);
 
     psFree (xVec);
@@ -234 +226 @@
 {
     //make sure we want to plot this
-    if (!isVisual || !plotRefStars) return true;
+    if (!pmVisualIsVisual() || !plotRefStars) return true;
 
     //set up plotting variables
     if (!pmVisualInitWindow (&kapa, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -284 +275 @@
 
     // pause and wait for user input:
-    pmVisualAskUser(&plotRefStars, &isVisual);
+    pmVisualAskUser(&plotRefStars);
 
     psFree (xVec);
@@ -301 +292 @@
 
     // make sure we want to plot this
-    if ( !isVisual || !plotLumFunc ) return true;
+    if ( !pmVisualIsVisual() || !plotLumFunc ) return true;
 
     //set up plotting variables
     if ( !pmVisualInitWindow (&kapa, "psastro:plots")){
-        isVisual = false;
         return false;
     }
@@ -377 +367 @@
 
         // pause and wait for user input:
-        pmVisualAskUser (&plotLumFunc, &isVisual);
+        pmVisualAskUser (&plotLumFunc);
     }
     return true;
@@ -391 +381 @@
 
     //make sure we want to plot this
-    if (!isVisual || !plotRemoveClumps) return true;
+    if (!pmVisualIsVisual() || !plotRemoveClumps) return true;
 
     //set up plot variables
     if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -480 +469 @@
 
     //ask for user input and finish
-    pmVisualAskUser (&plotRemoveClumps, &isVisual);
+    pmVisualAskUser (&plotRemoveClumps);
     psFree (xVec);
     psFree (yVec);
@@ -496 +485 @@
 
     //make sure we want to plot this
-    if (!isVisual || !plotOneChipFit)
+    if (!pmVisualIsVisual() || !plotOneChipFit)
         return true;
 
     if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -511 +499 @@
 
     //ask for user input and finish
-    pmVisualAskUser(&plotOneChipFit, &isVisual);
+    pmVisualAskUser(&plotOneChipFit);
     return true;
 }
@@ -522 +510 @@
 {
     //make sure we want to plot this
-    if(!isVisual || !plotFixChips) return true;
+    if(!pmVisualIsVisual() || !plotFixChips) return true;
 
     if(!pmVisualInitWindow(&kapa, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -588 +575 @@
     KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y");
 
-    pmVisualAskUser(&plotFixChips, &isVisual);
+    pmVisualAskUser(&plotFixChips);
     psFree(xNew);
     psFree(yNew);
@@ -606 +593 @@
 
     //make sure we want to plot this
-    if (!isVisual || !plotAstromGuessCheck) return true;
+    if (!pmVisualIsVisual() || !plotAstromGuessCheck) return true;
 
     //set up graph window
     if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -686 +672 @@
     psFree(xResid);
     psFree(yResid);
-    pmVisualAskUser (&plotAstromGuessCheck, &isVisual);
+    pmVisualAskUser (&plotAstromGuessCheck);
     return true;
 }
@@ -696 +682 @@
 {
     //make sure we want to plot this
-    if (!isVisual || !plotCommonScale) return true;
+    if (!pmVisualIsVisual() || !plotCommonScale) return true;
 
     if (!pmVisualInitWindow(&kapa, "psastro:plots")){
-        isVisual = false;
         return false;
     }
@@ -736 +721 @@
     KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP);
 
-    pmVisualAskUser (&plotCommonScale, &isVisual);
+    pmVisualAskUser (&plotCommonScale);
 
     psFree(xVec);
@@ -750 +735 @@
 
     //make sure we want to plot this
-    if (!isVisual || !plotMosaicOneChip) return true;
+    if (!pmVisualIsVisual() || !plotMosaicOneChip) return true;
 
     if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -759 +743 @@
     //plot the residuals
     if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) {
-        isVisual = false;
+        pmVisualSetVisual(false);
         return false;
     }
 
     //ask for user input and finish
-    pmVisualAskUser(&plotMosaicOneChip, &isVisual);
+    pmVisualAskUser(&plotMosaicOneChip);
 
     return true;
@@ -775 +759 @@
 {
     //make sure we want to plot this
-    if (!isVisual || !plotMosaicMatches) return true;
+    if (!pmVisualIsVisual() || !plotMosaicMatches) return true;
 
     char title[60];
@@ -782 +766 @@
 
     if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
 
     if (!residPlot(rawstars, refstars, match, recipe, title)){
-        isVisual = false;
+        pmVisualSetVisual(false);
         return false;
     }
 
     //ask for user input
-    pmVisualAskUser (&plotMosaicMatches, &isVisual);
+    pmVisualAskUser (&plotMosaicMatches);
     return true;
 }
@@ -807 +790 @@
 {
     //make sure we want to plot this
-    if (!isVisual || !plotGridMatch) return true;
+    if (!pmVisualIsVisual() || !plotGridMatch) return true;
     if (!pmVisualInitWindow(&kapa, "psastro:plots")){
-        isVisual = false;
         return false;
     }
@@ -952 +934 @@
     KapaPlotVector (kapa, 2, yslice, "y");
 
-    pmVisualAskUser(&plotGridMatch, &isVisual);
+    pmVisualAskUser(&plotGridMatch);
     psFree(dXplot);
     psFree(dYplot);
@@ -966 +948 @@
 {
     //make sure we want to plot this
-    if (!isVisual || !plotTweak) return true;
+    if (!pmVisualIsVisual() || !plotTweak) return true;
     if (!pmVisualInitWindow(&kapa, "psastro:plots")) {
-        isVisual = false;
         return false;
     }
@@ -1051 +1032 @@
                    KAPA_LABEL_XP);
 
-    pmVisualAskUser(&plotTweak, &isVisual);
+    pmVisualAskUser(&plotTweak);
 
     psFree(xIndices);
@@ -1341 +1322 @@
 # else
 
-bool pmAstromSetVisual(bool mode) { return true; }
 bool pmAstromVisualClose() { return true; }
 bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset) { return true; }

branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.h

@@ -22 +22 @@
     int sPeak;                          ///< sum of stars to peak bin
 } pmLumFunc;
-
-
-/** Enable or disable visual plotting for psastro routines
- * @param mode true/false to enable/disable plotting
- * @return true for success */
-bool pmAstromSetVisual(bool mode);
 
 

branches/cnb_branch_20090215/psModules/src/camera/pmFPACopy.c

@@ -457 +457 @@
     psMetadataCopy(target->parent->concepts, source->parent->concepts);
 
+    // Update the astrometric parameters
+    target->toFPA = psMemIncrRefCounter (source->toFPA);
+    target->fromFPA = psMemIncrRefCounter (source->fromFPA);
+
+    // Update the parent fpa astrometry parameters, or check that they match
+    pmFPA *targetFPA = target->parent;
+    pmFPA *sourceFPA = source->parent;
+    // XXX should we require that both of these exist?
+
+    if (targetFPA && sourceFPA) {
+        if (targetFPA->toSky) {
+            psAssert (targetFPA->toSky == sourceFPA->toSky, "chips within FPA have inconsistent astrometry references");
+        } else {
+            targetFPA->toSky = psMemIncrRefCounter (sourceFPA->toSky);
+        }
+        if (targetFPA->toTPA) {
+            psAssert (targetFPA->toTPA == sourceFPA->toTPA, "chips within FPA have inconsistent astrometry references");
+        } else {
+            targetFPA->toTPA = psMemIncrRefCounter (sourceFPA->toTPA);
+        }
+        if (targetFPA->fromTPA) {
+            psAssert (targetFPA->fromTPA == sourceFPA->fromTPA, "chips within FPA have inconsistent astrometry references");
+        } else {
+            targetFPA->fromTPA = psMemIncrRefCounter (sourceFPA->fromTPA);
+        }
+    }
+
     target->data_exists = true;
     return status;

branches/cnb_branch_20090215/psModules/src/detrend/pmShutterCorrection.h

@@ -5 +5 @@
  * @author Paul Price, IfA
  *
- * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-12 19:33:30 $
+ * @version $Revision: 1.23.2.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  * Copyright 2006 Institute for Astronomy, University of Hawaii
  */

branches/cnb_branch_20090215/psModules/src/extras/pmVisual.c

@@ -28 +28 @@
 # define KAPAY 700
 
+static bool isVisual = false;
+
+
+bool pmVisualSetVisual(bool value) {
+    isVisual = value;
+    return true;
+}
+
+
+bool pmVisualIsVisual(void) {return isVisual;}
+
+bool pmVisualClose(void) {
+    pmAstrometryVisualClose();
+    pmSubtractionVisualClose();
+    pmStackVisualClose();
+    //XXX handle psphot
+    //    psphotVisualClose();
+    return true;
+}
 
 bool pmVisualInitWindow (int *kapid, char *name) {
@@ -34 +53 @@
         if (*kapid == -1) {
             fprintf (stderr, "Failure to open kapa.\n");
+            isVisual = false;
             return false;
         }
@@ -52 +72 @@
 
 
-bool pmVisualAskUser(bool *plotFlag, bool *packageFlag)
+bool pmVisualAskUser(bool *plotFlag)
 {
     char key[10];
@@ -63 +83 @@
     }
     if (key[0] == 'a') {
-        *packageFlag = false;
+        isVisual = false;
     }
     return true;
@@ -304 +324 @@
 
 #else
-
+bool pmVisualSetVisual(bool value) {return true;}
+bool pmVisualIsVisual(void) {return false;}
+bool pmVisualClose(void) {return true;}
 bool pmVisualInitWindow(int *kapid, char *name){return true;}
 bool pmVisualInitGraph (int kapa, void *section, void *graphdata){return true;}
 
-bool pmVisualAskUser(bool *plotFlag, bool *packageFlag){return true;}
+bool pmVisualAskUser(bool *plotFlag){return true;}
 bool pmVisualImStats(psImage *image, double *mean,
                      double *stdev, double *min, double *max){return true;}

branches/cnb_branch_20090215/psModules/src/extras/pmVisual.h

@@ -10 +10 @@
 
 #if (HAVE_KAPA)
+
+
+/** Globally enable or disable plotting
+ * @param value - true to enable plotting
+ * @return true
+ */
+bool pmVisualSetVisual(bool value);
+
+
+/** Check whether plotting is enabled
+ * @retrun true if plots should be generated
+ */
+bool pmVisualIsVisual(void);
+
+
+/** Destroy plotting windows at the end of a run
+ * @return true for success */
+bool pmVisualClose(void);
 
 
@@ -32 +50 @@
  * At the user's request, this will disable diagnostic plotting.
  * @param plotFlag, set to false if this plot should be disabled in the future
- * @param packageFlag, set to false if all plots from this package (e.g. psastro, pmSubtraction)
- *                      should be disabled.
  */
-bool pmVisualAskUser(bool *plotFlag, bool *packageFlag);
+bool pmVisualAskUser(bool *plotFlag);
 
 
@@ -115 +131 @@
 
 // kapa-specific data types are changed to void
+bool pmVisualSetVisual(bool value);
+bool pmVisualIsVisual(void);
+bool pmVisualClose(void);
 bool pmVisualInitWindow (int *kapid, char *name);
 bool pmVisualInitGraph (int kapa, void *section, void *graphdata);
-bool pmVisualAskUser(bool *plotFlag, bool *packageFlag);
+bool pmVisualAskUser(bool *plotFlag);
 bool pmVisualScaleImage(int kapaFD, psImage *inImage,
                         const char *name, int channel, bool clip);

branches/cnb_branch_20090215/psModules/src/imcombine/Makefile.am

@@ -19 +19 @@
         pmSubtractionThreads.c  \
         pmPSFEnvelope.c         \
-        pmSubtractionVisual.c
+        pmSubtractionVisual.c   \
+        pmStackVisual.c
 
 pkginclude_HEADERS = \
@@ -37 +38 @@
         pmSubtractionThreads.h  \
         pmPSFEnvelope.h         \
-        pmSubtractionVisual.h
+        pmSubtractionVisual.h   \
+        pmStackVisual.h
 
 CLEANFILES = *~

branches/cnb_branch_20090215/psModules/src/imcombine/pmPSFEnvelope.c

@@ -30 +30 @@
 
 #include "pmPSFEnvelope.h"
-
-
-
-//#define TESTING                         // Enable test output
+#include "pmStackVisual.h"
+
+
+#define TESTING                         // Enable test output
 #define PEAK_FLUX 1.0e4                 // Peak flux for each source
 #define SKY_VALUE 0.0e0                 // Sky value for fake image
@@ -162 +162 @@
             psStringAppend(&name, "psf_field_%03d.fits", i);
             psFits *fits = psFitsOpen(name, "w");
+            pmStackVisualPlotTestImage(fakeRO->image, name);
             psFitsWriteImage(fits, NULL, fakeRO->image, 0, NULL);
             psFitsClose(fits);
@@ -179 +180 @@
         // Write out the envelope
         psFits *fits = psFitsOpen("psf_field_envelope.fits", "w");
+        pmStackVisualPlotTestImage(envelope, "psf_field_envelope.fits");
         psFitsWriteImage(fits, NULL, envelope, 0, NULL);
         psFitsClose(fits);
@@ -231 +233 @@
         // Write out the envelope
         psFits *fits = psFitsOpen("psf_field_full.fits", "w");
+        pmStackVisualPlotTestImage(readout->image, "psf_field_full.fits");
         psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
         psFitsClose(fits);
@@ -312 +315 @@
         {
             psFits *fits = psFitsOpen("psf_field_model.fits", "w");
+            pmStackVisualPlotTestImage(generated->image, "psf_field_model.fits");
             psFitsWriteImage(fits, NULL, generated->image, 0, NULL);
             psFitsClose(fits);
@@ -318 +322 @@
         {
             psFits *fits = psFitsOpen("psf_field_resid.fits", "w");
+            pmStackVisualPlotTestImage(generated->image, "psf_field_resid.fits");
             psFitsWriteImage(fits, NULL, generated->image, 0, NULL);
             psFitsClose(fits);

branches/cnb_branch_20090215/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-13 23:52:14 $
+ *  @version $Revision: 1.48.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -33 +33 @@
 #define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
 
-//#define TESTING                         // Enable test output
-//#define TEST_X 2318                     // x coordinate to examine
-//#define TEST_Y 2306                     // y coordinate to examine
+//CNB change all of these to false
+#define TESTING                         // Enable test output
+#define TEST_X 2318                     // x coordinate to examine
+#define TEST_Y 2306                     // y coordinate to examine
 
 

branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionIO.c

@@ -16 +16 @@
 
 #include "pmSubtractionIO.h"
+#include "pmStackVisual.h"
 
 #define ARRAY_BUFFER 16                 // Number to add to array at a time
@@ -144 +145 @@
     // CVS tags, used to identify the version of this file (in case incompatibilities are introduced)
     psString cvsFile = psStringCopy("$RCSfile: pmSubtractionIO.c,v $");
-    psString cvsRev  = psStringCopy("$Revision: 1.9 $");
-    psString cvsDate = psStringCopy("$Date: 2008-10-24 21:21:06 $");
+    psString cvsRev  = psStringCopy("$Revision: 1.9.18.1 $");
+    psString cvsDate = psStringCopy("$Date: 2009-02-19 17:59:50 $");
     psStringSubstitute(&cvsFile, NULL, "RCSfile: ");
     psStringSubstitute(&cvsRev,  NULL, "Revision: ");
@@ -167 +168 @@
 
     psImage *image = psMetadataLookupPtr(NULL, ro->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL_IMAGE); // Image
+    pmStackVisualPlotTestImage(image, "Subtraction_kernels.fits");
+
     if (image && !psFitsWriteImage(fits, header, image, 0, EXTNAME_IMAGE)) {
         psError(PS_ERR_IO, false, "Unable to write subtraction kernel image.");

branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.c

@@ -28 +28 @@
 
 //variables to determine when things are plotted
-static bool isVisual             = false;
 static bool plotConvKernels      = true;
 static bool plotStamps           = true;
@@ -42 +41 @@
 // Initialization Routines
 
-/** start or stop plotting */
-bool pmSubtractionSetVisual (bool mode) {
-    isVisual = mode;
-    return true;
-}
 
 
@@ -65 +59 @@
  *    @return true for success */
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {
-    if (!isVisual || !plotConvKernels) return true;
+    if (!pmVisualIsVisual() || !plotConvKernels) return true;
     if (!pmVisualInitWindow(&kapa, "ppSub:Images")) {
-        isVisual = false;
         return false;
     }
     pmVisualScaleImage(kapa, convKernels, "Convolution_Kernels", 0, true);
-    pmVisualAskUser(&plotConvKernels, &isVisual);
+    pmVisualAskUser(&plotConvKernels);
     return true;
 }
@@ -80 +73 @@
     @return true for success */
 bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro) {
-    if (!isVisual || !plotStamps) return true;
+    if (!pmVisualIsVisual() || !plotStamps) return true;
     if (!pmVisualInitWindow (&kapa, "ppSub:Images")) {
-        isVisual = false;
         return false;
     }
     if (!pmVisualInitWindow (&kapa2, "ppSub:StampMasterImage")) {
-        isVisual = false;
         return false;
     }
@@ -145 +136 @@
     pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true);
 
-    pmVisualAskUser(&plotStamps, &isVisual);
+    pmVisualAskUser(&plotStamps);
     return true;
 }
@@ -152 +143 @@
 bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {
 
-    if (!isVisual || !plotLeastSquares) return true;
+    if (!pmVisualIsVisual() || !plotLeastSquares) return true;
     if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
-        isVisual = false;
         return false;
     }
@@ -209 +199 @@
     pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true);
 
-    pmVisualAskUser(&plotLeastSquares, &isVisual);
+    pmVisualAskUser(&plotLeastSquares);
     psFree(canvas);
     psFree(canvas32);
@@ -216 +206 @@
 
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) {
-    if (!isVisual || !plotImage) return true;
+    if (!pmVisualIsVisual() || !plotImage) return true;
     if (!pmVisualInitWindow (&kapa, "PPSub:Images")) {
-        isVisual = false;
         return false;
     }
@@ -225 +214 @@
     pmVisualScaleImage(kapa, ref, "Reference", 1, true);
     pmVisualScaleImage(kapa, sub, "Subtraction", 2, true);
-    pmVisualAskUser(&plotImage, &isVisual);
+    pmVisualAskUser(&plotImage);
     return true;
 }
@@ -267 +256 @@
 
 #else
-bool pmSubtractionSetVisual (bool mode) {return true;}
 bool pmSubtractionVisualClose() {return true;}
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {return true;}

branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.h

@@ -2 +2 @@
 #define PM_SUBTRACTION_VISUAL_H
 
-bool pmSubtractionSetVisual (bool mode);
 bool pmSubtractionVisualClose();
 bool pmSubtractionVisualPlotConvKernels(psImage *convKernels);

branches/cnb_branch_20090215/psModules/src/objects/Makefile.am

@@ -36 +36 @@
         pmSourceIO_PS1_DEV_1.c \
         pmSourceIO_PS1_CAL_0.c \
+        pmSourceIO_CMF_PS1_V1.c \
         pmSourcePlots.c \
         pmSourcePlotPSFModel.c \

branches/cnb_branch_20090215/psModules/src/objects/pmModel.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-06 20:54:14 $
+ *  @version $Revision: 1.27.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -59 +59 @@
     tmp->chisqNorm = 0.0;
     tmp->nDOF  = 0;
+    tmp->nPix  = 0;
     tmp->nIter = 0;
     tmp->radiusFit = 0;
@@ -234 +235 @@
 
     if (model->residuals) {
-        DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
-        DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
-        Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
-        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
-        Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
-        Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
-        if (Ro) {
-            NX = model->residuals->Ro->numCols;
-            NY = model->residuals->Ro->numRows;
-        }           
+        DX = xBin*(image->col0 - xCenter - dx) + model->residuals->xCenter + 0.5;
+        DY = yBin*(image->row0 - yCenter - dy) + model->residuals->yCenter + 0.5;
+        Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
+        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
+        Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
+        Rm = (model->residuals->mask) ? model->residuals->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
+        if (Ro) {
+            NX = model->residuals->Ro->numCols;
+            NY = model->residuals->Ro->numRows;
+        }
     }
 
@@ -255 +256 @@
 
             // XXX should we use using 0.5 pixel offset?
-            // Convert to coordinate in parent image, with offset (dx,dy)
+            // Convert to coordinate in parent image, with offset (dx,dy)
             imageCol = ix + image->col0 - dx;
             imageRow = iy + image->row0 - dy;
@@ -275 +276 @@
                 float rx = xBin*ix + DX;
 
-                int rx0 = rx - 0.5;
-                int rx1 = rx + 0.5;
-                int ry0 = ry - 0.5;
-                int ry1 = ry + 0.5;
-
-                if (rx0 < 0) goto skip;
-                if (ry0 < 0) goto skip;
-                if (rx1 >= NX) goto skip;
-                if (ry1 >= NY) goto skip;
-
-                // these go from 0.0 to 1.0 between the centers of the pixels 
-                float fx = rx - 0.5 - rx0;
-                float Fx = 1.0 - fx;
-                float fy = ry - 0.5 - ry0;
-                float Fy = 1.0 - fy;
-
-                // check the residual image mask (if set). give up if any of the 4 pixels are masked.
-                if (Rm) {
-                    if (Rm[ry0][rx0]) goto skip;
-                    if (Rm[ry0][rx1]) goto skip;
-                    if (Rm[ry1][rx0]) goto skip;
-                    if (Rm[ry1][rx1]) goto skip;
-                }
-
-                // a possible further optimization if we re-use these values
-                // XXX allow for masked pixels, and add pixel weights
-                float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
-                float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
-                float Vo = V0*Fy + V1*fy;
-                if (!isfinite(Vo)) goto skip;
-
-                float Vx = 0.0;
-                float Vy = 0.0;
-
-                // skip Rx,Ry if Ro is masked
-                if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
-                    V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
-                    V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
-                    Vx = V0*Fy + V1*fy;
-
-                    V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
-                    V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
-                    Vy = V0*Fy + V1*fy;
-                }
-                if (!isfinite(Vx)) goto skip;
-                if (!isfinite(Vy)) goto skip;
-
-                // 2D residual variations are set for the true source position
-                pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
+                int rx0 = rx - 0.5;
+                int rx1 = rx + 0.5;
+                int ry0 = ry - 0.5;
+                int ry1 = ry + 0.5;
+
+                if (rx0 < 0) goto skip;
+                if (ry0 < 0) goto skip;
+                if (rx1 >= NX) goto skip;
+                if (ry1 >= NY) goto skip;
+
+                // these go from 0.0 to 1.0 between the centers of the pixels
+                float fx = rx - 0.5 - rx0;
+                float Fx = 1.0 - fx;
+                float fy = ry - 0.5 - ry0;
+                float Fy = 1.0 - fy;
+
+                // check the residual image mask (if set). give up if any of the 4 pixels are masked.
+                if (Rm) {
+                    if (Rm[ry0][rx0]) goto skip;
+                    if (Rm[ry0][rx1]) goto skip;
+                    if (Rm[ry1][rx0]) goto skip;
+                    if (Rm[ry1][rx1]) goto skip;
+                }
+
+                // a possible further optimization if we re-use these values
+                // XXX allow for masked pixels, and add pixel weights
+                float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
+                float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
+                float Vo = V0*Fy + V1*fy;
+                if (!isfinite(Vo)) goto skip;
+
+                float Vx = 0.0;
+                float Vy = 0.0;
+
+                // skip Rx,Ry if Ro is masked
+                if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
+                    V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
+                    V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
+                    Vx = V0*Fy + V1*fy;
+
+                    V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
+                    V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
+                    Vy = V0*Fy + V1*fy;
+                }
+                if (!isfinite(Vx)) goto skip;
+                if (!isfinite(Vy)) goto skip;
+
+                // 2D residual variations are set for the true source position
+                pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
             }
 
-        skip:
+        skip:
             // add or subtract the value
             if (add) {

branches/cnb_branch_20090215/psModules/src/objects/pmModel.h

@@ -5 +5 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-27 06:39:38 $
+ * @version $Revision: 1.18.8.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  *
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -92 +92 @@
     float mag;                          ///< integrated model magnitude
     float magErr;                       ///< integrated model magnitude error
+    int nPix;                           ///< number of pixels used for fit
     int nDOF;                           ///< number of degrees of freedom
     int nIter;                          ///< number of iterations to reach min
@@ -150 +151 @@
     pmModel *model,                     ///< The input pmModel
     pmModelOpMode mode,                 ///< mode to control how the model is added into the image
-    psImageMaskType maskVal             ///< Value to mask
+    psImageMaskType maskVal             ///< Value to mask
 );
 
@@ -168 +169 @@
     pmModel *model,                     ///< The input pmModel
     pmModelOpMode mode,                 ///< mode to control how the model is added into the image
-    psImageMaskType maskVal             ///< Value to mask
+    psImageMaskType maskVal             ///< Value to mask
 );
 

branches/cnb_branch_20090215/psModules/src/objects/pmPeaks.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.25.14.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 AddPeak(): A private function which allocates a psArray, if the peaks
 argument is NULL, otherwise it adds the peak to that array.
-XXX EAM : row,col now refer to image coords, NOT parent (since this is private) 
+XXX EAM : row,col now refer to image coords, NOT parent (since this is private)
 XXX EAM : now also calculates fractional peak positions from 3x3 bicube region
 *****************************************************************************/
@@ -63 +63 @@
         peak->xf = min.x + ix + image->col0;
         peak->yf = min.y + iy + image->row0;
-        peak->dx = NAN;
-        peak->dy = NAN;
-        
-        // xf,yf must land on image with 0 pixel border
-        peak->xf = PS_MAX (PS_MIN (peak->xf, image->numCols - 1), image->col0);
-        peak->yf = PS_MAX (PS_MIN (peak->yf, image->numRows - 1), image->row0);
+
+        // These errors are fractional errors, and should be scaled by the
+        // error on the peak pixel (see, eg, psphotFindPeaks)
+        peak->dx = min.xErr;
+        peak->dy = min.yErr;
+
+        // xf,yf must land on image with 0 pixel border
+        peak->xf = PS_MAX (PS_MIN (peak->xf, image->numCols - 1), image->col0);
+        peak->yf = PS_MAX (PS_MIN (peak->yf, image->numRows - 1), image->row0);
     } else {
         peak->xf = ix;
         peak->yf = iy;
-        peak->dx = NAN;
-        peak->dy = NAN;
+        peak->dx = NAN;
+        peak->dy = NAN;
     }
 
@@ -86 +89 @@
 getRowVectorFromImage(): a private function which simply returns a
 psVector containing the specified row of data from the psImage.
- 
+
 XXX: Is there a better way to do this?
 XXX EAM: does this really need to alloc a new vector???
@@ -245 +248 @@
 above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
 the location (x value) of all peaks.
- 
+
 XXX: What types should be supported?  Only F32 is implemented.
- 
+
 XXX: We currently step through the input vector twice; once to determine the
 size of the output vector, then to set the values of the output vector.
@@ -253 +256 @@
 *****************************************************************************/
 psVector *pmPeaksInVector(const psVector *vector,
-                         psF32 threshold)
+                         psF32 threshold)
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
@@ -347 +350 @@
 above the given threshold.  Returns a psArray containing location (x/y value)
 of all peaks.
- 
+
 XXX: I'm not convinced the peak type definition in the SDRS is mutually
 exclusive.  Some peaks can have multiple types.  Edges for sure.  Also, a
 digonal line with the same value at each point will have a peak for every
 point on that line.
- 
+
 XXX: This does not work if image has either a single row, or a single column.
- 
+
 The peak is returned in the image parent coordinates
 

branches/cnb_branch_20090215/psModules/src/objects/pmSource.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-06 02:31:25 $
+ *  @version $Revision: 1.69.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -113 +113 @@
     source->type = PM_SOURCE_TYPE_UNKNOWN;
     source->mode = PM_SOURCE_MODE_DEFAULT;
+    source->tmpFlags = 0;
     source->extpars = NULL;
     source->region = psRegionSet(NAN, NAN, NAN, NAN);
@@ -534 +535 @@
     float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
 
+    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
+
     // XXX allow clump size to be scaled relative to sigmas?
     // make rough IDs based on clumpX,Y,DX,DY
@@ -540 +543 @@
         pmSource *source = (pmSource *) sources->data[i];
 
-        if (source->peak->x < region->x0) continue;
-        if (source->peak->x > region->x1) continue;
-        if (source->peak->y < region->y0) continue;
-        if (source->peak->y > region->y1) continue;
-
-        source->peak->type = 0;
+        // psf clumps are found for image subregions:
+        // skip sources not in this region
+        if (source->peak->x <  region->x0) continue;
+        if (source->peak->x >= region->x1) continue;
+        if (source->peak->y <  region->y0) continue;
+        if (source->peak->y >= region->y1) continue;
+
+        // should be set by pmSourceAlloc
+        psAssert (source->type == PM_SOURCE_TYPE_UNKNOWN, "source type was not init-ed?");
 
         // we are basically classifying by moments; use the default if not found
         if (!source->moments) {
             source->type = PM_SOURCE_TYPE_STAR;
+            psAssert (source->mode & noMoments, "why is this source missing moments?");
             Nstar++;
             continue;
@@ -1049 +1056 @@
 }
 
+// sort by Seq (ascending)
+int pmSourceSortBySeq (const void **a, const void **b)
+{
+    pmSource *A = *(pmSource **)a;
+    pmSource *B = *(pmSource **)b;
+
+    int iA = A->seq;
+    int iB = B->seq;
+
+    int diff = iA - iB;
+    if (diff > 0) return (+1);
+    if (diff < 0) return (-1);
+    return (0);
+}
+
 pmSourceMode pmSourceModeFromString (const char *name) {
   if (!strcasecmp (name, "DEFAULT"   )) return PM_SOURCE_MODE_DEFAULT;
@@ -1066 +1088 @@
   if (!strcasecmp (name, "CRLIMIT"   )) return PM_SOURCE_MODE_CR_LIMIT;
   if (!strcasecmp (name, "EXTLIMIT"  )) return PM_SOURCE_MODE_EXT_LIMIT;
-  if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED;
   return PM_SOURCE_MODE_DEFAULT;
 }
@@ -1088 +1109 @@
     case PM_SOURCE_MODE_CR_LIMIT   : return psStringCopy ("CRLIMIT"   );
     case PM_SOURCE_MODE_EXT_LIMIT  : return psStringCopy ("EXTLIMIT"  );
-    case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED");
     default:
       return NULL;

branches/cnb_branch_20090215/psModules/src/objects/pmSource.h

@@ -3 +3 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-02-06 02:31:25 $
+ * @version $Revision: 1.28.4.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -35 +35 @@
 } pmSourceType;
 
+// bit flags to distinguish analysis results
 typedef enum {
-    PM_SOURCE_MODE_DEFAULT    = 0x0000, ///<
-    PM_SOURCE_MODE_PSFMODEL   = 0x0001, ///< Source fitted with a psf model (linear or non-linear)
-    PM_SOURCE_MODE_EXTMODEL   = 0x0002, ///< Source fitted with an extended-source model
-    PM_SOURCE_MODE_FITTED     = 0x0004, ///< Source fitted with non-linear model (PSF or EXT; good or bad)
-    PM_SOURCE_MODE_FAIL       = 0x0008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero)
-    PM_SOURCE_MODE_POOR       = 0x0010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?)
-    PM_SOURCE_MODE_PAIR       = 0x0020, ///< Source fitted with a double psf
-    PM_SOURCE_MODE_PSFSTAR    = 0x0040, ///< Source used to define PSF model
-    PM_SOURCE_MODE_SATSTAR    = 0x0080, ///< Source model peak is above saturation
-    PM_SOURCE_MODE_BLEND      = 0x0100, ///< Source is a blend with other sourcers
-    PM_SOURCE_MODE_EXTERNAL   = 0x0200, ///< Source based on supplied input position
-    PM_SOURCE_MODE_BADPSF     = 0x0400, ///< Failed to get good estimate of object's PSF
-    PM_SOURCE_MODE_DEFECT     = 0x0800, ///< Source is thought to be a defect
-    PM_SOURCE_MODE_SATURATED  = 0x1000, ///< Source is thought to be saturated pixels (bleed trail)
-    PM_SOURCE_MODE_CR_LIMIT   = 0x2000, ///< Source has crNsigma above limit
-    PM_SOURCE_MODE_EXT_LIMIT  = 0x4000, ///< Source has extNsigma above limit
-    PM_SOURCE_MODE_SUBTRACTED = 0x8000, ///< XXX this flag is actually only used internally (move)
+    PM_SOURCE_MODE_DEFAULT          = 0x00000000, ///< Initial value: resets all bits
+    PM_SOURCE_MODE_PSFMODEL         = 0x00000001, ///< Source fitted with a psf model (linear or non-linear)
+    PM_SOURCE_MODE_EXTMODEL         = 0x00000002, ///< Source fitted with an extended-source model
+    PM_SOURCE_MODE_FITTED           = 0x00000004, ///< Source fitted with non-linear model (PSF or EXT; good or bad)
+    PM_SOURCE_MODE_FAIL             = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero)
+    PM_SOURCE_MODE_POOR             = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?)
+    PM_SOURCE_MODE_PAIR             = 0x00000020, ///< Source fitted with a double psf
+    PM_SOURCE_MODE_PSFSTAR          = 0x00000040, ///< Source used to define PSF model
+    PM_SOURCE_MODE_SATSTAR          = 0x00000080, ///< Source model peak is above saturation
+    PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sourcers
+    PM_SOURCE_MODE_EXTERNAL         = 0x00000200, ///< Source based on supplied input position
+    PM_SOURCE_MODE_BADPSF           = 0x00000400, ///< Failed to get good estimate of object's PSF
+    PM_SOURCE_MODE_DEFECT           = 0x00000800, ///< Source is thought to be a defect
+    PM_SOURCE_MODE_SATURATED        = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail)
+    PM_SOURCE_MODE_CR_LIMIT         = 0x00002000, ///< Source has crNsigma above limit
+    PM_SOURCE_MODE_EXT_LIMIT        = 0x00004000, ///< Source has extNsigma above limit
+    PM_SOURCE_MODE_MOMENTS_FAILURE  = 0x00008000, ///< could not measure the moments
+    PM_SOURCE_MODE_SKY_FAILURE      = 0x00010000, ///< could not measure the local sky
+    PM_SOURCE_MODE_SKYVAR_FAILURE   = 0x00020000, ///< could not measure the local sky variance
+    PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N
+    PM_SOURCE_MODE_BIG_RADIUS       = 0x00100000, ///< poor moments for small radius, try large radius
+    PM_SOURCE_MODE_AP_MAGS          = 0x00200000, ///< source has an aperture magnitude
+    PM_SOURCE_MODE_BLEND_FIT        = 0x00400000, ///< source was fitted as a blend
+    PM_SOURCE_MODE_EXTENDED_FIT     = 0x00800000, ///< full extended fit was used
+    PM_SOURCE_MODE_EXTENDED_STATS   = 0x01000000, ///< extended aperture stats calculated
+    PM_SOURCE_MODE_LINEAR_FIT       = 0x02000000, ///< source fitted with the linear fit
+    PM_SOURCE_MODE_NONLINEAR_FIT    = 0x04000000, ///< source fitted with the non-linear fit
+    PM_SOURCE_MODE_RADIAL_FLUX      = 0x08000000, ///< radial flux measurements calculated
+    PM_SOURCE_MODE_SIZE_SKIPPED     = 0x10000000, ///< size could not be determined
 } pmSourceMode;
+
+typedef enum {
+    PM_SOURCE_TMPF_MODEL_GUESS = 0x0001,
+    PM_SOURCE_TMPF_SUBTRACTED  = 0x0002,
+} pmSourceTmpF;
 
 /** pmSource data structure
@@ -68 +86 @@
  */
 struct pmSource {
-    const int id;                       ///< Unique ID for object
-    int seq;                            ///< ID for output (generated on write)
+    const int id;                       ///< Unique ID for object (generated on alloc)
+    int seq;                            ///< ID for output (generated on write OR set on read)
     pmPeak  *peak;                      ///< Description of peak pixel.
     psImage *pixels;                    ///< Rectangular region including object pixels.
@@ -83 +101 @@
     pmSourceType type;                  ///< Best identification of object.
     pmSourceMode mode;                  ///< analysis flags set for object.
+    pmSourceTmpF tmpFlags;              ///< internal-only flags
     psArray *blends;                    ///< collection of sources thought to be confused with object
     float psfMag;                       ///< calculated from flux in modelPSF
@@ -240 +259 @@
 bool pmSourceCachePSF (pmSource *source, psImageMaskType maskVal);
 
-int             pmSourceSortBySN (const void **a, const void **b);
-int             pmSourceSortByY (const void **a, const void **b);
+int  pmSourceSortBySN (const void **a, const void **b);
+int  pmSourceSortByY (const void **a, const void **b);
+int  pmSourceSortBySeq (const void **a, const void **b);
 
 pmSourceMode pmSourceModeFromString (const char *name);

branches/cnb_branch_20090215/psModules/src/objects/pmSourceFitModel.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-06 02:31:25 $
+ *  @version $Revision: 1.30.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -184 +184 @@
     model->chisq = myMin->value;
     model->nIter = myMin->iter;
+    model->nPix  = y->n;
     model->nDOF  = y->n - nParams;
     model->flags |= PM_MODEL_STATUS_FITTED;

branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.69.14.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -488 +488 @@
             }
             if (!strcmp (exttype, "PS1_DEV_1")) {
-                status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname, xsrcname);
+                status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_CAL_0")) {
-                status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname, xsrcname);
+                status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
+            }
+            if (!strcmp (exttype, "PS1_V1")) {
+                status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
             }
             if (xsrcname) {
@@ -500 +503 @@
                   status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
               }
+              if (!strcmp (exttype, "PS1_V1")) {
+                  status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
+              }
             }
             if (xfitname) {
@@ -507 +513 @@
               if (!strcmp (exttype, "PS1_CAL_0")) {
                   status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
+              }
+              if (!strcmp (exttype, "PS1_V1")) {
+                  status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
               }
             }
@@ -929 +938 @@
                 sources = pmSourcesRead_PS1_DEV_1 (file->fits, hdu->header);
             }
-            if (!strcmp (exttype, "PS1_CAL_)")) {
-                sources = pmSourcesRead_PS1_CAL_0 (file->fits, hdu->header);
+            if (!strcmp (exttype, "PS1_V1")) {
+                sources = pmSourcesRead_CMF_PS1_V1 (file->fits, hdu->header);
             }
         }

branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.h

@@ -4 +4 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-11-14 02:11:45 $
+ * @version $Revision: 1.19.16.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  *
@@ -26 +26 @@
 bool pmSourcesWrite_SMPDATA (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
 bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
-bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname);
+
+bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname);
 
-bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname);
+bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
 bool pmSourcesWrite_PS1_CAL_0_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe);
 bool pmSourcesWrite_PS1_CAL_0_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname);
+
+bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname);
 
 bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
@@ -42 +47 @@
 psArray *pmSourcesRead_PS1_DEV_1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_PS1_CAL_0 (psFits *fits, psMetadata *header);
+psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
 
 bool pmSourcesWritePSFs (psArray *sources, char *filename);
@@ -67 +73 @@
 bool pmReadoutCheckDataStatusForSources (const pmReadout *readout);
 
+bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos);
+
 /// @}
 # endif /* PM_SOURCE_IO_H */

branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.2.14.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49 +49 @@
 // XXX how do I generate the source tables which I need to send to PSPS?
 
-bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader,
-                               psMetadata *tableHeader, char *extname, char *xsrcname)
+bool pmSourcesWrite_PS1_CAL_0 (psFits *fits, pmReadout *readout, psArray *sources,
+                               psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
@@ -68 +68 @@
     pmFPA  *fpa  = chip->parent;
     if (!chip->toFPA || !fpa->toTPA || !fpa->toSky) {
-        psWarning ("astrometry calibration is missing, no calibrated coords");
+        psWarning ("astrometry calibration is missing, no calibrated coords");
     }
 
@@ -77 +77 @@
     float zeropt  = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
     if (!status1 || !status2 || (exptime == 0.0)) {
-        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
+        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
     } else {
-        magOffset = zeropt + 2.5*log10(exptime);
+        magOffset = zeropt + 2.5*log10(exptime);
     }
 
@@ -90 +90 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        if (source->seq == -1) {
-            source->seq = i;
-        }
+        if (source->seq == -1) {
+            source->seq = i;
+        }
 
         // no difference between PSF and non-PSF model
@@ -104 +104 @@
             xErr = dPAR[PM_PAR_XPOS];
             yErr = dPAR[PM_PAR_YPOS];
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-
-            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
+            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
+                axes = pmPSF_ModelToAxes (PAR, 20.0);
+            } else {
+                axes.major = NAN;
+                axes.minor = NAN;
+                axes.theta = NAN;
+            }
+            chisq = model->chisq;
+
+            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
+            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
         } else {
             xPos = source->peak->xf;
@@ -123 +123 @@
             axes.minor = NAN;
             axes.theta = NAN;
-            chisq = NAN;
-            errMag = NAN;
-        }
-
-        float calMag = source->psfMag + magOffset;
+            chisq = NAN;
+            errMag = NAN;
+        }
+
+        float calMag = source->psfMag + magOffset;
         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
         psS16 nImageOverlap = 1;
 
-        // generate RA,DEC
-        psPlane ptCH, ptFP, ptTP;
-        psSphere ptSky;
-
-        ptCH.x = xPos;
-        ptCH.y = yPos;
-        if (chip->toFPA && fpa->toTPA && fpa->toSky) {
-            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
-            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
-            psDeproject (&ptSky, &ptTP, fpa->toSky);
-        } else {
-            ptSky.r = NAN;
-            ptSky.d = NAN;
-        }
+        // generate RA,DEC
+        psPlane ptCH, ptFP, ptTP;
+        psSphere ptSky;
+
+        ptCH.x = xPos;
+        ptCH.y = yPos;
+        if (chip->toFPA && fpa->toTPA && fpa->toSky) {
+            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+            psDeproject (&ptSky, &ptTP, fpa->toSky);
+        } else {
+            ptSky.r = NAN;
+            ptSky.d = NAN;
+        }
 
         row = psMetadataAlloc ();
@@ -150 +150 @@
         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
-        // XXX need to do the error propagation correctly..
+        // XXX need to do the error propagation correctly..
         // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF_SIG",       PS_DATA_F32, "Sigma of PSF fit RA",                     dRA);
         // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF_SIG",      PS_DATA_F32, "Sigma of PSF fit DEC",                    dDEC);
@@ -170 +170 @@
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32,  "Nsigma deviations from PSF to EXT",         source->extNsigma);
 
-        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
-        // subtracted
-
-        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
-        // edge; 3) any pixels in the 3x3 peak region are masked; 
+        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
+        // subtracted
+
+        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
+        // edge; 3) any pixels in the 3x3 peak region are masked;
 
         // XXX these should be major and minor, not 'x' and 'y'
@@ -277 +277 @@
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
         source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
-        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
 
         pmPSF_AxesToModel (PAR, axes);
@@ -289 +289 @@
         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
 
-        // note that some older versions used PSF_PROBABILITY: this was not well defined.
+        // note that some older versions used PSF_PROBABILITY: this was not well defined.
         model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
         source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
@@ -321 +321 @@
     pmFPA  *fpa  = chip->parent;
     if (!chip->toFPA || !fpa->toTPA || !fpa->toSky) {
-        psWarning ("astrometry calibration is missing, no calibrated coords");
+        psWarning ("astrometry calibration is missing, no calibrated coords");
     }
 
@@ -331 +331 @@
     float zeropt  = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
     if (!status1 || !status2 || (exptime == 0.0)) {
-        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
+        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
     } else {
-        magOffset = zeropt + 2.5*log10(exptime);
-        calMags = true;
+        magOffset = zeropt + 2.5*log10(exptime);
+        calMags = true;
     }
 
@@ -360 +360 @@
     // we write out all sources, regardless of quality.  the source flags tell us the state
     for (int i = 0; i < sources->n; i++) {
-        // skip source if it is not a ext sourc
-        // XXX we have two places that extended source parameters are measured:
-        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
-        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
-        // should we require both?
-
-        pmSource *source = sources->data[i];
-
-        // skip sources without measurements
-        if (source->extpars == NULL) continue;
-
-        // we require a PSF model fit (ignore the real crud)
-        pmModel *model = source->modelPSF;
-        if (model == NULL) continue;
-
-        // XXX I need to split the extended models from the extended aperture measurements
-        PAR = model->params->data.F32;
-        dPAR = model->dparams->data.F32;
-        xPos = PAR[PM_PAR_XPOS];
-        yPos = PAR[PM_PAR_YPOS];
-        xErr = dPAR[PM_PAR_XPOS];
-        yErr = dPAR[PM_PAR_YPOS];
-
-        // generate RA,DEC
-        psPlane ptCH, ptFP, ptTP;
-        psSphere ptSky;
-
-        ptCH.x = xPos;
-        ptCH.y = yPos;
-        if (chip->toFPA && fpa->toTPA && fpa->toSky) {
-            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
-            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
-            psDeproject (&ptSky, &ptTP, fpa->toSky);
-        } else {
-            ptSky.r = NAN;
-            ptSky.d = NAN;
-        }
+        // skip source if it is not a ext sourc
+        // XXX we have two places that extended source parameters are measured:
+        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
+        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
+        // should we require both?
+
+        pmSource *source = sources->data[i];
+
+        // skip sources without measurements
+        if (source->extpars == NULL) continue;
+
+        // we require a PSF model fit (ignore the real crud)
+        pmModel *model = source->modelPSF;
+        if (model == NULL) continue;
+
+        // XXX I need to split the extended models from the extended aperture measurements
+        PAR = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
+        xPos = PAR[PM_PAR_XPOS];
+        yPos = PAR[PM_PAR_YPOS];
+        xErr = dPAR[PM_PAR_XPOS];
+        yErr = dPAR[PM_PAR_YPOS];
+
+        // generate RA,DEC
+        psPlane ptCH, ptFP, ptTP;
+        psSphere ptSky;
+
+        ptCH.x = xPos;
+        ptCH.y = yPos;
+        if (chip->toFPA && fpa->toTPA && fpa->toSky) {
+            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+            psDeproject (&ptSky, &ptTP, fpa->toSky);
+        } else {
+            ptSky.r = NAN;
+            ptSky.d = NAN;
+        }
 
         row = psMetadataAlloc ();
@@ -410 +410 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
 
-        // Petrosian measurements
-        // XXX insert header data: petrosian ref radius, flux ratio
-        if (doPetrosian) {
-            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
-            if (petrosian) {
-                if (calMags) {
-                    psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_CAL",    PS_DATA_F32, "Petrosian Magnitude (calibrated)", petrosian->mag + magOffset);
-                } else {
-                    psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_INST",   PS_DATA_F32, "Petrosian Magnitude (instrumental)", petrosian->mag);
-                }
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
-            }
-        } 
-
-        // Kron measurements
-        if (doKron) {
-            pmSourceKronValues *kron = source->extpars->kron;
-            if (kron) {
-                if (calMags) {
-                    psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_CAL",  PS_DATA_F32, "Kron Magnitude",     kron->mag + magOffset);
-                } else {
-                    psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_INST", PS_DATA_F32, "Kron Magnitude",     kron->mag);
-                }
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
-            }
-        }
-
-        // Isophot measurements
-        // XXX insert header data: isophotal level
-        if (doIsophotal) {
-            pmSourceIsophotalValues *isophot = source->extpars->isophot;
-            if (isophot) {
-                if (calMags) {
-                    psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_CAL",    PS_DATA_F32, "Isophot Magnitude (calibrated)",   isophot->mag + magOffset);
-                } else {
-                    psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_INST",   PS_DATA_F32, "Isophot Magnitude (uncalibrated)", isophot->mag);
-                }
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
-            }
-        }
-
-        // Flux Annuli
-        if (doAnnuli) {
-            pmSourceAnnuli *annuli = source->extpars->annuli;
-            if (annuli) {
-                psVector *fluxVal = annuli->flux;
-                psVector *fluxErr = annuli->fluxErr;
-                psVector *fluxVar = annuli->fluxVar;
-
-                for (int j = 0; j < fluxVal->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
-                } 
-            } else {
-                for (int j = 0; j < radialBinsLower->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
-                } 
-            }
-        }
-
-        psArrayAdd (table, 100, row);
-        psFree (row);
+        // Petrosian measurements
+        // XXX insert header data: petrosian ref radius, flux ratio
+        if (doPetrosian) {
+            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
+            if (petrosian) {
+                if (calMags) {
+                    psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_CAL",    PS_DATA_F32, "Petrosian Magnitude (calibrated)", petrosian->mag + magOffset);
+                } else {
+                    psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_INST",   PS_DATA_F32, "Petrosian Magnitude (instrumental)", petrosian->mag);
+                }
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
+            }
+        }
+
+        // Kron measurements
+        if (doKron) {
+            pmSourceKronValues *kron = source->extpars->kron;
+            if (kron) {
+                if (calMags) {
+                    psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_CAL",  PS_DATA_F32, "Kron Magnitude",     kron->mag + magOffset);
+                } else {
+                    psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_INST", PS_DATA_F32, "Kron Magnitude",     kron->mag);
+                }
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
+            }
+        }
+
+        // Isophot measurements
+        // XXX insert header data: isophotal level
+        if (doIsophotal) {
+            pmSourceIsophotalValues *isophot = source->extpars->isophot;
+            if (isophot) {
+                if (calMags) {
+                    psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_CAL",    PS_DATA_F32, "Isophot Magnitude (calibrated)",   isophot->mag + magOffset);
+                } else {
+                    psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_INST",   PS_DATA_F32, "Isophot Magnitude (uncalibrated)", isophot->mag);
+                }
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
+            }
+        }
+
+        // Flux Annuli
+        if (doAnnuli) {
+            pmSourceAnnuli *annuli = source->extpars->annuli;
+            if (annuli) {
+                psVector *fluxVal = annuli->flux;
+                psVector *fluxErr = annuli->fluxErr;
+                psVector *fluxVar = annuli->fluxVar;
+
+                for (int j = 0; j < fluxVal->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
+                }
+            } else {
+                for (int j = 0; j < radialBinsLower->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
+                }
+            }
+        }
+
+        psArrayAdd (table, 100, row);
+        psFree (row);
     }
 
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
-        psFree (outhead);
-        psFree (table);
-        return true;
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
     }
 
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
     }
     psFree (outhead);
@@ -540 +540 @@
     pmFPA  *fpa  = chip->parent;
     if (!chip->toFPA || !fpa->toTPA || !fpa->toSky) {
-        psWarning ("astrometry calibration is missing, no calibrated coords");
+        psWarning ("astrometry calibration is missing, no calibrated coords");
     }
 
@@ -549 +549 @@
     float zeropt  = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
     if (!status1 || !status2 || (exptime == 0.0)) {
-        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
+        psWarning ("exposure time or measured zero point not found for a readout, no calibrated mags");
     } else {
-        magOffset = zeropt + 2.5*log10(exptime);
+        magOffset = zeropt + 2.5*log10(exptime);
     }
 
@@ -566 +566 @@
     int nParamMax = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->modelFits == NULL) continue;
-        for (int j = 0; j < source->modelFits->n; j++) {
-            pmModel *model = source->modelFits->data[j];
-            assert (model);
-            nParamMax = PS_MAX (nParamMax, model->params->n);
-        }
+        pmSource *source = sources->data[i];
+        if (source->modelFits == NULL) continue;
+        for (int j = 0; j < source->modelFits->n; j++) {
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+            nParamMax = PS_MAX (nParamMax, model->params->n);
+        }
     }
 
@@ -580 +580 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        // XXX if no model fits are saved, write out modelEXT?
-        if (source->modelFits == NULL) continue;
-
-        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
-        for (int j = 0; j < source->modelFits->n; j++) {
-
-            // choose the convolved EXT model, if available, otherwise the simple one
-            pmModel *model = source->modelFits->data[j];
-            assert (model);
-
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            xErr = dPAR[PM_PAR_XPOS];
-            yErr = dPAR[PM_PAR_YPOS];
-
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
-
-            // generate RA,DEC
-            psPlane ptCH, ptFP, ptTP;
-            psSphere ptSky;
-
-            ptCH.x = xPos;
-            ptCH.y = yPos;
-            if (chip->toFPA && fpa->toTPA && fpa->toSky) {
-                psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
-                psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
-                psDeproject (&ptSky, &ptTP, fpa->toSky);
-            } else {
-                ptSky.r = NAN;
-                ptSky.d = NAN;
-            }
-
-            row = psMetadataAlloc ();
-
-            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
-            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model RA (degrees)",                     ptSky.r);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model DEC (degrees)",                    ptSky.d);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_CAL_MAG",      0, "EXT fit calibrated magnitude",               model->mag + magOffset);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
-
-            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
-            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
-
-            // XXX these should be major and minor, not 'x' and 'y'
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
-
-            // write out the other generic parameters
-            for (int k = 0; k < nParamMax; k++) {
-                if (k == PM_PAR_I0) continue;
-                if (k == PM_PAR_SKY) continue;
-                if (k == PM_PAR_XPOS) continue;
-                if (k == PM_PAR_YPOS) continue;
-                if (k == PM_PAR_SXX) continue;
-                if (k == PM_PAR_SXY) continue;
-                if (k == PM_PAR_SYY) continue;
-
-                snprintf (name, 64, "EXT_PAR_%02d", k);
-
-                if (k < model->params->n) {
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
-                } else {
-                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
-                }
-            }
-
-            // XXX other parameters which may be set.
-            // XXX flag / value to define the model
-            // XXX write out the model type, fit status flags
-
-            psArrayAdd (table, 100, row);
-            psFree (row);
-        }
+        pmSource *source = sources->data[i];
+
+        // XXX if no model fits are saved, write out modelEXT?
+        if (source->modelFits == NULL) continue;
+
+        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
+        for (int j = 0; j < source->modelFits->n; j++) {
+
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+
+            PAR = model->params->data.F32;
+            dPAR = model->dparams->data.F32;
+            xPos = PAR[PM_PAR_XPOS];
+            yPos = PAR[PM_PAR_YPOS];
+            xErr = dPAR[PM_PAR_XPOS];
+            yErr = dPAR[PM_PAR_YPOS];
+
+            axes = pmPSF_ModelToAxes (PAR, 20.0);
+
+            // generate RA,DEC
+            psPlane ptCH, ptFP, ptTP;
+            psSphere ptSky;
+
+            ptCH.x = xPos;
+            ptCH.y = yPos;
+            if (chip->toFPA && fpa->toTPA && fpa->toSky) {
+                psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+                psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+                psDeproject (&ptSky, &ptTP, fpa->toSky);
+            } else {
+                ptSky.r = NAN;
+                ptSky.d = NAN;
+            }
+
+            row = psMetadataAlloc ();
+
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model RA (degrees)",                     ptSky.r);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model DEC (degrees)",                    ptSky.d);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_CAL_MAG",      0, "EXT fit calibrated magnitude",               model->mag + magOffset);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
+
+            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
+            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
+
+            // XXX these should be major and minor, not 'x' and 'y'
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
+
+            // write out the other generic parameters
+            for (int k = 0; k < nParamMax; k++) {
+                if (k == PM_PAR_I0) continue;
+                if (k == PM_PAR_SKY) continue;
+                if (k == PM_PAR_XPOS) continue;
+                if (k == PM_PAR_YPOS) continue;
+                if (k == PM_PAR_SXX) continue;
+                if (k == PM_PAR_SXY) continue;
+                if (k == PM_PAR_SYY) continue;
+
+                snprintf (name, 64, "EXT_PAR_%02d", k);
+
+                if (k < model->params->n) {
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
+                } else {
+                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
+                }
+            }
+
+            // XXX other parameters which may be set.
+            // XXX flag / value to define the model
+            // XXX write out the model type, fit status flags
+
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
     }
 
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
-        psFree (outhead);
-        psFree (table);
-        return true;
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
     }
 
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
     }
     psFree (outhead);

branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.14.14.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49 +49 @@
 // XXX how do I generate the source tables which I need to send to PSPS?
 
-bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader,
-                               psMetadata *tableHeader, char *extname, char *xsrcname)
+bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources,
+                               psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
@@ -73 +73 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        if (source->seq == -1) {
-            source->seq = i;
-        }
+        if (source->seq == -1) {
+            source->seq = i;
+        }
 
         // no difference between PSF and non-PSF model
@@ -87 +87 @@
             xErr = dPAR[PM_PAR_XPOS];
             yErr = dPAR[PM_PAR_YPOS];
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-
-            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
+            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
+                axes = pmPSF_ModelToAxes (PAR, 20.0);
+            } else {
+                axes.major = NAN;
+                axes.minor = NAN;
+                axes.theta = NAN;
+            }
+            chisq = model->chisq;
+
+            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
+            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
         } else {
             xPos = source->peak->xf;
@@ -106 +106 @@
             axes.minor = NAN;
             axes.theta = NAN;
-            chisq = NAN;
-            errMag = NAN;
+            chisq = NAN;
+            errMag = NAN;
         }
 
@@ -130 +130 @@
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32,  "Nsigma deviations from PSF to EXT",         source->extNsigma);
 
-        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
-        // subtracted
-
-        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
-        // edge; 3) any pixels in the 3x3 peak region are masked; 
-
-        // CR_NSIGMA and 
+        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
+        // subtracted
+
+        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
+        // edge; 3) any pixels in the 3x3 peak region are masked;
+
+        // CR_NSIGMA and
 
         // XXX these should be major and minor, not 'x' and 'y'
@@ -247 +247 @@
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
         source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
-        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
 
         pmPSF_AxesToModel (PAR, axes);
@@ -259 +259 @@
         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
 
-        // note that some older versions used PSF_PROBABILITY: this was not well defined.
+        // note that some older versions used PSF_PROBABILITY: this was not well defined.
         model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
         source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
@@ -311 +311 @@
     // we write out all sources, regardless of quality.  the source flags tell us the state
     for (int i = 0; i < sources->n; i++) {
-        // skip source if it is not a ext sourc
-        // XXX we have two places that extended source parameters are measured:
-        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
-        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
-        // should we require both?
-
-        pmSource *source = sources->data[i];
-
-        // skip sources without measurements
-        if (source->extpars == NULL) continue;
-
-        // we require a PSF model fit (ignore the real crud)
-        pmModel *model = source->modelPSF;
-        if (model == NULL) continue;
-
-        // XXX I need to split the extended models from the extended aperture measurements
-        PAR = model->params->data.F32;
-        dPAR = model->dparams->data.F32;
-        xPos = PAR[PM_PAR_XPOS];
-        yPos = PAR[PM_PAR_YPOS];
-        xErr = dPAR[PM_PAR_XPOS];
-        yErr = dPAR[PM_PAR_YPOS];
+        // skip source if it is not a ext sourc
+        // XXX we have two places that extended source parameters are measured:
+        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
+        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
+        // should we require both?
+
+        pmSource *source = sources->data[i];
+
+        // skip sources without measurements
+        if (source->extpars == NULL) continue;
+
+        // we require a PSF model fit (ignore the real crud)
+        pmModel *model = source->modelPSF;
+        if (model == NULL) continue;
+
+        // XXX I need to split the extended models from the extended aperture measurements
+        PAR = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
+        xPos = PAR[PM_PAR_XPOS];
+        yPos = PAR[PM_PAR_YPOS];
+        xErr = dPAR[PM_PAR_XPOS];
+        yErr = dPAR[PM_PAR_YPOS];
 
         row = psMetadataAlloc ();
@@ -343 +343 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
 
-        // Petrosian measurements
-        // XXX insert header data: petrosian ref radius, flux ratio
-        if (doPetrosian) {
-            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
-            if (petrosian) {
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
-            }
-        } 
-
-        // Kron measurements
-        if (doKron) {
-            pmSourceKronValues *kron = source->extpars->kron;
-            if (kron) {
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
-            }
-        }
-
-        // Isophot measurements
-        // XXX insert header data: isophotal level
-        if (doIsophotal) {
-            pmSourceIsophotalValues *isophot = source->extpars->isophot;
-            if (isophot) {
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
-            } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
-            }
-        }
-
-        // Flux Annuli
-        if (doAnnuli) {
-            pmSourceAnnuli *annuli = source->extpars->annuli;
-            if (annuli) {
-                psVector *fluxVal = annuli->flux;
-                psVector *fluxErr = annuli->fluxErr;
-                psVector *fluxVar = annuli->fluxVar;
-
-                for (int j = 0; j < fluxVal->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
-                } 
-            } else {
-                for (int j = 0; j < radialBinsLower->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
-                } 
-            }
-        }
-
-        psArrayAdd (table, 100, row);
-        psFree (row);
+        // Petrosian measurements
+        // XXX insert header data: petrosian ref radius, flux ratio
+        if (doPetrosian) {
+            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
+            if (petrosian) {
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
+            }
+        }
+
+        // Kron measurements
+        if (doKron) {
+            pmSourceKronValues *kron = source->extpars->kron;
+            if (kron) {
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
+            }
+        }
+
+        // Isophot measurements
+        // XXX insert header data: isophotal level
+        if (doIsophotal) {
+            pmSourceIsophotalValues *isophot = source->extpars->isophot;
+            if (isophot) {
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
+            }
+        }
+
+        // Flux Annuli
+        if (doAnnuli) {
+            pmSourceAnnuli *annuli = source->extpars->annuli;
+            if (annuli) {
+                psVector *fluxVal = annuli->flux;
+                psVector *fluxErr = annuli->fluxErr;
+                psVector *fluxVar = annuli->fluxVar;
+
+                for (int j = 0; j < fluxVal->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
+                }
+            } else {
+                for (int j = 0; j < radialBinsLower->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
+                }
+            }
+        }
+
+        psArrayAdd (table, 100, row);
+        psFree (row);
     }
 
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
-        psFree (outhead);
-        psFree (table);
-        return true;
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
     }
 
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
     }
     psFree (outhead);
@@ -470 +470 @@
     int nParamMax = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->modelFits == NULL) continue;
-        for (int j = 0; j < source->modelFits->n; j++) {
-            pmModel *model = source->modelFits->data[j];
-            assert (model);
-            nParamMax = PS_MAX (nParamMax, model->params->n);
-        }
+        pmSource *source = sources->data[i];
+        if (source->modelFits == NULL) continue;
+        for (int j = 0; j < source->modelFits->n; j++) {
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+            nParamMax = PS_MAX (nParamMax, model->params->n);
+        }
     }
 
@@ -484 +484 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
-
-        // XXX if no model fits are saved, write out modelEXT?
-        if (source->modelFits == NULL) continue;
-
-        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
-        for (int j = 0; j < source->modelFits->n; j++) {
-
-            // choose the convolved EXT model, if available, otherwise the simple one
-            pmModel *model = source->modelFits->data[j];
-            assert (model);
-
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            xErr = dPAR[PM_PAR_XPOS];
-            yErr = dPAR[PM_PAR_YPOS];
-
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
-
-            row = psMetadataAlloc ();
-
-            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
-            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
-
-            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
-            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
-
-            // XXX these should be major and minor, not 'x' and 'y'
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
-
-            // write out the other generic parameters
-            for (int k = 0; k < nParamMax; k++) {
-                if (k == PM_PAR_I0) continue;
-                if (k == PM_PAR_SKY) continue;
-                if (k == PM_PAR_XPOS) continue;
-                if (k == PM_PAR_YPOS) continue;
-                if (k == PM_PAR_SXX) continue;
-                if (k == PM_PAR_SXY) continue;
-                if (k == PM_PAR_SYY) continue;
-
-                snprintf (name, 64, "EXT_PAR_%02d", k);
-
-                if (k < model->params->n) {
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
-                } else {
-                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
-                }
-            }
-
-            // XXX other parameters which may be set.
-            // XXX flag / value to define the model
-            // XXX write out the model type, fit status flags
-
-            psArrayAdd (table, 100, row);
-            psFree (row);
-        }
+        pmSource *source = sources->data[i];
+
+        // XXX if no model fits are saved, write out modelEXT?
+        if (source->modelFits == NULL) continue;
+
+        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
+        for (int j = 0; j < source->modelFits->n; j++) {
+
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+
+            PAR = model->params->data.F32;
+            dPAR = model->dparams->data.F32;
+            xPos = PAR[PM_PAR_XPOS];
+            yPos = PAR[PM_PAR_YPOS];
+            xErr = dPAR[PM_PAR_XPOS];
+            yErr = dPAR[PM_PAR_YPOS];
+
+            axes = pmPSF_ModelToAxes (PAR, 20.0);
+
+            row = psMetadataAlloc ();
+
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
+
+            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
+            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
+
+            // XXX these should be major and minor, not 'x' and 'y'
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
+
+            // write out the other generic parameters
+            for (int k = 0; k < nParamMax; k++) {
+                if (k == PM_PAR_I0) continue;
+                if (k == PM_PAR_SKY) continue;
+                if (k == PM_PAR_XPOS) continue;
+                if (k == PM_PAR_YPOS) continue;
+                if (k == PM_PAR_SXX) continue;
+                if (k == PM_PAR_SXY) continue;
+                if (k == PM_PAR_SYY) continue;
+
+                snprintf (name, 64, "EXT_PAR_%02d", k);
+
+                if (k < model->params->n) {
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
+                } else {
+                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
+                }
+            }
+
+            // XXX other parameters which may be set.
+            // XXX flag / value to define the model
+            // XXX write out the model type, fit status flags
+
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
     }
 
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
-        psFree (outhead);
-        psFree (table);
-        return true;
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
     }
 
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
     }
     psFree (outhead);

branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-06 02:31:25 $
+ *  @version $Revision: 1.49.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-19 17:59:50 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -145 +145 @@
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal);
+        pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
     }
 
@@ -164 +164 @@
     // XXX full model or just analytical?
     // XXX use pmSourceAdd instead?
-    if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
+    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
         pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
     }
@@ -170 +170 @@
     // if we are measuring aperture photometry and applying the growth correction,
     // we need to shift the flux in the selected pixels (but not the mask)
-    // psImageShift ();
     psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer
     if (mode & PM_SOURCE_PHOT_INTERP) {
@@ -186 +185 @@
         }
 
-        // XXX this is test code to verify the shift is doing the right thing (seems to be)
-        # if (0)
-            // measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-            psImage *image = source->pixels;
-            float xo = 0.0;
-            float yo = 0.0;
-            float xo2 = 0.0;
-            float yo2 = 0.0;
-            float no = 0.0;
-            for (int j = 0; j < image->numRows; j++)
-            {
-                for (int i = 0; i < image->numCols; i++) {
-                    xo += i*image->data.F32[j][i];
-                    yo += j*image->data.F32[j][i];
-                    xo2 += i*i*image->data.F32[j][i];
-                    yo2 += j*j*image->data.F32[j][i];
-                    no += image->data.F32[j][i];
-                }
-            }
-            xo /= no;
-            yo /= no;
-            xo2 = sqrt (xo2/no - xo*xo);
-            yo2 = sqrt (yo2/no - yo*yo);
-            fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-
-        // measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-            psImage *image = flux;
-            float xo = 0.0;
-            float yo = 0.0;
-            float xo2 = 0.0;
-            float yo2 = 0.0;
-            float no = 0.0;
-            for (int j = 0; j < image->numRows; j++)
-            {
-                for (int i = 0; i < image->numCols; i++) {
-                    xo += i*image->data.F32[j][i];
-                    yo += j*image->data.F32[j][i];
-                    xo2 += i*i*image->data.F32[j][i];
-                    yo2 += j*j*image->data.F32[j][i];
-                    no += image->data.F32[j][i];
-                }
-            }
-            xo /= no;
-            yo /= no;
-            xo2 = sqrt (xo2/no - xo*xo);
-            yo2 = sqrt (yo2/no - yo*yo);
-            fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-        # endif
-
     } else {
         flux = source->pixels;
@@ -269 +215 @@
     // if source was originally subtracted, re-subtract object, leave local sky
     // XXX replace with pmSourceSub...
-    if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
+    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
         pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
     }
@@ -344 +290 @@
 
 // return source aperture magnitude
-bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal)
+bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(pixWeight, false);
-    PS_ASSERT_PTR_NON_NULL(image, false);
     PS_ASSERT_PTR_NON_NULL(mask, false);
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -377 +322 @@
     Yo = params->data.F32[PM_PAR_YPOS];
 
-    dX = Xo - image->col0;
-    dP = image->numCols - dX;
+    dX = Xo - mask->col0;
+    dP = mask->numCols - dX;
     DX = PS_MAX(dX, dP);
-    NX = image->numCols;
-
-    dY = Yo - image->row0;
-    dP = image->numRows - dY;
+    NX = mask->numCols;
+
+    dY = Yo - mask->row0;
+    dP = mask->numRows - dY;
     DY = PS_MAX(dY, dP);
-    NY = image->numRows;
+    NY = mask->numRows;
 
     // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
@@ -612 +557 @@
 # endif
 
+// determine chisq, etc for linear normalization-only fit
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance,
                     psImageMaskType maskVal)
@@ -632 +578 @@
         }
     }
+    model->nPix = Npix;
     model->nDOF = Npix - 1;
     model->chisq = dC;
@@ -817 +764 @@
 }
 
+// XXX this is test code to verify the shift is doing the right thing (seems to be)
+# if (0)
+// measure centroid of unshifted gaussian (should be 16.0,16.0)
+        {
+          psImage *image = source->pixels;
+          float xo = 0.0;
+          float yo = 0.0;
+          float xo2 = 0.0;
+          float yo2 = 0.0;
+          float no = 0.0;
+          for (int j = 0; j < image->numRows; j++)
+          {
+            for (int i = 0; i < image->numCols; i++) {
+              xo += i*image->data.F32[j][i];
+              yo += j*image->data.F32[j][i];
+              xo2 += i*i*image->data.F32[j][i];
+              yo2 += j*j*image->data.F32[j][i];
+              no += image->data.F32[j][i];
+            }
+          }
+          xo /= no;
+          yo /= no;
+          xo2 = sqrt (xo2/no - xo*xo);
+          yo2 = sqrt (yo2/no - yo*yo);
+          fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
+        }
+
+// measure centroid of unshifted gaussian (should be 16.0,16.0)
+        {
+          psImage *image = flux;
+          float xo = 0.0;
+          float yo = 0.0;
+          float xo2 = 0.0;
+          float yo2 = 0.0;
+          float no = 0.0;
+          for (int j = 0; j < image->numRows; j++)
+          {
+            for (int i = 0; i < image->numCols; i++) {
+              xo += i*image->data.F32[j][i];
+              yo += j*image->data.F32[j][i];
+              xo2 += i*i*image->data.F32[j][i];
+              yo2 += j*j*image->data.F32[j][i];
+              no += image->data.F32[j][i];
+            }
+          }
+          xo /= no;
+          yo /= no;
+          xo2 = sqrt (xo2/no - xo*xo);
+          yo2 = sqrt (yo2/no - yo*yo);
+          fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
+        }
+# endif

branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.h

@@ -4 +4 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- * @date $Date: 2009-01-27 06:39:38 $
+ * @version $Revision: 1.12.8.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -46 +46 @@
     psImage *image,                     ///< image pixels to be used
     psImage *mask,                      ///< mask of pixels to ignore
-    psImageMaskType maskVal             ///< Value to mask
+    psImageMaskType maskVal             ///< Value to mask
 );
 
 bool pmSourceMagnitudesInit (psMetadata *config);
 bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal);
-bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask, psImageMaskType maskVal);
+bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *mask, psImageMaskType maskVal);
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal);
 

branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.c

@@ -9 +9 @@
 # if (HAVE_KAPA)
 # include <kapa.h>
+#include "pmVisual.h"
 
 // functions used to visualize the analysis as it goes
 // these are invoked by the -visual options
 
-static bool isVisual = false;
 static int kapa1 = -1;
 // static int kapa2 = -1;
@@ -20 +20 @@
 bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi);
 
-bool pmSourceSetVisual (bool mode) {
-
-    isVisual = mode;
-    return true;
-}
 
 bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {
 
-    KapaSection section;  // put the positive profile in one and the residuals in another? 
+    KapaSection section;  // put the positive profile in one and the residuals in another?
 
     Graphdata graphdata;
 
-    if (!isVisual) return true;
+    if (!pmVisualIsVisual()) return true;
 
     if (kapa1 == -1) {
         kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
-        if (kapa1 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            isVisual = false;
-            return false;
-        }
-    }  
+        if (kapa1 == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            pmVisualSetVisual(false);
+            return false;
+        }
+    }
 
     KapaClearPlots (kapa1);
@@ -55 +50 @@
 
     for (int i = 0; i < x->n; i++) {
-        model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
-        resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
-        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
-        min = PS_MIN (min, resid->data.F32[i]);
-        max = PS_MAX (max, resid->data.F32[i]);
-        Min = PS_MIN (min, param->data.F32[i]);
-        Max = PS_MAX (max, param->data.F32[i]);
+        model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
+        resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
+        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+        min = PS_MIN (min, resid->data.F32[i]);
+        max = PS_MAX (max, resid->data.F32[i]);
+        Min = PS_MIN (min, param->data.F32[i]);
+        Max = PS_MAX (max, param->data.F32[i]);
     }
 
@@ -70 +65 @@
     psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
     for (int i = 0; i < x->n; i++) {
-        xn->data.F32[i] = x->data.F32[i] / 5000.0;
-        yn->data.F32[i] = y->data.F32[i] / 5000.0;
-        zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
-        Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
-        Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
+        xn->data.F32[i] = x->data.F32[i] / 5000.0;
+        yn->data.F32[i] = y->data.F32[i] / 5000.0;
+        zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
+        Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
+        Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
     }
 
@@ -155 +150 @@
     fprintf (stdout, "[c]ontinue? ");
     if (!fgets(key, 8, stdin)) {
-        psWarning("Unable to read option");
+        psWarning("Unable to read option");
     }
     return true;
@@ -173 +168 @@
 
     for (int i = 0; i < xn->n; i++) {
-        xv->data.F32[2*i+0] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi) + zn->data.F32[i]*sin(theta)*sin(phi);
-        yv->data.F32[2*i+0] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi) + zn->data.F32[i]*cos(theta)*sin(phi);
-        zv->data.F32[2*i+0] = -yn->data.F32[i]*sin(phi)   + zn->data.F32[i]*cos(phi);
-        xv->data.F32[2*i+1] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi);
-        yv->data.F32[2*i+1] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi);
-        zv->data.F32[2*i+1] = -yn->data.F32[i]*sin(phi);
-        graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+0]);
-        graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+0]);
-        graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+0]);
-        graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+0]);
-        graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+1]);
-        graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+1]);
-        graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+1]);
-        graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+1]);
+        xv->data.F32[2*i+0] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi) + zn->data.F32[i]*sin(theta)*sin(phi);
+        yv->data.F32[2*i+0] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi) + zn->data.F32[i]*cos(theta)*sin(phi);
+        zv->data.F32[2*i+0] = -yn->data.F32[i]*sin(phi)   + zn->data.F32[i]*cos(phi);
+        xv->data.F32[2*i+1] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi);
+        yv->data.F32[2*i+1] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi);
+        zv->data.F32[2*i+1] = -yn->data.F32[i]*sin(phi);
+        graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+0]);
+        graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+0]);
+        graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+0]);
+        graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+0]);
+        graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+1]);
+        graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+1]);
+        graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+1]);
+        graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+1]);
     }
     xv->n = xn->n;

branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.h

@@ -1 +1 @@
 /* @file  pmKapaPlots.h
- * @brief functions to make plots with the external program 'kapa' 
+ * @brief functions to make plots with the external program 'kapa'
  *
  * @author EAM, IfA
  *
- * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-11-08 01:52:34 $
+ * @version $Revision: 1.1.16.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-19 17:59:50 $
  * Copyright 2006 Institute for Astronomy, University of Hawaii
  */
@@ -15 +15 @@
 /// @{
 
-bool pmSourceSetVisual (bool mode);
 bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);