Index: branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.c (revision 21536) @@ -23,5 +23,4 @@ //variables to determine when things are plotted -static bool isVisual = false; static bool plotGridMatch = true; static bool plotTweak = true; @@ -48,10 +47,4 @@ /* Initialization Routines */ -bool pmAstromSetVisual (bool mode) { - isVisual = mode; - return true; -} - - bool pmAstromVisualClose() { @@ -68,9 +61,8 @@ { // 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,5 +214,5 @@ // pause and wait for user input: - pmVisualAskUser(&plotRawStars, &isVisual); + pmVisualAskUser(&plotRawStars); psFree (xVec); @@ -234,9 +226,8 @@ { //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,5 +275,5 @@ // pause and wait for user input: - pmVisualAskUser(&plotRefStars, &isVisual); + pmVisualAskUser(&plotRefStars); psFree (xVec); @@ -301,9 +292,8 @@ // 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,5 +367,5 @@ // pause and wait for user input: - pmVisualAskUser (&plotLumFunc, &isVisual); + pmVisualAskUser (&plotLumFunc); } return true; @@ -391,9 +381,8 @@ //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,5 +469,5 @@ //ask for user input and finish - pmVisualAskUser (&plotRemoveClumps, &isVisual); + pmVisualAskUser (&plotRemoveClumps); psFree (xVec); psFree (yVec); @@ -496,9 +485,8 @@ //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,5 +499,5 @@ //ask for user input and finish - pmVisualAskUser(&plotOneChipFit, &isVisual); + pmVisualAskUser(&plotOneChipFit); return true; } @@ -522,8 +510,7 @@ { //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,5 +575,5 @@ KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y"); - pmVisualAskUser(&plotFixChips, &isVisual); + pmVisualAskUser(&plotFixChips); psFree(xNew); psFree(yNew); @@ -606,9 +593,8 @@ //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,5 +672,5 @@ psFree(xResid); psFree(yResid); - pmVisualAskUser (&plotAstromGuessCheck, &isVisual); + pmVisualAskUser (&plotAstromGuessCheck); return true; } @@ -696,8 +682,7 @@ { //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,5 +721,5 @@ KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP); - pmVisualAskUser (&plotCommonScale, &isVisual); + pmVisualAskUser (&plotCommonScale); psFree(xVec); @@ -750,8 +735,7 @@ //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,10 +743,10 @@ //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,5 +759,5 @@ { //make sure we want to plot this - if (!isVisual || !plotMosaicMatches) return true; + if (!pmVisualIsVisual() || !plotMosaicMatches) return true; char title[60]; @@ -782,15 +766,14 @@ 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,7 +790,6 @@ { //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,5 +934,5 @@ KapaPlotVector (kapa, 2, yslice, "y"); - pmVisualAskUser(&plotGridMatch, &isVisual); + pmVisualAskUser(&plotGridMatch); psFree(dXplot); psFree(dYplot); @@ -966,7 +948,6 @@ { //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,5 +1032,5 @@ KAPA_LABEL_XP); - pmVisualAskUser(&plotTweak, &isVisual); + pmVisualAskUser(&plotTweak); psFree(xIndices); @@ -1341,5 +1322,4 @@ # 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; } Index: branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/astrom/pmAstrometryVisual.h (revision 21536) @@ -22,10 +22,4 @@ 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); Index: branches/cnb_branch_20090215/psModules/src/camera/pmFPACopy.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/camera/pmFPACopy.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/camera/pmFPACopy.c (revision 21536) @@ -457,4 +457,31 @@ 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; Index: branches/cnb_branch_20090215/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/detrend/pmShutterCorrection.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/detrend/pmShutterCorrection.h (revision 21536) @@ -5,6 +5,6 @@ * @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 */ Index: branches/cnb_branch_20090215/psModules/src/extras/pmVisual.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/extras/pmVisual.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/extras/pmVisual.c (revision 21536) @@ -28,4 +28,23 @@ # 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,4 +53,5 @@ if (*kapid == -1) { fprintf (stderr, "Failure to open kapa.\n"); + isVisual = false; return false; } @@ -52,5 +72,5 @@ -bool pmVisualAskUser(bool *plotFlag, bool *packageFlag) +bool pmVisualAskUser(bool *plotFlag) { char key[10]; @@ -63,5 +83,5 @@ } if (key[0] == 'a') { - *packageFlag = false; + isVisual = false; } return true; @@ -304,9 +324,11 @@ #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;} Index: branches/cnb_branch_20090215/psModules/src/extras/pmVisual.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/extras/pmVisual.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/extras/pmVisual.h (revision 21536) @@ -10,4 +10,22 @@ #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,8 +50,6 @@ * 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,7 +131,10 @@ // 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); Index: branches/cnb_branch_20090215/psModules/src/imcombine/Makefile.am =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/Makefile.am (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/Makefile.am (revision 21536) @@ -19,5 +19,6 @@ pmSubtractionThreads.c \ pmPSFEnvelope.c \ - pmSubtractionVisual.c + pmSubtractionVisual.c \ + pmStackVisual.c pkginclude_HEADERS = \ @@ -37,5 +38,6 @@ pmSubtractionThreads.h \ pmPSFEnvelope.h \ - pmSubtractionVisual.h + pmSubtractionVisual.h \ + pmStackVisual.h CLEANFILES = *~ Index: branches/cnb_branch_20090215/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/pmPSFEnvelope.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/pmPSFEnvelope.c (revision 21536) @@ -30,8 +30,8 @@ #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,4 +162,5 @@ 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,4 +180,5 @@ // 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,4 +233,5 @@ // 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,4 +315,5 @@ { 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,4 +322,5 @@ { psFits *fits = psFitsOpen("psf_field_resid.fits", "w"); + pmStackVisualPlotTestImage(generated->image, "psf_field_resid.fits"); psFitsWriteImage(fits, NULL, generated->image, 0, NULL); psFitsClose(fits); Index: branches/cnb_branch_20090215/psModules/src/imcombine/pmStack.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/pmStack.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/pmStack.c (revision 21536) @@ -8,6 +8,6 @@ * @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,7 +33,8 @@ #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 Index: branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionIO.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionIO.c (revision 21536) @@ -16,4 +16,5 @@ #include "pmSubtractionIO.h" +#include "pmStackVisual.h" #define ARRAY_BUFFER 16 // Number to add to array at a time @@ -144,6 +145,6 @@ // 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,4 +168,6 @@ 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."); Index: branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.c (revision 21536) @@ -28,5 +28,4 @@ //variables to determine when things are plotted -static bool isVisual = false; static bool plotConvKernels = true; static bool plotStamps = true; @@ -42,9 +41,4 @@ // Initialization Routines -/** start or stop plotting */ -bool pmSubtractionSetVisual (bool mode) { - isVisual = mode; - return true; -} @@ -65,11 +59,10 @@ * @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,11 +73,9 @@ @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,5 +136,5 @@ pmVisualScaleImage(kapa, canvas, "Subtraction_Stamps", 0, true); - pmVisualAskUser(&plotStamps, &isVisual); + pmVisualAskUser(&plotStamps); return true; } @@ -152,7 +143,6 @@ 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,5 +199,5 @@ pmVisualScaleImage(kapa, canvas32, "Least_Squares", 0, true); - pmVisualAskUser(&plotLeastSquares, &isVisual); + pmVisualAskUser(&plotLeastSquares); psFree(canvas); psFree(canvas32); @@ -216,7 +206,6 @@ 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,5 +214,5 @@ pmVisualScaleImage(kapa, ref, "Reference", 1, true); pmVisualScaleImage(kapa, sub, "Subtraction", 2, true); - pmVisualAskUser(&plotImage, &isVisual); + pmVisualAskUser(&plotImage); return true; } @@ -267,5 +256,4 @@ #else -bool pmSubtractionSetVisual (bool mode) {return true;} bool pmSubtractionVisualClose() {return true;} bool pmSubtractionVisualPlotConvKernels(psImage *convKernels) {return true;} Index: branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/imcombine/pmSubtractionVisual.h (revision 21536) @@ -2,5 +2,4 @@ #define PM_SUBTRACTION_VISUAL_H -bool pmSubtractionSetVisual (bool mode); bool pmSubtractionVisualClose(); bool pmSubtractionVisualPlotConvKernels(psImage *convKernels); Index: branches/cnb_branch_20090215/psModules/src/objects/Makefile.am =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/Makefile.am (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/Makefile.am (revision 21536) @@ -36,4 +36,5 @@ pmSourceIO_PS1_DEV_1.c \ pmSourceIO_PS1_CAL_0.c \ + pmSourceIO_CMF_PS1_V1.c \ pmSourcePlots.c \ pmSourcePlotPSFModel.c \ Index: branches/cnb_branch_20090215/psModules/src/objects/pmModel.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmModel.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmModel.c (revision 21536) @@ -6,6 +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,4 +59,5 @@ tmp->chisqNorm = 0.0; tmp->nDOF = 0; + tmp->nPix = 0; tmp->nIter = 0; tmp->radiusFit = 0; @@ -234,14 +235,14 @@ 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,5 +256,5 @@ // 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,56 +276,56 @@ 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) { Index: branches/cnb_branch_20090215/psModules/src/objects/pmModel.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmModel.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmModel.h (revision 21536) @@ -5,6 +5,6 @@ * @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,4 +92,5 @@ 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,5 +151,5 @@ 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,5 +169,5 @@ 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 ); Index: branches/cnb_branch_20090215/psModules/src/objects/pmPeaks.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmPeaks.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmPeaks.c (revision 21536) @@ -6,6 +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,5 +28,5 @@ 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,15 +63,18 @@ 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,5 +89,5 @@ 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,7 +248,7 @@ 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,5 +256,5 @@ *****************************************************************************/ psVector *pmPeaksInVector(const psVector *vector, - psF32 threshold) + psF32 threshold) { psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__); @@ -347,12 +350,12 @@ 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 Index: branches/cnb_branch_20090215/psModules/src/objects/pmSource.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSource.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSource.c (revision 21536) @@ -6,6 +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,4 +113,5 @@ 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,4 +535,6 @@ 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,14 +543,18 @@ 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,4 +1056,19 @@ } +// 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,5 +1088,4 @@ 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,5 +1109,4 @@ 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; Index: branches/cnb_branch_20090215/psModules/src/objects/pmSource.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSource.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSource.h (revision 21536) @@ -3,6 +3,6 @@ * @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,23 +35,41 @@ } 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,6 +86,6 @@ */ 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,4 +101,5 @@ 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,6 +259,7 @@ 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); Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceFitModel.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceFitModel.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceFitModel.c (revision 21536) @@ -6,6 +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,4 +184,5 @@ model->chisq = myMin->value; model->nIter = myMin->iter; + model->nPix = y->n; model->nDOF = y->n - nParams; model->flags |= PM_MODEL_STATUS_FITTED; Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.c (revision 21536) @@ -3,6 +3,6 @@ * @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,8 +488,11 @@ } 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,4 +503,7 @@ 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,4 +513,7 @@ 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,6 +938,6 @@ 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); } } Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO.h (revision 21536) @@ -4,6 +4,6 @@ * @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,11 +26,16 @@ 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,4 +47,5 @@ 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,4 +73,6 @@ bool pmReadoutCheckDataStatusForSources (const pmReadout *readout); +bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos); + /// @} # endif /* PM_SOURCE_IO_H */ Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 21536) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 21536) @@ -0,0 +1,675 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.3.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-19 17:59:50 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmResiduals.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmModelClass.h" +#include "pmSourceIO.h" + +// panstarrs-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data + +bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, + psMetadata *imageHeader, psMetadata *tableHeader, char *extname) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(extname, false); + + psArray *table; + psMetadata *row; + int i; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + psF32 errMag, chisq, apRadius; + psS32 nPix, nDOF; + + pmChip *chip = readout->parent->parent; + pmFPA *fpa = chip->parent; + + bool status1 = false; + bool status2 = false; + float magOffset = NAN; + float exptime = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE"); + float zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); + float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR"); + if (status1 && status2 && (exptime > 0.0)) { + magOffset = zeropt + 2.5*log10(exptime); + } + + // if the sequence is defined, write these in seq order; otherwise + // write them in S/N order: + if (sources->n > 0) { + pmSource *source = (pmSource *) sources->data[0]; + if (source->seq == -1) { + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state + // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // If source->seq is -1, source was generated in this analysis. If source->seq is + // not -1, source was read from elsewhere: in the latter case, preserve the source + // ID. source.seq is used instead of source.id since the latter is a const + // generated on Alloc, and would thus be wrong for read in sources. + if (source->seq == -1) { + source->seq = i; + } + + // no difference between PSF and non-PSF model + pmModel *model = source->modelPSF; + + if (model != NULL) { + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // in linear-fit mode, there is no error on the centroid + xErr = source->peak->dx; + yErr = source->peak->dy; + } + 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; + nDOF = model->nDOF; + nPix = model->nPix; + apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + xErr = source->peak->dx; + yErr = source->peak->dy; + axes.major = NAN; + axes.minor = NAN; + axes.theta = NAN; + chisq = NAN; + nDOF = 0; + nPix = 0; + apRadius = NAN; + errMag = NAN; + } + + float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN; + float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN; + psS16 nImageOverlap = 1; + + psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + float posAngle = 0.0; + float pltScale = 0.0; + pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", xErr); // XXX this is only measured for non-linear fits + psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", yErr); // XXX this is only measured for non-linear fits + 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); + psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", posAngle*PS_DEG_RAD); + psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", pltScale*PS_DEG_RAD*3600.0); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", errMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", apRadius); + psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", peakMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", calMag); + psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); + psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", chisq); + psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeight); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", nDOF); + psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", nPix); + + // distinguish moments measure from window vs S/N > XX ?? + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", source->moments->Mxx); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", source->moments->Mxy); + psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", source->moments->Myy); + + psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + + // XXX not sure how to get this : need to load Nimages with weight? + psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap); + psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + + psArrayAdd (table, 100, row); + psFree (row); + + // 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; + } + + if (table->n == 0) { + psFitsWriteBlank (fits, tableHeader, extname); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTable (fits, tableHeader, table, extname)) { + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree(table); + return false; + } + psFree (table); + + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(header, false); + + bool status; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + + // define PSF model type + // XXX need to carry the extra model parameters + int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + assert (modelType > -1); + + // We get the size of the table, and allocate the array of sources first because the table + // is large and ephemeral --- when the table gets blown away, whatever is allocated after + // the table is read blocks the free. In fact, it's better to read the table row by row. + long numSources = psFitsTableSize(fits); // Number of sources in table + psArray *sources = psArrayAlloc(numSources); // Array of sources, to return + + // convert the table to the pmSource entriesa + for (int i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags + + // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); + axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); + axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); + axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); + + PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); + dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); + source->sky = PAR[PM_PAR_SKY]; + source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); + source->errMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); + source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + 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); + + float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); + float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + source->peak->flux = peakFlux; + source->peak->dx = dPAR[PM_PAR_XPOS]; + source->peak->dy = dPAR[PM_PAR_YPOS]; + + source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF"); + source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); + source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); + model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); + model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); + model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); + + source->moments = pmMomentsAlloc (); + source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); + source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); + source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); + + source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + assert (status); + + sources->data[i] = source; + psFree(row); + } + + return sources; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + + bool status; + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psF32 xPos, yPos; + psF32 xErr, yErr; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + table = psArrayAllocEmpty (sources->n); + + // which extended source analyses should we perform? + bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + + psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + assert (radialBinsLower->n == radialBinsUpper->n); + + // 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]; + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); + psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); + psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); + 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); + } + + if (table->n == 0) { + 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; + } + psFree (outhead); + psFree (table); + + return true; +} + +// XXX this layout is still the same as PS1_DEV_1 +bool pmSourcesWrite_CMF_PS1_V1_XFIT (psFits *fits, psArray *sources, char *extname) +{ + + psArray *table; + psMetadata *row; + psF32 *PAR, *dPAR; + psEllipseAxes axes; + psF32 xPos, yPos; + psF32 xErr, yErr; + char name[64]; + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortBySN); + + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + 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); + } + } + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + 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); + } + } + + if (table->n == 0) { + 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; + } + psFree (outhead); + psFree (table); + return true; +} + +bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) { + + pmFPA *fpa = chip->parent; + + if (!chip->toFPA) goto escape; + if (!fpa->toTPA) goto escape; + if (!fpa->toSky) goto escape; + + // generate RA,DEC + psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y; + + // calculate the astrometry for the coordinate of interest + ptCH.x = xPos; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP); + psDeproject (ptSky, &ptTP_o, fpa->toSky); + + // calculate the astrometry for the coordinate + 1pix in X + ptCH.x = xPos + 1.0; + ptCH.y = yPos; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP); + + // calculate the astrometry for the coordinate + 1pix in Y + ptCH.x = xPos; + ptCH.y = yPos + 1.0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP); + + // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane + // degrees + + float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x); + float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y); + + float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x); + float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y); + + float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx); + float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy); + *pltScale = 0.5*(pltScale_x + pltScale_y); + + float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx); + float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy); + *posAngle = 0.5*(posAngle_x + posAngle_y); + + return true; + +escape: + // no astrometry calibration, give up + ptSky->r = NAN; + ptSky->d = NAN; + *posAngle = NAN; + *pltScale = NAN; + + return false; +} Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 21536) @@ -3,6 +3,6 @@ * @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,6 +49,6 @@ // 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,5 +68,5 @@ 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,7 +77,7 @@ 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,7 +90,7 @@ 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,15 +104,15 @@ 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,26 +123,26 @@ 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,5 +150,5 @@ 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,9 +170,9 @@ 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,6 +277,6 @@ 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,5 +289,5 @@ 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,5 +321,5 @@ 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,8 +331,8 @@ 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,41 +360,41 @@ // 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,113 +410,113 @@ 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,5 +540,5 @@ 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,7 +549,7 @@ 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,11 +566,11 @@ 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,103 +580,103 @@ 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); Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 21536) @@ -3,6 +3,6 @@ * @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,6 +49,6 @@ // 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,7 +73,7 @@ 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,15 +87,15 @@ 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,6 +106,6 @@ axes.minor = NAN; axes.theta = NAN; - chisq = NAN; - errMag = NAN; + chisq = NAN; + errMag = NAN; } @@ -130,11 +130,11 @@ 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,6 +247,6 @@ 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,5 +259,5 @@ 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,26 +311,26 @@ // 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,101 +343,101 @@ 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,11 +470,11 @@ 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,85 +484,85 @@ 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); Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.c (revision 21536) @@ -3,6 +3,6 @@ * @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,5 +145,5 @@ // 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,5 +164,5 @@ // 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,5 +170,4 @@ // 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,57 +185,4 @@ } - // 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,5 +215,5 @@ // 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,8 +290,7 @@ // 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,13 +322,13 @@ 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,4 +557,5 @@ # endif +// determine chisq, etc for linear normalization-only fit bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal) @@ -632,4 +578,5 @@ } } + model->nPix = Npix; model->nDOF = Npix - 1; model->chisq = dC; @@ -817,2 +764,54 @@ } +// 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 Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourcePhotometry.h (revision 21536) @@ -4,6 +4,6 @@ * @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,10 +46,10 @@ 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); Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.c =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.c (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.c (revision 21536) @@ -9,9 +9,9 @@ # if (HAVE_KAPA) # include +#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,26 +20,21 @@ 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,11 +50,11 @@ 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,9 +65,9 @@ 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,5 +150,5 @@ fprintf (stdout, "[c]ontinue? "); if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); + psWarning("Unable to read option"); } return true; @@ -173,18 +168,18 @@ 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; Index: branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.h =================================================================== --- branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.h (revision 21495) +++ branches/cnb_branch_20090215/psModules/src/objects/pmSourceVisual.h (revision 21536) @@ -1,9 +1,9 @@ /* @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,5 +15,4 @@ /// @{ -bool pmSourceSetVisual (bool mode); bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);