Index: /trunk/psModules/src/astrom/pmAstrometryDistortion.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryDistortion.c (revision 41892) @@ -29,4 +29,5 @@ #include "pmAstrometryRegions.h" #include "pmAstrometryDistortion.h" +#include "pmAstrometryUtils.h" #include "pmKapaPlots.h" @@ -183,4 +184,6 @@ PS_ASSERT_PTR_NON_NULL(fpa, false); PS_ASSERT_ARRAY_NON_NULL(gradients, false); + + int ExtraOrders = pmAstrometryGetExtraOrders(); psPolynomial2D *localX = NULL; @@ -340,12 +343,7 @@ psRegion *region = pmAstromFPAExtent (fpa); - // XXX psFree (fpa->fromTPA); - // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4); - // XXX fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, *region, 50); - // XXX psFree (myPT); - // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear) psFree (fpa->fromTPA); - fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 100, 6); + fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50, ExtraOrders); psFree (region); Index: /trunk/psModules/src/astrom/pmAstrometryModel.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 41892) @@ -461,4 +461,6 @@ bool status; + int ExtraOrders = pmAstrometryGetExtraOrders(); + // set FITS cursor if (!psFitsMoveExtName (file->fits, "CHIPS")) { @@ -498,5 +500,5 @@ int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER"); if (chip->toFPA == NULL) { - chip->toFPA = psPlaneTransformAlloc(nX, nY); + chip->toFPA = psPlaneTransformAlloc(nX, nY, PS_POLYNOMIAL_ORD); // chip->fpa uses ordinary poly } else { REQUIRE (chip->toFPA->x->nX == nX, "mismatch in chip order"); @@ -525,6 +527,5 @@ // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear) psFree (chip->fromFPA); - chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6); - + chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders); psFree (region); } @@ -539,4 +540,6 @@ bool status; + + int ExtraOrders = pmAstrometryGetExtraOrders(); if (!psFitsMoveExtName (file->fits, "FP")) { @@ -565,5 +568,5 @@ if (file->fpa->toTPA == NULL) { // allocate the new transformation - file->fpa->toTPA = psPlaneTransformAlloc(nX, nY); + file->fpa->toTPA = psPlaneTransformAlloc(nX, nY, PS_POLYNOMIAL_ORD); // fpa->tpa uses ORD } else { REQUIRE (file->fpa->toTPA->x->nX == nX, "mismatch in chip order"); @@ -585,12 +588,7 @@ psRegion *region = pmAstromFPAExtent (file->fpa); - // XXX psFree (file->fpa->fromTPA); - // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(file->fpa->toTPA->x->nX+4, file->fpa->toTPA->x->nY+4); - // XXX file->fpa->fromTPA = psPlaneTransformInvert(myPT, file->fpa->toTPA, *region, 50); - // XXX psFree (myPT); - // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear) psFree (file->fpa->fromTPA); - file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 100, 6); + file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50, ExtraOrders); psFree (model); @@ -689,4 +687,6 @@ bool status; + int ExtraOrders = pmAstrometryGetExtraOrders(); + // these externally supplied values are used to set the final transformation terms double RA = psMetadataLookupF64 (&status, concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA"); @@ -728,8 +728,4 @@ psLogMsg ("psModules.astrom", 4, "Position Angle: %f, Model Position Angle Zero Point: %f\n", POS, PosZero); -// psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, rotatorParity * (PosZero - POS)); -// psFree (file->fpa->fromTPA); -// file->fpa->fromTPA = fromTPA; - psPlaneTransform *toTPA = psPlaneTransformRotate (NULL, file->fpa->toTPA, rotatorParity * (POS - PosZero)); psFree (file->fpa->toTPA); @@ -739,11 +735,6 @@ psRegion *region = pmAstromFPAExtent (file->fpa); - // XXX psFree (file->fpa->fromTPA); - // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(file->fpa->toTPA->x->nX+4, file->fpa->toTPA->x->nY+4); - // XXX file->fpa->fromTPA = psPlaneTransformInvert(myPT, file->fpa->toTPA, *region, 50); - // XXX psFree (myPT); - psFree (file->fpa->fromTPA); - file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region,100, 10); + file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50, ExtraOrders); psFree (region); Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 41892) @@ -40,4 +40,5 @@ // XXX this is defined in pmPSFtry.h, which makes no sense float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction); +float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value); #define PM_ASTROMETRYOBJECTS_DEBUG 1 @@ -339,4 +340,11 @@ psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32); + // we measure the stdev of the median residual in NxN bins. + // use only valid (not NAN) measurements + psVector *xPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *yPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *xResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32); + psVector *yResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32); + for (int i = 0; i < match->n; i++) { if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; @@ -344,4 +352,18 @@ pmAstromObj *rawStar = raw->data[pair->raw]; if (!isfinite(rawStar->dMag)) continue; + + bool isValid = true; + isValid = isValid & isfinite (x->data.F32[i]); + isValid = isValid & isfinite (y->data.F32[i]); + isValid = isValid & isfinite (xRes->data.F32[i]); + isValid = isValid & isfinite (yRes->data.F32[i]); + if (isValid) { + psVectorAppend (xPosValid, x->data.F32[i]); + psVectorAppend (yPosValid, y->data.F32[i]); + psVectorAppend (xResValid, xRes->data.F32[i]); + psVectorAppend (yResValid, yRes->data.F32[i]); + } + + // for the systematic error, use only high S/N stars if (rawStar->dMag > 0.02) continue; @@ -369,4 +391,5 @@ results->dXsys = results->dYsys = NAN; results->dXrange = results->dYrange = NAN; + results->dXstdev = results->dYstdev = NAN; } else { results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05); @@ -375,8 +398,11 @@ results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev); results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev); - } - - psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange); - psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange); + + results->dXstdev = pmAstrom2DSystematics (xPosValid, yPosValid, xResValid); + results->dYstdev = pmAstrom2DSystematics (xPosValid, yPosValid, yResValid); + } + + psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f, dXstdev: %f\n", results->dXsys, results->dXrange, results->dXstdev); + psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f, dYstdev: %f\n", results->dYsys, results->dYrange, results->dYstdev); psFree (xErr); @@ -386,4 +412,8 @@ psFree (xResGood); psFree (yResGood); + psFree (xPosValid); + psFree (yPosValid); + psFree (xResValid); + psFree (yResValid); psFree (x); @@ -395,4 +425,97 @@ return (results); +} + +# define VAL_COUNT 9.0 +# define MIN_COUNT 5.0 + +float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value) { + + // pre-filter the values to ensure no NANs, other invalid + if (xPos->n < VAL_COUNT*2*2) { + return NAN; + } + + // generate a grid covering the full range of x,y and measure the median value in each + // grid cell. calculate the stdev of those median values. + + // VAL_COUNT (9) is the (min) goal density, but a single cell may have only MIN_COUNT (7) + // we have N points. require a min of 9 pts per cell (configurable?). grid is square. + // Ncell*Ncell*9 = Npts, Ncell = MIN(sqrt(Npts/9), 5) + + int Ncell = PS_MIN(sqrt((float)(xPos->n / VAL_COUNT)), 2); // have to at least have a 2x2 grid + + // find the range of x,y values + float xMin = 1e9, xMax = -1e9, yMin = 1e9, yMax = 1e9; + for (int i = 0; i < xPos->n; i++) { + xMin = PS_MIN(xPos->data.F32[i],xMin); + xMax = PS_MAX(xPos->data.F32[i],xMax); + yMin = PS_MIN(yPos->data.F32[i],yMin); + yMax = PS_MAX(yPos->data.F32[i],yMax); + } + + float xStep = Ncell / (xMax - xMin); + float yStep = Ncell / (yMax - yMin); + + if (isnan(xStep)) { + return NAN; + } + if (isnan(yStep)) { + return NAN; + } + + psArray *xBin = psArrayAlloc (Ncell); + for (int ix = 0; ix < Ncell; ix++) { + psArray *yBin = psArrayAlloc (Ncell); + xBin->data[ix] = yBin; + for (int iy = 0; iy < Ncell; iy++) { + yBin->data[iy] = psVectorAllocEmpty(128, PS_TYPE_F32); + } + } + + // xValue = ix/xStep + xMin -> ix = (xValue - xMin) * xStep; + + for (int i = 0; i < xPos->n; i++) { + int ix = PS_MIN(PS_MAX(0, (xPos->data.F32[i] - xMin) * xStep), Ncell - 1); + int iy = PS_MIN(PS_MAX(0, (yPos->data.F32[i] - yMin) * yStep), Ncell - 1); + + psArray *yBin = xBin->data[ix]; + psVector *Bin = yBin->data[iy]; + + psVectorAppend(Bin, value->data.F32[i]); + } + + psVector *sample = psVectorAllocEmpty (Ncell*Ncell, PS_TYPE_F32); + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + + // calculate the median for each vector and save on a vector + for (int ix = 0; ix < Ncell; ix++) { + psArray *yBin = xBin->data[ix]; + for (int iy = 0; iy < Ncell; iy++) { + psVector *Bin = yBin->data[iy]; + if (Bin->n < MIN_COUNT) { continue; } + + // psVectorStats resets stats so we can call this repeatedly + psVectorStats (stats, Bin, NULL, NULL, 0); + if (isfinite(stats->sampleMedian)) { + psVectorAppend (sample, stats->sampleMedian); + } + } + } + + stats->options = PS_STAT_SAMPLE_STDEV; + psVectorStats (stats, sample, NULL, NULL, 0); + float result = stats->sampleStdev; + + if (!isfinite(stats->sampleStdev)) { + fprintf (stderr, "*** bad solution ***\n"); + } + + // NOTE: the elements of xBin are freed automatically, which extends down to the vectors + psFree (xBin); + psFree (stats); + psFree (sample); + + return result; } Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 41892) @@ -102,4 +102,6 @@ double dXrange; ///< 10% - 90% range X residuals (unmasked, high S/N) double dYrange; ///< 10% - 90% range Y residuals (unmasked, high S/N) + double dXstdev; ///< stdev of median residual in NxN bins + double dYstdev; ///< stdev of median residual in NxN bins } pmAstromFitResults; Index: /trunk/psModules/src/astrom/pmAstrometryRegions.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryRegions.h (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryRegions.h (revision 41892) @@ -10,4 +10,7 @@ #ifndef PM_ASTROMETRY_REGIONS_H #define PM_ASTROMETRY_REGIONS_H + +// used by pmAstrometryDistortion.c, pmAstrometryWCS.c, pmAstrometryModel.c +// # define EXTRA_ORDERS 6 /// @addtogroup Astrometry Index: /trunk/psModules/src/astrom/pmAstrometryUtils.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryUtils.c (revision 41892) @@ -23,4 +23,15 @@ #include "pmAstrometryUtils.h" +static int transform_extra_orders = 3; + +int pmAstrometryGetExtraOrders (void) { + return transform_extra_orders; +} + +int pmAstrometrySetExtraOrders (int orders) { + transform_extra_orders = orders; + return transform_extra_orders; +} + // this is used by the test output block static int Nout = 0; @@ -120,11 +131,15 @@ // convert a transformation L(x,y) to L'(x-xo,y-yo) +// is this used for an upward (e.g., chip->fpa) or a downward (fpa->chip) transform? psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo) { - // validate fit order + // validate fit type: + // polynomial in input transforms must match type -- and for now be ORD + psAssert (input->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); + psAssert (input->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); if (output == NULL) { - output = psPlaneTransformAlloc(input->x->nX, input->x->nY); + output = psPlaneTransformAlloc(input->x->nX, input->x->nY, PS_POLYNOMIAL_ORD); } @@ -169,5 +184,5 @@ // rotate a transformation L(x,y) by theta -psPlaneTransform *psPlaneTransformRotate (psPlaneTransform *output, psPlaneTransform *input, double theta) +psPlaneTransform *psPlaneTransformRotate (psPlaneTransform *output, psPlaneTransform *input, double theta) { /* given the polynomial transformations: @@ -181,5 +196,6 @@ if (output == NULL) { - output = psPlaneTransformAlloc(input->x->nX, input->x->nY); + // generate a new transform using the same order and type as the input + output = psPlaneTransformAlloc(input->x->nX, input->x->nY, input->x->type); } @@ -216,5 +232,5 @@ // all coeffs and masks initially set to 0 - transform = psPlaneTransformAlloc (order, order); + transform = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD); for (int i = 0; i <= order; i++) { Index: /trunk/psModules/src/astrom/pmAstrometryUtils.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryUtils.h (revision 41892) @@ -25,4 +25,7 @@ bool psPlaneDistortIsDiagonal (psPlaneDistort *distort); +int pmAstrometryGetExtraOrders (void); +int pmAstrometrySetExtraOrders (int orders); + /// @} #endif Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 41891) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 41892) @@ -501,4 +501,6 @@ psPlaneTransform *toFPA; + int ExtraOrders = pmAstrometryGetExtraOrders(); + // create transformation with 0,0 reference pixel and units of degrees/pixel toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); @@ -552,5 +554,5 @@ // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter // XXX this only works if toTPA is at most a linear transformation - psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY); + psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY, PS_POLYNOMIAL_ORD); for (int i = 0; i <= toFPA->x->nX; i++) { for (int j = 0; j <= toFPA->x->nY; j++) { @@ -608,5 +610,5 @@ // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear) psFree (chip->fromFPA); - chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6); + chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders); psFree (region); @@ -648,9 +650,11 @@ */ + int ExtraOrders = pmAstrometryGetExtraOrders(); + psFree (chip->toFPA); if ((fabs(wcs->crpix1) > 0.01) || (fabs(wcs->crpix2) > 0.01)) { chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2); } else { - chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY); + chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY, PS_POLYNOMIAL_ORD); // copy the toFPA x,y, transformations to the wcs version @@ -668,5 +672,5 @@ // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear) psFree (chip->fromFPA); - chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6); + chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders); psFree (region); @@ -700,5 +704,5 @@ // the region defines the FPA pixels covered by the tranformation psFree (fpa->fromTPA); - int additional_orders = 4; // This is the number of orders that should be added. + int additional_orders = pmAstrometryGetExtraOrders(); // This is the number of orders that should be added. bool status = false; int config_additional_orders = psMetadataLookupS32(&status,fpa->analysis, "ADDITIONAL_WCS_ORDERS"); @@ -706,5 +710,5 @@ additional_orders = config_additional_orders; } - fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 100, additional_orders); + fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50, additional_orders); return true; } @@ -739,5 +743,5 @@ // XXX require fpa->toTPA->nX == 1 - psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY); + psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY, PS_POLYNOMIAL_ORD); for (int i = 0; i <= toTPA->x->nX; i++) { @@ -958,5 +962,5 @@ } - psPlaneTransform *newTrans = psPlaneTransformAlloc(1, 1); + psPlaneTransform *newTrans = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD); if (!psPlaneTransformFit(newTrans, src, dst, 0, 0)) { @@ -997,4 +1001,6 @@ PS_ASSERT_PTR_NON_NULL(inChip, NULL); + int ExtraOrders = pmAstrometryGetExtraOrders(); + if (outFPA == NULL) { outFPA = inFPA; @@ -1032,5 +1038,5 @@ // NOTE: the extraOrders value (4) should be ignored since outToFPA is specified to be linear - psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 100, 6); + psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50, ExtraOrders); if (!outFromFPA) { psFree(outToFPA); @@ -1116,5 +1122,5 @@ } - psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1); + psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD); newToFPA->x->coeffMask[1][1] = 1; newToFPA->y->coeffMask[1][1] = 1; @@ -1148,5 +1154,5 @@ psFree(dst); - // this is a linear transformation + // this is a linear transformation, no extra orders are needed psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1, 0); if (!newFromFPA) { @@ -1187,5 +1193,6 @@ psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree); - wcs->trans = psPlaneTransformAlloc (nXorder, nYorder); + // note: WCS transforms are always defined as chip to sky + wcs->trans = psPlaneTransformAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD); wcs->toSky = NULL; wcs->wcsCDkeys = 0; Index: /trunk/psModules/src/camera/pmFPAConstruct.c =================================================================== --- /trunk/psModules/src/camera/pmFPAConstruct.c (revision 41891) +++ /trunk/psModules/src/camera/pmFPAConstruct.c (revision 41892) @@ -853,4 +853,8 @@ } else { const char *content = getContent(fileInfo, phdu->header, contents); // The chip type + if (!content) { + psError(PS_ERR_UNKNOWN, false, "Unable to find CONTENT entry in header"); + return false; + } int chipNum = -1; // Chip number Index: /trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.c (revision 41891) +++ /trunk/psModules/src/camera/pmFPAfile.c (revision 41892) @@ -551,4 +551,7 @@ if (!strcasecmp(type, "KH.CORRECT")) { return PM_FPA_FILE_KH_CORRECT; + } + if (!strcasecmp(type, "PATTERN.ROW.AMP")) { + return PM_FPA_FILE_PATTERN_ROW_AMP; } if (!strcasecmp(type, "SUBKERNEL")) { @@ -606,4 +609,6 @@ case PM_FPA_FILE_KH_CORRECT: return ("KH.CORRECT"); + case PM_FPA_FILE_PATTERN_ROW_AMP: + return ("PATTERN.ROW.AMP"); case PM_FPA_FILE_SUBKERNEL: return ("SUBKERNEL"); Index: /trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.h (revision 41891) +++ /trunk/psModules/src/camera/pmFPAfile.h (revision 41892) @@ -52,4 +52,5 @@ PM_FPA_FILE_SRCTEXT, PM_FPA_FILE_PATTERN, + PM_FPA_FILE_PATTERN_ROW_AMP, PM_FPA_FILE_LINEARITY, PM_FPA_FILE_EXPNUM, Index: /trunk/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 41891) +++ /trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 41892) @@ -152,4 +152,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: { pmHDU *hdu = pmFPAviewThisHDU(view, fpa); Index: /trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileIO.c (revision 41891) +++ /trunk/psModules/src/camera/pmFPAfileIO.c (revision 41892) @@ -145,4 +145,25 @@ } +// list all defined pmFPAfiles +bool pmFPAfileIOList (pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(config, false); + PS_ASSERT_PTR_NON_NULL(config->files, false); + + psMetadata *files = config->files; + + // attempt to perform all create, read, write, close operations + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (files, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + pmFPAfile *file = item->data.V; + + fprintf (stderr, "%s : %d %d %d : %d %d %d %d\n", file->name, file->type, file->mode, file->state, + file->fileLevel, file->dataLevel, file->freeLevel, file->mosaicLevel); + } + psFree (iter); + return true; +} + // read the file, if necessary and possible bool pmFPAfileRead(pmFPAfile *file, const pmFPAview *view, pmConfig *config) @@ -231,4 +252,7 @@ case PM_FPA_FILE_PATTERN: status = pmPatternRead(view, file, config); + break; + case PM_FPA_FILE_PATTERN_ROW_AMP: + status = pmPatternRowAmpRead(view, file, config); break; case PM_FPA_FILE_SX: @@ -342,4 +366,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -359,5 +384,5 @@ if (file->state & PM_FPA_STATE_INACTIVE) { - psTrace("psModules.camera", 6, "skip write for %s, files is inactive", file->name); + psTrace("psModules.camera", 6, "skip write for %s, file is inactive", file->name); return true; } @@ -430,4 +455,8 @@ return true; } + if (file->type == PM_FPA_FILE_PATTERN_ROW_AMP) { + psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name); + return true; + } // open the file if not yet opened @@ -526,4 +555,8 @@ case PM_FPA_FILE_KH_CORRECT: psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name); + break; + + case PM_FPA_FILE_PATTERN_ROW_AMP: + psError(PS_ERR_IO, true, "cannot write type PATTERN.ROW.AMP (%s)", file->name); break; @@ -603,4 +636,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -681,4 +715,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_EXPNUM: psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); @@ -844,4 +879,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -1048,4 +1084,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: /trunk/psModules/src/camera/pmFPAfileIO.h =================================================================== --- /trunk/psModules/src/camera/pmFPAfileIO.h (revision 41891) +++ /trunk/psModules/src/camera/pmFPAfileIO.h (revision 41892) @@ -55,4 +55,6 @@ bool pmFPAfileReadPHU (pmFPAfile *file, const pmFPAview *view, pmConfig *config); +bool pmFPAfileIOList (pmConfig *config); + /// @} # endif Index: /trunk/psModules/src/detrend/pmDetrendDB.c =================================================================== --- /trunk/psModules/src/detrend/pmDetrendDB.c (revision 41891) +++ /trunk/psModules/src/detrend/pmDetrendDB.c (revision 41892) @@ -112,4 +112,5 @@ DETREND_STRING_CASE(AUXMASK); DETREND_STRING_CASE(KH_CORRECT); + DETREND_STRING_CASE(PATTERN_ROW_AMP); default: return NULL; Index: /trunk/psModules/src/detrend/pmDetrendDB.h =================================================================== --- /trunk/psModules/src/detrend/pmDetrendDB.h (revision 41891) +++ /trunk/psModules/src/detrend/pmDetrendDB.h (revision 41892) @@ -41,4 +41,5 @@ PM_DETREND_TYPE_AUXMASK, PM_DETREND_TYPE_KH_CORRECT, + PM_DETREND_TYPE_PATTERN_ROW_AMP, } pmDetrendType; Index: /trunk/psModules/src/detrend/pmPattern.c =================================================================== --- /trunk/psModules/src/detrend/pmPattern.c (revision 41891) +++ /trunk/psModules/src/detrend/pmPattern.c (revision 41892) @@ -9,4 +9,107 @@ #define PATTERN_ROW_BKG_FIX 1 +/* some in-line notes: + + patternMaskRow sets the data value to NAN, inconsistent with new plan + + here is the outline of pmPatternRow + + * at this point, we have already done overscan subtraction, right? + + * measure stats on the full cell (MEDIAN, STDEV) + ** subsample? + ** if it fails, it masks the entire cell and set the value to NAN + + * calculate an upper and lower threshold (median +/- T * sigma) + * define a normalized x-coordinate ('index') : + ** see note below on chebys + + * each row is treated independently + * pixels are masked for the fit if they are out-of-range + or if they are already masked + + ** note that the clipping threshold will be larger if there + are pixels which have astronomical structures + + a possible better option would be to set the threshold based on the median + and a sigma calculated from Poisson stats (do we know the gain?) + + ** fit is allowed to proceed if even N+1 pixels exist, which is clearly too low + + ** Remaining pixels are fitted with clip-fit + + ** solution is subtracted from the data + (this is implemented with psPolynomial1DEvalVector) + perhaps faster if we fixed the order to 2 and hardwired the result + + * after each row is fitted, the intercept (A value) is fitted + as a function of the y-coordinate and the result is subtracted + + * the slope value is also fitted as a function of the column + and added back in -- I'm not sure I understand this step. + + ***************** + + ** what we calculate are related to chebychevs (domain is -1 : +1) + *** T0(x) = 1 + *** T1(x) = x + *** T2(x) = 2x^2 - 1 + + *** we calculate y = A + Bx + Cx^2 + + a_0 + a_1 x + a_2 (2x^2 - 1) = A + B + Cx^2 + + a_1 = B + a_0 - a_2 = A + 2 a_2 = C + + a_0 = A + C/2 + a_1 = B + a_2 = C/2 + + ***************** + + I have 3 goals in re-working the code: + + 1) improve overall speed + 2) improve reliability of the fit + 3) skip fit if we can + + Let's assume the signal in the cell is light + bias drift + + The bias drift has an amplitude of ~5 - 10 DN + + That makes a detectable source with ~N * a few counts (multiple pixels in a row) + + So, the effective flux is ~10 * 5 = 50 DN + for which sky level is this value - 3 sigma? + + 50 / sqrt(sky sigma^2 * effective area) + + area ~ 5pixels, sky sigma^2 = sky + + 10 * Npix / sqrt(sky * Npix) = 3 + + Npix = sky * (S/N)^2 / (peak^2) + + sky = Npix * peak^2 / SN^2 + + if (sky < Npix * peak^2 / SN^2), we should skip: + + Npix ~ 5 + peak ~ 10 + SN ~ 3 (or even less) + + sky < 5 * 100 / 9 = 55 or so + + + To address these in order: + + 1) speed: + * the analysis threaded is not threaded: thread across cells + + 2) + + */ // Mask a row as bad @@ -33,9 +136,47 @@ } +// Comparison and swap functions for sorting values directly +#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B]) +#define SORT_SWAP(TYPE,A,B) { \ + if (A != B) { \ + TYPE temp = sampleArray[A]; \ + sampleArray[A] = sampleArray[B]; \ + sampleArray[B] = temp; \ + } \ +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Measurement and application ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad); + + +bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad); + + +// XXX allow user choice of binned vs unbinned analysis? bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) { + + bool status = false; + if (true) { + status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + } else { + status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + } + return status; +} + +// USE_BACKGROUND_STDEV: if TRUE, the analysis will use the measured robust stdev to clip the out-of-range pixles +// if FALSE, the stdev will be estimated based on the Poisson statistics of the median (sky level). +# define USE_BACKGROUND_STDEV 0 + +bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh, psStatsOptions clipMean, psStatsOptions clipStdev, psImageMaskType maskVal, psImageMaskType maskBad) @@ -48,4 +189,11 @@ PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + psImage *image = ro->image; // Image to correct psImage *mask = ro->mask; // Mask for image @@ -55,20 +203,57 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { - psWarning("Unable to calculate statistics on readout; masking entire readout."); + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); psErrorClear(); psFree(stats); psFree(rng); - psImageInit(image, NAN); - if (mask) { - psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); - } - if (ro->variance) { - psImageInit(image, NAN); - } return true; } + +# if (USE_BACKGROUND_STDEV) float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data float background = stats->robustMedian; +# else + // the signal we are looking for is a small variation on top of the background. if + // the background is uniform with only read noise + sky noise, then the pixel-to-pixel + // stdev should only be due to known noise sources and predictable. If the + // pixel-to-pixel variations are from other features, then those variations will + // probably dominate the row-by-row bias variations. + + // instead of using the image pixel statistics to measure the stdev, lets assume only + // dark noise plus poisson sky noise. we are not carrying in the read noise, but it is + // fairly modest for GPC1 (~10 DN) + + // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have + // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17 + + // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2) + // so smaller by a factor of 1.4 than what we predict, which is not very large + + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + + // XXX retrieve noise and gain from 'concepts' if possible +# define READNOISE 10 /* arbitrary number */ + float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data + float background = stats->robustMedian; + + // if the noise from the background is too large + float significance = nominalAmplitude / sigma; + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + return true; + } + fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); +# endif + psFree(stats); psFree(rng); @@ -109,18 +294,35 @@ psVectorInit(yaxisMask, 0); #endif + + // we really need more than order + 1 points (= 4). + // this should be tunable, but let's try 5 - 10% + int validNmin = numCols * 0.1; + for (int y = 0; y < numRows; y++) { psVectorInit(clipMask, 0); data = psImageRow(data, image, y); int num = 0; // Number of good pixels + + // if the unmasked pixels only span a small range in x then we cannot fit the + // 2nd order polynomial variations very well. Require a minimum fractional range + float validXmin = +1; + float validXmax = -1; + + // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY) + // (x % 3) || for (int x = 0; x < numCols; x++) { - if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || - data->data.F32[x] < lower || data->data.F32[x] > upper) { - clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || + data->data.F32[x] < lower || data->data.F32[x] > upper) { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; } else { clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0; num++; + validXmin = PS_MIN(indices->data.F32[x], validXmin); + validXmax = PS_MAX(indices->data.F32[x], validXmax); } } - if (num < order + 1) { + + // XXX how much time is spent in the fitting + if (num < validNmin) { // Not enough points to fit patternMaskRow(ro, y, maskBad); @@ -131,4 +333,5 @@ continue; } + // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise? if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) { psWarning("Unable to fit polynomial to row %d", y); @@ -215,5 +418,5 @@ for (int x = 0; x < numCols; x++) { image->data.F32[y][x] += solution->data.F32[y]; - psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]); + psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]); } corr->data.F64[y][0] -= solution->data.F32[y]; @@ -241,5 +444,6 @@ for (int x = 0; x < numCols; x++) { image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x]; - psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]); + // XXX EAM : this was [x] which is wrong + psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]); } corr->data.F64[y][1] -= solution->data.F32[y] ; @@ -262,4 +466,339 @@ psFree(indices); + psFree(clip); + psFree(clipMask); + psFree(poly); + psFree(data); + + return true; +} + +# define NPIX 15 + +// bin by NPIX in the x-direction to reduce the number of calculations needed to measure +// the pattern +bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PS_ASSERT_INT_NONNEGATIVE(order, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + + psImage *image = ro->image; // Image to correct + psImage *mask = ro->mask; // Mask for image + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); + psErrorClear(); + psFree(stats); + psFree(rng); + + // EAM 20211011 : we used to mask cells which fail the above, but this seems excessive + // psImageInit(image, NAN); + // if (mask) { + // psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); + // } + // if (ro->variance) { + // psImageInit(image, NAN); + // } + return true; + } + + // if USE_BACKGROUND_STDEV is TRUE, the observed standard deviation is used to set the + // thresholds. this is going to be an overestimate if there is any structure in the + // image. If FALSE, the thresholds are set based on poisson stats for the background + // level. We assume the gain is 1, so this is an overestimate if the gain is > 1 + +# if (USE_BACKGROUND_STDEV) + float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data + float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data + float background = stats->robustMedian; +# else + // the signal we are looking for is a small variation on top of the background. if + // the background is uniform with only read noise + sky noise, then the pixel-to-pixel + // stdev should only be due to known noise sources and predictable. If the + // pixel-to-pixel variations are from other features, then those variations will + // probably dominate the row-by-row bias variations. + + // instead of using the image pixel statistics to measure the stdev, lets assume only + // dark noise plus poisson sky noise. we are not carrying in the read noise, but it is + // fairly modest for GPC1 (~10 DN) + + // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have + // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17 + + // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2) + // so smaller by a factor of 1.4 than what we predict, which is not very large + + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + + // retrieve noise and gain from 'concepts' if possible +# define READNOISE 10 /* arbitrary number */ + float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data + float background = stats->robustMedian; + + // if the noise from the background is too large + float significance = nominalAmplitude / sigma; + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psFree(stats); + psFree(rng); + return true; + } + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); +# endif + + psFree(stats); + psFree(rng); + + // the vector 'indices' maps the x-coordinate to a range [-1:1]. the element number (i) of indices + // related to the x-coordinate (column number) by x = (i + 0.5) * NPIX + + int nSamples = numCols / NPIX; + + psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fit solutions + psVector *xFit = psVectorAlloc(nSamples, PS_TYPE_F32); // x-coordinate for fitting + psVector *yFit = psVectorAlloc(nSamples, PS_TYPE_F32); // flux values for fitting + + // xFit elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX + + float norm = 2.0 / (float)numCols; // Normalisation for indices + for (int sample = 0; sample < nSamples; sample ++) { + int x = (sample + 0.5)*NPIX; + xFit->data.F32[sample] = x * norm - 1.0; + } + for (int x = 0; x < numCols; x ++) { + indices->data.F32[x] = x * norm - 1.0; + } + + psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics + clip->clipIter = iter; + clip->clipSigma = rej; + psVector *clipMask = psVectorAlloc(nSamples, PS_TYPE_VECTOR_MASK); // Mask for clipping + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit + psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit + + psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied + psImageInit(corr, NAN); + + // CZW: 2011-11-30 + // Define the vectors to hold the "x" and "y" slope trends. + // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term + // of the PATTERN.ROW fit between individual rows across the cell. Similarly, the 1-st + // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we + // are fitting with PATTERN.ROW in the first place). However, the thing we're trying to + // fix with PATTERN.ROW is the detector level bias wiggles. These should be overlaid on + // the true sky level. Therefore, simply applying the PATTERN.ROW correction will + // introduce cell-to-cell sky variations as these two trends are removed. To avoid this, + // We store the 0th and 1st order values used for each row, and then fit a polynomial to + // these results. By re-adding these systematic trends back, we can remove the row-to-row + // variations without improperly removing the real sky trend. + psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term + psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit + psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term + psVectorInit(yaxisMask, 0); + + // validNmin is the minimum number of samples needed to measure the trend. + // this should be tunable, but let's try 5 - 10% + int validNmin = PS_MAX (nSamples * 0.1, order + 2); + + for (int y = 0; y < numRows; y++) { + psVectorInit(clipMask, 0); + data = psImageRow(data, image, y); + int num = 0; // Number of good pixels + + // if the unmasked pixels only span a small range in x then we cannot fit the + // 2nd order polynomial variations very well. Require a minimum fractional range + float validXmin = +1; + float validXmax = -1; + + for (int sample = 0; sample < nSamples; sample ++) { + + // store valid samples in the array to be sorted + float sampleArray[NPIX]; + int seq = 0; + for (int j = 0; j < NPIX; j++) { + int pix = sample * NPIX + j; // real pixel elements in x-dir + psAssert (pix >= 0, "invalid pix value"); + psAssert (pix < numCols, "invalid pix value"); + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][pix] & maskVal)) continue; + if (data->data.F32[pix] < lower || data->data.F32[pix] > upper) continue; + sampleArray[seq] = data->data.F32[pix]; // store the value to be sorted + seq ++; + } + if (seq < 1) { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF; + yFit->data.F32[sample] = NAN; + continue; + } + // note that we are treating the x-coordinate as the center + // of the binned pixel group, even if some or most pixels have + // been masked. compared to the amplitude of the slope, this + // error is small + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0; + validXmin = PS_MIN(xFit->data.F32[sample], validXmin); + validXmax = PS_MAX(xFit->data.F32[sample], validXmax); + num++; + + // PSSORT operates on sampleArray (see define of macro SORT_SWAP above) + PSSORT (seq, SORT_COMPARE, SORT_SWAP, float); + + int midPt = 0.5 * seq; + if (seq % 2 == 1) { psAssert (midPt >= 0, "invalid midPt"); } + if (seq % 2 == 0) { psAssert (midPt >= 1, "invalid midPt"); } + psAssert (midPt < NPIX, "invalid midPt"); + + float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]); + yFit->data.F32[sample] = medValue; + } + + // If not enough points are valid, skip + if (num < validNmin) { + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise? + if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, yFit, NULL, xFit)) { + psWarning("Unable to fit polynomial to row %d", y); + psErrorClear(); + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + // Store the results we found for this row. + yaxisData->data.F32[y] = poly->coeff[0]; + xaxisData->data.F32[y] = poly->coeff[1]; + psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]); + + memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector + if (!solution) { + psWarning("Unable to evaluate polynomial for row %d", y); + psErrorClear(); + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + + psAssert (solution->n == numCols, "oops"); + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] -= solution->data.F32[x]; + psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]); + } + psFree(solution); + } + + // Put the global trends back that were removed by the PATTERN.ROW correction. + // Set up the indices for the polynomial + psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32); + norm = 2.0 / (float)numRows; + for (int y = 0; y < numRows; y++) { + yaxisIndices->data.F32[y] = y * norm - 1.0; + psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]); + } + + // Fit the trend of the constant term, producing the y-axis global trend + psStatsInit(clip); + psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to y-axis trend"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + } + corr->data.F64[y][0] -= background; + } + } else { + psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + } + corr->data.F64[y][0] -= background; + } + } else { + psAssert (solution->n == numRows, "oops"); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y]; + // XXX EAM : this was [x], which is wrong + psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]); + } + corr->data.F64[y][0] -= solution->data.F32[y]; + } + } + psFree(solution); + } + + // Fit the trend of the linear term, producing the x-axis global trend + // We can use the same mask vector, as the same rows failed the row-fit earlier. + psStatsInit(clip); + psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to x-axis trend"); + psErrorClear(); + } else { + psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + } else { + psAssert (solution->n == numRows, "oops"); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x]; + // XXX EAM : this was set to [x] which is wrong (numCols > numRows) + psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]); + } + corr->data.F64[y][1] -= solution->data.F32[y] ; + } + } + psFree(solution); + } + psFree(yaxisPoly); + psFree(xaxisPoly); + psFree(yaxisIndices); + psFree(yaxisMask); + psFree(yaxisData); + psFree(xaxisData); + + psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, + "Pattern row correction", corr); + psFree(corr); + + psFree(indices); + psFree(xFit); + psFree(yFit); psFree(clip); psFree(clipMask); @@ -1316,2 +1855,3 @@ return true; } + Index: /trunk/psModules/src/detrend/pmPatternIO.c =================================================================== --- /trunk/psModules/src/detrend/pmPatternIO.c (revision 41891) +++ /trunk/psModules/src/detrend/pmPatternIO.c (revision 41892) @@ -8,4 +8,5 @@ #include "pmFPAview.h" #include "pmFPAfile.h" +#include "pmFPAUtils.h" #include "pmFPAfileFitsIO.h" #include "pmFPAHeader.h" @@ -452,2 +453,109 @@ return pmReadoutReadPattern(readout, file->fits); } + +/**************** PatternRowAmp(litude) I/O *************************/ + +bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) + { + // read the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "Pattern Row Amplitude must be read at the FPA level"); + return false; + } + + if (!pmPatternRowAmpReadFPA (file)) { + psError(PS_ERR_IO, false, "Failed to read Pattern Row Amplitude for fpa"); + return false; + } + return true; + } + +// read in all chip-level Pattern Row Amplitude data for this FPA +bool pmPatternRowAmpReadFPA (pmFPAfile *file) { + + if (!pmPatternRowAmpReadChips (file)) { + psError(PS_ERR_IO, false, "Failed to read Pattern Row Amplitude for chips"); + return false; + } + + return true; +} + +// read the set of tables, one for each chip +bool pmPatternRowAmpReadChips (pmFPAfile *file) { + + bool haveData, status; + + // loop over the extensions + // for each extension, use the extname (eg, XY01.ptn) to assign to a chip + + // move to the start of the file + haveData = psFitsMoveExtNum (file->fits, 1, false); + if (!haveData) { + psError(PS_ERR_IO, false, "Failed to read even the first extension?"); + return false; + } + + int nGood = 0; + int nTotal = 0; + while (haveData) { + + // load the header + psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!header) psAbort("cannot read model header"); + + // load the full model in one shot + psArray *model = psFitsReadTable (file->fits); + if (!model) psAbort("cannot read model"); + + // determine the chip: + char *extname = psMetadataLookupStr (&status, header, "EXTNAME"); + psLogMsg ("psModules.detrend", 8, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname); + nTotal += model->n; + + // I expect to find a name of the form: chipName.ptn (eg, XY01.ptn) + // where chipName like 'XY01' + psAssert (strlen(extname) == 8, "invalid extension %s", extname); + psAssert (extname[5] == 'p', "invalid extension %s", extname); + psAssert (extname[6] == 't', "invalid extension %s", extname); + psAssert (extname[7] == 'n', "invalid extension %s", extname); + + char chipName[5]; + strncpy (chipName, extname, 4); + chipName[4] = 0; + + pmChip *chip = pmConceptsChipFromName (file->fpa, chipName); + if (!chip) psAbort ("invalid chip?"); + + // parse the model entries + for (int i = 0; i < model->n; i++) { + psMetadata *row = model->data[i]; + psAssert (row, "missing model row"); + + char *cellName = psMetadataLookupStr(&status, row, "CELL_NAME"); + if (!cellName) continue; + + float amplitude = psMetadataLookupF32(&status, row, "VALUE_MEDIAN"); + + int cellNumber = pmChipFindCell (chip, cellName); + psAssert ((cellNumber >=0) && (cellNumber < chip->cells->n), "invalid cell number"); + + pmCell *cell = chip->cells->data[cellNumber]; + if (!cell) continue; + + psAssert (cell->analysis, "oops"); + + psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude); + nGood ++; + } + + psFree (model); + psFree (header); + + // move to the next extension + haveData = psFitsMoveExtNum (file->fits, 1, true); + } + psLogMsg ("psModules.detrend", 4, "read %d of %d rows from Pattern Row Amplitude file\n", nGood, nTotal); + + return true; +} Index: /trunk/psModules/src/detrend/pmPatternIO.h =================================================================== --- /trunk/psModules/src/detrend/pmPatternIO.h (revision 41891) +++ /trunk/psModules/src/detrend/pmPatternIO.h (revision 41892) @@ -35,3 +35,8 @@ ); + +bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmPatternRowAmpReadFPA (pmFPAfile *file); +bool pmPatternRowAmpReadChips (pmFPAfile *file); + #endif Index: /trunk/psModules/src/imcombine/pmStack.c =================================================================== --- /trunk/psModules/src/imcombine/pmStack.c (revision 41891) +++ /trunk/psModules/src/imcombine/pmStack.c (revision 41892) @@ -756,4 +756,7 @@ // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used // check for set bits and increment counter as appropriate + // count the number of times a given mask bit is set in the input pixels. + // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only + // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs' { psImageMaskType value = 0x0001; @@ -961,5 +964,5 @@ // KMM values; - float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN; + float Punimodal = 1.0, KMMmean = NAN, KMMsigma = NAN, KMMpi = NAN; int KMM_MINIMUM_INPUTS = 6; bool useKMM = false; @@ -1510,4 +1513,130 @@ } +# define MIN_GOOD_PERCENTILE 2 +bool pmStackCombineByPercentile( + pmReadout *combined, + psArray *stackData, + psF64 minRange, + psF64 maxRange, + psImageMaskType badMaskBits, // treat these bits as 'bad' + psImageMaskType suspectMaskBits, // treat these bits as 'suspect' + psImageMaskType blankMaskBits // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) +) { + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + + // we need to copy the readouts to their own array for validation + psArray *stackReadouts = psArrayAlloc(stackData->n); + for (int i = 0; i < stackData->n; i++) { + pmStackData *data = stackData->data[i]; // Stack data for this input + pmReadout *ro = data->readout; // Readout of interest + stackReadouts->data[i] = psMemIncrRefCounter(ro); + } + + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) { + psError(psErrorCodeLast(), false, "Input stack is not valid."); + psFree(stackReadouts); + return false; + } + psFree(stackReadouts); + + psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32); + psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK); + + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { + + int nGood = 0; + for (int i = 0; i < stackData->n; i++) { + + pmStackData *data = stackData->data[i]; // Stack data for this input + pmReadout *ro = data->readout; // Readout of interest + psImage *image = ro->image; + psImage *mask = ro->mask; + + int xIn = x - data->readout->col0; + int yIn = y - data->readout->row0; // Coordinates on input readout + + if (!isfinite(image->data.F32[yIn][xIn])) continue; + if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue; + // XXX need to test the input mask as well here + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits); + + // XXX save the count of suspect bits + +# if (0) + // count the number of times a given mask bit is set in the input pixels. + // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only + // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs' + if (1) { + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) { + // XXX nGoodBits[nbit] ++; + } + value <<= 1; + } + } +# endif + + pixelData->data.F32[nGood] = image->data.F32[yIn][xIn]; + nGood ++; + } + pixelData->n = nGood; + psVectorSortInPlace (pixelData); + + if (nGood < MIN_GOOD_PERCENTILE) { + combined->image->data.F32[y][x] = NAN; + combined->variance->data.F32[y][x] = NAN; + combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 1; + continue; + } + +# if (0) +# define SUSPECT_FRACTION 0.65 + // set the mask bits if nGoodBits[i] > f*numGood + if (0) { + psImageMaskType value = 0x0001; + for (int nbit = 0; nbit < 16; nbit ++) { + if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) { + *goodMask |= value; + } + value <<= 1; + } + } +# endif + + int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5 + int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 + int Npt = Ne - Ns; + + float sum = 0.0; + for (int n = Ns; n < Ne; n++) { + sum += pixelData->data.F32[n]; + } + float mean = sum / (float) Npt; + + float varSum = 0.0; + for (int n = Ns; n < Ne; n++) { + varSum += SQ(pixelData->data.F32[n] - mean); + } + // variance on the mean (stdev / sqrt(N))^2 + float varValue = varSum / (Npt - 1) / Npt; + + // XXX this is probably not the correct variance to save + // the predicted variance should be the 1 / sum (1/var) + + combined->image->data.F32[y][x] = mean; + combined->variance->data.F32[y][x] = varValue; + combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + } + } + + psFree(pixelData); + psFree(pixelMask); + + return (true); +} + /// Stack input images bool pmStackCombine( @@ -1644,10 +1773,4 @@ expweight = expmaps->variance; } - -/* // Pre-reject inputs using KMM bimodality test. */ -/* if (1) { */ -/* /\* KMMRejectUnpopular(input,x,y); *\/ */ -/* rejection = true; */ -/* } */ // Set up rejection list @@ -1660,20 +1783,6 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { + // this is only for TESTING: CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); - -/* // Pre-reject inputs using KMM bimodality test. */ -/* if (1) { */ -/* KMMRejectUnpopular(input,x,y); */ -/* /\* rejection = true; *\/ */ -/* } */ -/* else { */ -/* KMMRejectBright(input,x,y); */ -/* } */ - -#ifdef TESTING - if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { - fprintf (stderr, "combine\n"); - } -# endif psVector *reject = NULL; // Images to reject for this pixel Index: /trunk/psModules/src/imcombine/pmStack.h =================================================================== --- /trunk/psModules/src/imcombine/pmStack.h (revision 41891) +++ /trunk/psModules/src/imcombine/pmStack.h (revision 41892) @@ -47,4 +47,14 @@ ); +bool pmStackCombineByPercentile( + pmReadout *combined, + psArray *input, + psF64 minRange, + psF64 maxRange, + psImageMaskType badMaskBits, // treat these bits as 'bad' + psImageMaskType suspectMaskBits, // treat these bits as 'suspect' + psImageMaskType blankMaskBits // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) + ); + /// Stack input images bool pmStackCombine(pmReadout *combined,///< Combined readout (output) Index: /trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 41891) +++ /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 41892) @@ -57,12 +57,31 @@ // followed by a zero-size matrix, followed by the table data -bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +/** this file is used to generate c-code for a number of different formats. the perl program + ** 'mksource.pl' reads this template file and generates an output file for one of the named + ** formats, e.g., PS1_V1 or PS1_DV4 (see list in mksource.pl). + + ** Any instances of the word @CMFMODE@ are replaced by the named version of the format. + + ** Any line which starts with a command of the form: @COMMAND[,COMMAND...]@ is kept or removed + ** from the output c-code depending on the commands. These may be: ALL (lines are kept for + ** all formats) or the name of a format. A line with just a plain format will be kept only + ** for that format. The format name may be optionally proceded by >, <, <=, >=, ! or =. In + ** these cases, the line is kept if the format is greater / lesser / not the same, etc as the + ** listed format. The sequence of formats is defined for series (V, SV, DV) in mksource.pl so + ** that, e.g., >PS1_V3 would include e.g., PS1_V4 and PS1_V5, but not PS1_DV4. The number in + ** a series may be replaced with ? to stand for the full set. The format + ** commands may be grouped on a line separated by commas. Note that only ! (not) is + ** restrictive: only this rule will reduce the number of matches. So for example, + ** @ALL,= PS1_DV3 will be matched by + ** ALL. But @ALL,!PS1_DV3@ will generate a line for any format except PS1_DV3. + + **/ + +bool pmSourcesWrite_CMF_@CMFMODE@_New (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) { + // fprintf (stderr, "writing with %s\n", __func__); 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; pmChip *chip = readout->parent->parent; @@ -73,12 +92,10 @@ pmSource *source = (pmSource *) sources->data[0]; if (source->seq == -1) { - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); } else { - sources = psArraySort (sources, pmSourceSortBySeq); - } - } - - table = psArrayAllocEmpty (sources->n); + sources = psArraySort (sources, pmSourceSortBySeq); + } + } float magOffset; @@ -88,12 +105,206 @@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + // before we can generate the table structure, we need to know if certain fields were measured. + // this information is not available a priori, so we need to check all sources + + bool haveLensingOBJsmear = false; + bool haveLensingOBJshear = false; + bool haveLensingPSFsmear = false; + bool haveLensingPSFshear = false; + bool haveLensingPSF = false; + + for (int i = 0; (i < sources->n) && !haveLensingOBJsmear && !haveLensingOBJshear && !haveLensingPSFsmear && !haveLensingPSFshear && !haveLensingPSF; i++) { + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + if (source->lensingOBJ && source->lensingOBJ->smear) haveLensingOBJsmear = true; + if (source->lensingOBJ && source->lensingOBJ->shear) haveLensingOBJshear = true; + if (source->lensingPSF && source->lensingPSF->smear) haveLensingPSFsmear = true; + if (source->lensingPSF && source->lensingPSF->shear) haveLensingPSFshear = true; + if (source->lensingPSF ) haveLensingPSF = true; + } + + /************ generate the table columns *****************/ + + // before we allocate the table, generate an empty array of table columns and generate them + psArray *tableColumns = psArrayAllocEmpty (100); + + // add the named / typed columns to the collection of columns + @ALL@ psFitsTableColumnAdd (tableColumns, "IPP_IDET", PS_DATA_U32); // "IPP detection identifier index" + @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF", PS_DATA_F32); // "PSF x coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF", PS_DATA_F32); // "PSF y coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "X_PSF_SIG", PS_DATA_F32); // "Sigma in PSF x coordinate" + @ALL@ psFitsTableColumnAdd (tableColumns, "Y_PSF_SIG", PS_DATA_F32); // "Sigma in PSF y coordinate" + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @PS1_V1@ psFitsTableColumnAdd (tableColumns, "RA_PSF", PS_DATA_F32); // "PSF RA coordinate (degrees)" + @PS1_V1@ psFitsTableColumnAdd (tableColumns, "DEC_PSF", PS_DATA_F32); // "PSF DEC coordinate (degrees)" + + @ALL@ psFitsTableColumnAdd (tableColumns, "POSANGLE", PS_DATA_F32); // "position angle at source (degrees)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PLTSCALE", PS_DATA_F32); // "plate scale at source (arcsec/pixel)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG", PS_DATA_F32); // "PSF fit instrumental magnitude" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG_SIG", PS_DATA_F32); // "Sigma of PSF instrumental magnitude" + + @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX", PS_DATA_F32); // "PSF fit instrumental flux (counts)" + @ALL,!PS1_V1,!PS1_V2@ psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX_SIG", PS_DATA_F32); // "Sigma of PSF instrumental flux" + + @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG", PS_DATA_F32); // "magnitude in standard aperture" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RAW", PS_DATA_F32); // "magnitude in reported aperture" + @ALL@ psFitsTableColumnAdd (tableColumns, "AP_MAG_RADIUS", PS_DATA_F32); // "radius used for aperture mags" + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX", PS_DATA_F32); // "instrumental flux in standard aperture" + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX_SIG", PS_DATA_F32); // "aperture flux error" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "AP_NPIX", PS_DATA_S32); // "aperture unmasked pixels" + + @=PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude" + @ALL@ psFitsTableColumnAdd (tableColumns, "SKY", PS_DATA_F32); // "Sky level" + @ALL@ psFitsTableColumnAdd (tableColumns, "SKY_SIGMA", PS_DATA_F32); // "Sigma of sky level" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_CHISQ", PS_DATA_F32); // "Chisq of PSF-fit" + @ALL@ psFitsTableColumnAdd (tableColumns, "CR_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to CF" + @ALL@ psFitsTableColumnAdd (tableColumns, "EXT_NSIGMA", PS_DATA_F32); // "Nsigma deviations from PSF to EXT" + + // PSF shape parameters: + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MAJOR", PS_DATA_F32); // "PSF width (major axis)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_MINOR", PS_DATA_F32); // "PSF width (minor axis)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_THETA", PS_DATA_F32); // "PSF orientation angle" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_CORE", PS_DATA_F32); // "k term if defined" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MAJ", PS_DATA_F32); // "PSF FWHM (major axis)" + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MIN", PS_DATA_F32); // "PSF FWHM (minor axis)" + + // psf data quality + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_QF", PS_DATA_F32); // "PSF coverage/quality factor (bad)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "PSF_QF_PERFECT", PS_DATA_F32); // "PSF coverage/quality factor (poor)" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NDOF", PS_DATA_S32); // "degrees of freedom" + @ALL@ psFitsTableColumnAdd (tableColumns, "PSF_NPIX", PS_DATA_S32); // "number of pixels in fit" + + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XX", PS_DATA_F32); // "second moments (X^2)" + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_XY", PS_DATA_F32); // "second moments (X*Y)" + @ALL@ psFitsTableColumnAdd (tableColumns, "MOMENTS_YY", PS_DATA_F32); // "second moments (Y*Y)" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3C", PS_DATA_F32); // "third momemt cos theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M3S", PS_DATA_F32); // "third momemt sin theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4C", PS_DATA_F32); // "fourth momemt cos theta" + @>PS1_V2,PS1_SV?@ psFitsTableColumnAdd (tableColumns, "MOMENTS_M4S", PS_DATA_F32); // "fourth momemt sin theta" + + // Lensing parameters are only written if they are measured + if (haveLensingOBJsmear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_OBJ", PS_DATA_F32); // "smear polarizability element (objects)" + } + if (haveLensingOBJshear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + } + if (false) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_OBJ", PS_DATA_F32); // "shear polarizability element (objects)" + } + + if (haveLensingPSFsmear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_PSF", PS_DATA_F32); // "smear polarizability element (PSFs)" + } + if (haveLensingPSFshear) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + } + if (haveLensingPSF) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E1_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "LENS_E2_PSF", PS_DATA_F32); // "shear polarizability element (PSFs)" + } + + // if lensing params exist also include the backmapped chipID and chip coordinates + if (haveLensingPSF) { + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" + @>PS1_V4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_R1", PS_DATA_F32); // "first radial moment" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "MOMENTS_RH", PS_DATA_F32); // "half radial moment" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX", PS_DATA_F32); // "Kron Flux (in 2.5 R1)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_ERR", PS_DATA_F32); // "Kron Flux Error" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_INNER", PS_DATA_F32); // "Kron Flux (in 1.0 R1)" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "KRON_FLUX_OUTER", PS_DATA_F32); // "Kron Flux (in 4.0 R1)" + + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_RAD", PS_DATA_F32); // "Radius where object hits sky" + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_FLUX", PS_DATA_F32); // "Flux / pix where object hits sky" + @>PS1_V3@ psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_SLOPE", PS_DATA_F32); // "d(Flux/pix)/dRadius where object hits sky" + + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM", PS_DATA_S16); // "id of warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X", PS_DATA_S16); // "x coord in warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y", PS_DATA_S16); // "y coord in warp input chip" + @>PS1_DV4@ psFitsTableColumnAdd (tableColumns, "PADDING3", PS_DATA_S16); // "more padding" + + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NPOS", PS_DATA_S32); // "nPos (n pix > 3 sigma)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_FRATIO", PS_DATA_F32); // "fPos / (fPos + fNeg)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_BAD", PS_DATA_F32); // "nPos / (nPos + nNeg)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_MASK", PS_DATA_F32); // "nPos / (nPos + nMask)" + @PS1_DV?@ psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_ALL", PS_DATA_F32); // "nPos / (nGood + nMask + nBad)" + + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_P", PS_DATA_F32); // "distance to positive match source" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_P", PS_DATA_F32); // "signal-to-noise of pos match src" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_R_M", PS_DATA_F32); // "distance to negative match source" + @>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "DIFF_SN_M", PS_DATA_F32); // "signal-to-noise of neg match src" + + @ALL@ psFitsTableColumnAdd (tableColumns, "FLAGS", PS_DATA_U32); // "psphot analysis flags" + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableColumnAdd (tableColumns, "FLAGS2", PS_DATA_U32); // "psphot analysis flags" + @PS1_V3,PS1_SV2@ psFitsTableColumnAdd (tableColumns, "PADDING2", PS_DATA_S32); // "more padding" + @ALL@ psFitsTableColumnAdd (tableColumns, "N_FRAMES", PS_DATA_U16); // "Number of frames overlapping source center" + + @ALL@ psFitsTableColumnAdd (tableColumns, "PADDING", PS_DATA_S16); // "more padding" + + // note that this definition is inconsistent with the definition in + // Ohana/src/libautocode/def. This version creates a table with data not + // properly aligned with the 8-byte boundaries. The structure defined by + // libautocode does this, but has a different order of elements (and adds + // padding2 to fix things). We have generated may files with PS1_SV1 as is, so + // I'll leave it. But in future a PS1_SV2 should be forced to match + // libautocode. Note that addstar knows to detect the alternate version of + // PS1_SV1 and correctly interpret its fields. + + // 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; + + // generate an FITS table using the array of columns defined above, with space for all rows + psFitsTable *table = psFitsTableCreate (tableColumns, sources->n); + + /************ write the data to the table *****************/ + // 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 (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // this is the source associated with this image pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // If source->seq is -1, source was generated in this analysis. If source->seq is @@ -102,133 +313,361 @@ // generated on Alloc, and would thus be wrong for read in sources. if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - @PS1_DV?@ pmSourceDiffStats diffStats; - @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); - @PS1_DV?@ if (source->diffStats) { - @PS1_DV?@ diffStats = *source->diffStats; - @PS1_DV?@ } - - row = psMetadataAlloc (); + source->seq = i; + } + + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); + + @PS1_DV?@ pmSourceDiffStats diffStats; + @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); + @PS1_DV?@ if (source->diffStats) { diffStats = *source->diffStats;} // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); - - // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision - @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); - @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - - @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", source->apRadius); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); - - @=PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + @ALL@ psFitsTableSetU32 (table, i, "IPP_IDET", source->seq); // "IPP detection identifier index", + @ALL@ psFitsTableSetF32 (table, i, "X_PSF", outputs.xPos); // "PSF x coordinate", + @ALL@ psFitsTableSetF32 (table, i, "Y_PSF", outputs.yPos); // "PSF y coordinate", + @ALL@ psFitsTableSetF32 (table, i, "X_PSF_SIG", outputs.xErr); // "Sigma in PSF x coordinate", + @ALL@ psFitsTableSetF32 (table, i, "Y_PSF_SIG", outputs.yErr); // "Sigma in PSF y coordinate", + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @PS1_V1@ psFitsTableSetF32 (table, i, "RA_PSF", outputs.ra); // "PSF RA coordinate (degrees)", + @PS1_V1@ psFitsTableSetF32 (table, i, "DEC_PSF", outputs.dec); // "PSF DEC coordinate (degrees)", + + @ALL@ psFitsTableSetF32 (table, i, "POSANGLE", outputs.posAngle); // "position angle at source (degrees)", + @ALL@ psFitsTableSetF32 (table, i, "PLTSCALE", outputs.pltScale); // "plate scale at source (arcsec/pixel)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG", source->psfMag); // "PSF fit instrumental magnitude", + @ALL@ psFitsTableSetF32 (table, i, "PSF_INST_MAG_SIG", source->psfMagErr); // "Sigma of PSF instrumental magnitude", + + @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX", source->psfFlux); // "PSF fit instrumental flux (counts)", + @ALL,!PS1_V1,!PS1_V2@ psFitsTableSetF32 (table, i, "PSF_INST_FLUX_SIG", source->psfFluxErr); // "Sigma of PSF instrumental flux", + + @ALL@ psFitsTableSetF32 (table, i, "AP_MAG", source->apMag); // "magnitude in standard aperture", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "AP_MAG_RAW", source->apMagRaw); // "magnitude in reported aperture", + @ALL@ psFitsTableSetF32 (table, i, "AP_MAG_RADIUS", source->apRadius); // "radius used for aperture mags", + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX", source->apFlux); // "instrumental flux in standard aperture", + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX_SIG", source->apFluxErr); // "aperture flux error", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetS32 (table, i, "AP_NPIX", source->apNpixels); // "aperture unmasked pixels", + + @=PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG", outputs.peakMag); // "Peak flux expressed as magnitude", + @ALL@ psFitsTableSetF32 (table, i, "SKY", source->sky); // "Sky level", + @ALL@ psFitsTableSetF32 (table, i, "SKY_SIGMA", source->skyErr); // "Sigma of sky level", + + @ALL@ psFitsTableSetF32 (table, i, "PSF_CHISQ", outputs.chisq); // "Chisq of PSF-fit", + @ALL@ psFitsTableSetF32 (table, i, "CR_NSIGMA", source->crNsigma); // "Nsigma deviations from PSF to CF", + @ALL@ psFitsTableSetF32 (table, i, "EXT_NSIGMA", source->extNsigma); // "Nsigma deviations from PSF to EXT", // PSF shape parameters: - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); + @ALL@ psFitsTableSetF32 (table, i, "PSF_MAJOR", outputs.psfMajor); // "PSF width (major axis)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_MINOR", outputs.psfMinor); // "PSF width (minor axis)", + @ALL@ psFitsTableSetF32 (table, i, "PSF_THETA", outputs.psfTheta); // "PSF orientation angle", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_CORE", outputs.psfCore); // "k term if defined", + + // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MAJ", outputs.psfMajorFWHM); // "PSF FWHM (major axis)", + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MIN", outputs.psfMinorFWHM); // "PSF FWHM (minor axis)", + + // psf data quality + @ALL@ psFitsTableSetF32 (table, i, "PSF_QF", source->pixWeightNotBad); // "PSF coverage/quality factor (bad)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "PSF_QF_PERFECT", source->pixWeightNotPoor); // "PSF coverage/quality factor (poor)", + @ALL@ psFitsTableSetS32 (table, i, "PSF_NDOF", outputs.nDOF); // "degrees of freedom", + @ALL@ psFitsTableSetS32 (table, i, "PSF_NPIX", outputs.nPix); // "number of pixels in fit", + + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XX", moments.Mxx); // "second moments (X^2)", + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_XY", moments.Mxy); // "second moments (X*Y)", + @ALL@ psFitsTableSetF32 (table, i, "MOMENTS_YY", moments.Myy); // "second moments (Y*Y)", + + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3C", moments.M_c3); // "third momemt cos theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M3S", moments.M_s3); // "third momemt sin theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4C", moments.M_c4); // "fourth momemt cos theta", + @>PS1_V2,PS1_SV?@ psFitsTableSetF32 (table, i, "MOMENTS_M4S", moments.M_s4); // "fourth momemt sin theta", + + // Lensing parameters: + if (source->lensingOBJ && source->lensingOBJ->smear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_OBJ", source->lensingOBJ->smear->X11); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_OBJ", source->lensingOBJ->smear->X12); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_OBJ", source->lensingOBJ->smear->X22); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_OBJ", source->lensingOBJ->smear->e1); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_OBJ", source->lensingOBJ->smear->e2); // "smear polarizability element (objects)", + } + if (source->lensingOBJ && source->lensingOBJ->shear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_OBJ", source->lensingOBJ->shear->X11); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_OBJ", source->lensingOBJ->shear->X12); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_OBJ", source->lensingOBJ->shear->X22); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_OBJ", source->lensingOBJ->shear->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_OBJ", source->lensingOBJ->shear->e2); // "shear polarizability element (objects)", + } + if (false && source->lensingOBJ) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingOBJ->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingOBJ->e2); // "shear polarizability element (objects)", + } + + if (source->lensingPSF && source->lensingPSF->smear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SM_PSF", source->lensingPSF->smear->X11); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SM_PSF", source->lensingPSF->smear->X12); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SM_PSF", source->lensingPSF->smear->X22); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SM_PSF", source->lensingPSF->smear->e1); // "smear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SM_PSF", source->lensingPSF->smear->e2); // "smear polarizability element (objects)", + } + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X11_SH_PSF", source->lensingPSF->shear->X11); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X12_SH_PSF", source->lensingPSF->shear->X12); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_X22_SH_PSF", source->lensingPSF->shear->X22); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_SH_PSF", source->lensingPSF->shear->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_SH_PSF", source->lensingPSF->shear->e2); // "shear polarizability element (objects)", + } + if (source->lensingPSF) { + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E1_PSF", source->lensingPSF->e1); // "shear polarizability element (objects)", + @>PS1_V4@ psFitsTableSetF32 (table, i, "LENS_E2_PSF", source->lensingPSF->e2); // "shear polarizability element (objects)", + } + + // if lensing params exist also include the backmapped chipID and chip coordinates + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", + @>PS1_V4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_R1", moments.Mrf); // "first radial moment", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "MOMENTS_RH", moments.Mrh); // "half radial moment", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX", moments.Krf); // "Kron Flux (in 2.5 R1)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_ERR", moments.dKrf); // "Kron Flux Error", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_INNER", moments.Kinner); // "Kron Flux (in 1.0 R1)", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetF32 (table, i, "KRON_FLUX_OUTER", moments.Kouter); // "Kron Flux (in 4.0 R1)", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_RAD", source->skyRadius); // "Radius where object hits sky", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_FLUX", source->skyFlux); // "Flux / pix where object hits sky", + @>PS1_V3@ psFitsTableSetF32 (table, i, "SKY_LIMIT_SLOPE", source->skySlope); // "d(Flux/pix)/dRadius where object hits sky", + + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_NUM", source->chipNum); // "id of warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_X", source->chipX); // "x coord in warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "SRC_CHIP_Y", source->chipY); // "y coord in warp input chip", + @>PS1_DV4@ psFitsTableSetS16 (table, i, "PADDING3", 0); // "more padding", + + @PS1_DV?@ psFitsTableSetS32 (table, i, "DIFF_NPOS", diffStats.nGood); // "nPos (n pix > 3 sigma)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_FRATIO", diffStats.fRatio); // "fPos / (fPos + fNeg)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_BAD", diffStats.nRatioBad); // "nPos / (nPos + nNeg)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_MASK", diffStats.nRatioMask); // "nPos / (nPos + nMask)", + @PS1_DV?@ psFitsTableSetF32 (table, i, "DIFF_NRATIO_ALL", diffStats.nRatioAll); // "nPos / (nGood + nMask + nBad)", + + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_P", diffStats.Rp); // "distance to positive match source", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_P", diffStats.SNp); // "signal-to-noise of pos match src", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_R_M", diffStats.Rm); // "distance to negative match source", + @>PS1_DV1@ psFitsTableSetF32 (table, i, "DIFF_SN_M", diffStats.SNm); // "signal-to-noise of neg match src", + + @ALL@ psFitsTableSetU32 (table, i, "FLAGS", source->mode); // "psphot analysis flags", + @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetU32 (table, i, "FLAGS2", source->mode2); // "psphot analysis flags", + @PS1_V3,PS1_SV2@ psFitsTableSetS32 (table, i, "PADDING2", 0); // "padding", + + @ALL@ psFitsTableSetU16 (table, i, "N_FRAMES", source->nFrames); // "Number of frames overlapping source center", + @ALL@ psFitsTableSetS16 (table, i, "PADDING", 0); // "padding", + } + + // XXX why do we make a copy here to be supplemented with the masks? why not do this in the calling function? + psMetadata *header = psMetadataCopy(NULL, tableHeader); + pmSourceMasksHeader(header); + + psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); + if (!psFitsWriteTableNew(fits, header, table, extname)) { + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree(table); + psFree(header); + return false; + } + psFree(table); + psFree(header); + + return true; +} + +bool pmSourcesWrite_CMF_@CMFMODE@_Old (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + // fprintf (stderr, "writing with %s\n", __func__); + 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; + + pmChip *chip = readout->parent->parent; + + // 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, pmSourceSortByFlux); + } else { + sources = psArraySort (sources, pmSourceSortBySeq); + } + } + + table = psArrayAllocEmpty (sources->n); + + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + + // 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 (int i = 0; i < sources->n; i++) { + // this is the source associated with this image + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // 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; + } + + // set the 'best' values for various output fields: + pmSourceOutputs outputs; + pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); + + pmSourceOutputsMoments moments; + pmSourceOutputsSetMoments (&moments, source); + + @PS1_DV?@ pmSourceDiffStats diffStats; + @PS1_DV?@ pmSourceDiffStatsInit(&diffStats); + @PS1_DV?@ if (source->diffStats) { + @PS1_DV?@ diffStats = *source->diffStats; + @PS1_DV?@ } + + row = psMetadataAlloc (); + + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment" + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); + + // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision + @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F32, "PSF RA coordinate (degrees)", outputs.ra); + @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F32, "PSF DEC coordinate (degrees)", outputs.dec); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); + + @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); + @ALL,!PS1_V1,!PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); + @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", source->apRadius); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX", PS_DATA_S32, "aperture unmasked pixels", source->apNpixels); + + @=PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); + + // PSF shape parameters: + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE", PS_DATA_F32, "k term if defined", outputs.psfCore); - // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values + // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ", PS_DATA_F32, "PSF FWHM (major axis)", outputs.psfMajorFWHM); @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN", PS_DATA_F32, "PSF FWHM (minor axis)", outputs.psfMinorFWHM); // psf data quality - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - // Lensing parameters: - if (source->lensingOBJ && source->lensingOBJ->smear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e2); - } - if (source->lensingOBJ && source->lensingOBJ->shear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e2); - } - if (false && source->lensingOBJ) { - // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e2); - } - - if (source->lensingPSF && source->lensingPSF->smear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e2); - } - if (source->lensingPSF && source->lensingPSF->shear) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X22); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); - } - if (source->lensingPSF) { - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e1); - @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e2); - } + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); + @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); + + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + @>PS1_V2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + // Lensing parameters: + if (source->lensingOBJ && source->lensingOBJ->smear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ", PS_DATA_F32, "smear polarizability element (objects)", source->lensingOBJ->smear->e2); + } + if (source->lensingOBJ && source->lensingOBJ->shear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->shear->e2); + } + if (false && source->lensingOBJ) { + // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingOBJ->e2); + } + + if (source->lensingPSF && source->lensingPSF->smear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF", PS_DATA_F32, "smear polarizability element (objects)", source->lensingPSF->smear->e2); + } + if (source->lensingPSF && source->lensingPSF->shear) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X11); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X12); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->X22); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->shear->e2); + } + if (source->lensingPSF) { + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e1); + @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF", PS_DATA_F32, "shear polarizability element (objects)", source->lensingPSF->e2); + } // if lensing params exist also include the backmapped chipID and chip coordinates - if (source->lensingPSF && source->lensingPSF->shear) { + if (source->lensingPSF && source->lensingPSF->shear) { @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); @>PS1_V4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); @@ -244,40 +683,40 @@ @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD", PS_DATA_F32, "Radius where object hits sky", source->skyRadius); - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX", PS_DATA_F32, "Flux / pix where object hits sky", source->skyFlux); - @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE", PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky", source->skySlope); - - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); - @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); - - @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); - - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); - @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - @PS1_V3,PS1_SV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); - @PS1_SV?@ - - // note that this definition is inconsistent with the definition in - // Ohana/src/libautocode/def. This version creates a table with data not - // properly aligned with the 8-byte boundaries. The structure defined by - // libautocode does this, but has a different order of elements (and adds - // padding2 to fix things). We have generated may files with PS1_SV1 as is, so - // I'll leave it. But in future a PS1_SV2 should be forced to match - // libautocode. Note that addstar knows to detect the alternate version of - // PS1_SV1 and correctly interpret its fields. - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD", PS_DATA_F32, "Radius where object hits sky", source->skyRadius); + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX", PS_DATA_F32, "Flux / pix where object hits sky", source->skyFlux); + @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE", PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky", source->skySlope); + + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM", PS_DATA_S16, "id of warp input chip", source->chipNum); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X", PS_DATA_S16, "x coord in warp input chip", source->chipX); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y", PS_DATA_S16, "y coord in warp input chip", source->chipY); + @>PS1_DV4@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING3", PS_DATA_S16, "more padding", 0); + + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); + @PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); + + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); + @>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); + @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + @PS1_V3,PS1_SV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); + @PS1_SV?@ + + // note that this definition is inconsistent with the definition in + // Ohana/src/libautocode/def. This version creates a table with data not + // properly aligned with the 8-byte boundaries. The structure defined by + // libautocode does this, but has a different order of elements (and adds + // padding2 to fix things). We have generated may files with PS1_SV1 as is, so + // I'll leave it. But in future a PS1_SV2 should be forced to match + // libautocode. Note that addstar knows to detect the alternate version of + // PS1_SV1 and correctly interpret its fields. + + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); + @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); psArrayAdd (table, 100, row); @@ -320,7 +759,16 @@ } +bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) +{ + // bool status = pmSourcesWrite_CMF_@CMFMODE@_Old (fits, readout, sources, imageHeader, tableHeader, extname, recipe); + bool status = pmSourcesWrite_CMF_@CMFMODE@_New (fits, readout, sources, imageHeader, tableHeader, extname, recipe); + return status; +} + // read in a readout from the fits file -psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) +psArray *pmSourcesRead_CMF_@CMFMODE@_New (psFits *fits, psMetadata *header) { + // fprintf (stderr, "reading with %s\n", __func__); + PS_ASSERT_PTR_NON_NULL(fits, false); PS_ASSERT_PTR_NON_NULL(header, false); @@ -355,5 +803,216 @@ psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - // convert the table to the pmSource entriesa + // reads full table into memory + psFitsTable *table = psFitsReadTableNew (fits); + + // convert the table to the pmSource entries + for (int i = 0; i < numSources; i++) { + 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; + + @ALL@ source->seq = psFitsTableGetU32 (&status, table, i, "IPP_IDET"); + @ALL@ PAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF"); + @ALL@ PAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF"); + @ALL@ dPAR[PM_PAR_XPOS] = psFitsTableGetF32 (&status, table, i, "X_PSF_SIG"); + @ALL@ dPAR[PM_PAR_YPOS] = psFitsTableGetF32 (&status, table, i, "Y_PSF_SIG"); + @ALL@ axes.major = psFitsTableGetF32 (&status, table, i, "PSF_MAJOR"); + @ALL@ axes.minor = psFitsTableGetF32 (&status, table, i, "PSF_MINOR"); + @ALL@ axes.theta = psFitsTableGetF32 (&status, table, i, "PSF_THETA"); + @ALL@ axes.theta = axes.theta * PS_RAD_DEG; + + @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { + @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psFitsTableGetF32 (&status, table, i, "PSF_CORE"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ } + + @ALL@ PAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY"); + @ALL@ dPAR[PM_PAR_SKY] = psFitsTableGetF32 (&status, table, i, "SKY_SIGMA"); + @ALL@ source->sky = PAR[PM_PAR_SKY]; + @ALL@ source->skyErr = dPAR[PM_PAR_SKY]; + + // XXX use these to determine PAR[PM_PAR_I0]? + @ALL@ source->psfMag = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG"); + @ALL@ source->psfMagErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG_SIG"); + @ALL@ source->apMag = psFitsTableGetF32 (&status, table, i, "AP_MAG"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psFitsTableGetF32 (&status, table, i, "AP_MAG_RAW"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psFitsTableGetF32 (&status, table, i, "AP_FLUX"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psFitsTableGetF32 (&status, table, i, "AP_FLUX_SIG"); + + // XXX use these to determine PAR[PM_PAR_I0] if they exist? + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX_SIG"); + + // XXX this scaling is incorrect: does not include the 2 \pi AREA factor + @ALL@ PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; + @ALL@ dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; + + pmPSF_AxesToModel (PAR, axes, model->class->useReff); + + @ALL@ float peakMag = psFitsTableGetF32 (&status, table, i, "PEAK_FLUX_AS_MAG"); + @ALL@ float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; + + // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) + @ALL@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); + @ALL@ source->peak->rawFlux = peakFlux; + @ALL@ source->peak->smoothFlux = peakFlux; + @ALL@ source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position + @ALL@ source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position + @ALL@ source->peak->dx = dPAR[PM_PAR_XPOS]; + @ALL@ source->peak->dy = dPAR[PM_PAR_YPOS]; + + @ALL@ source->pixWeightNotBad = psFitsTableGetF32 (&status, table, i, "PSF_QF"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psFitsTableGetF32 (&status, table, i, "PSF_QF_PERFECT"); + @ALL@ source->crNsigma = psFitsTableGetF32 (&status, table, i, "CR_NSIGMA"); + @ALL@ source->extNsigma = psFitsTableGetF32 (&status, table, i, "EXT_NSIGMA"); + @ALL@ source->apRadius = psFitsTableGetF32 (&status, table, i, "AP_MAG_RADIUS"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psFitsTableGetS32 (&status, table, i, "AP_NPIX"); + + // note that some older versions used PSF_PROBABILITY: this was not well defined. + @ALL@ model->chisq = psFitsTableGetF32 (&status, table, i, "PSF_CHISQ"); + @ALL@ model->nDOF = psFitsTableGetS32 (&status, table, i, "PSF_NDOF"); + @ALL@ model->nPix = psFitsTableGetS32 (&status, table, i, "PSF_NPIX"); + + @ALL@ source->moments = pmMomentsAlloc (); + @ALL@ source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf + @ALL@ source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf + + @ALL@ source->moments->Mxx = psFitsTableGetF32 (&status, table, i, "MOMENTS_XX"); + @ALL@ source->moments->Mxy = psFitsTableGetF32 (&status, table, i, "MOMENTS_XY"); + @ALL@ source->moments->Myy = psFitsTableGetF32 (&status, table, i, "MOMENTS_YY"); + + // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, + // we are storing enough information so the output will be consistent with the input + @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3C"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxyy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Myyy = -1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3S"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxxx = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4C"); + @>PS1_V2,PS1_SV?@ source->moments->Mxxxy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxxyy = 0.00; + @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4S"); + @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.00; + + // Lensing parameters (on read if PS1_V5+) + if (haveLensOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + source->lensingOBJ->smear = pmLensingParsAlloc(); + source->lensingOBJ->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingOBJ->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_OBJ"); + } + + @>PS1_V4@ source->chipNum = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_NUM"); + @>PS1_V4@ source->chipX = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_X"); + @>PS1_V4@ source->chipY = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_Y"); + + if (haveLensPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + source->lensingPSF->smear = pmLensingParsAlloc(); + source->lensingPSF->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingPSF->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_PSF"); + @>PS1_V4@ source->lensingPSF->e1 = psFitsTableGetF32 (&status, table, i, "LENS_E1_PSF"); + @>PS1_V4@ source->lensingPSF->e2 = psFitsTableGetF32 (&status, table, i, "LENS_E2_PSF"); + } + + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psFitsTableGetF32 (&status, table, i, "MOMENTS_R1"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrh = psFitsTableGetF32 (&status, table, i, "MOMENTS_RH"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFlux = psFitsTableGetF32 (&status, table, i, "KRON_FLUX"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFluxErr = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_ERR"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFinner = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_INNER"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFouter = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_OUTER"); + + @>PS1_V3@ source->skyRadius = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_RAD"); + @>PS1_V3@ source->skyFlux = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_FLUX"); + @>PS1_V3@ source->skySlope = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_SLOPE"); + + @PS1_DV?@ int nPos = psFitsTableGetS32 (&status, table, i, "DIFF_NPOS"); + @PS1_DV?@ if (nPos) { + @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); + @PS1_DV?@ source->diffStats->nGood = nPos; + @PS1_DV?@ source->diffStats->fRatio = psFitsTableGetF32 (&status, table, i, "DIFF_FRATIO"); + @PS1_DV?@ source->diffStats->nRatioBad = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_BAD"); + @PS1_DV?@ source->diffStats->nRatioMask = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_MASK"); + @PS1_DV?@ source->diffStats->nRatioAll = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_ALL"); + + @>PS1_DV1@ source->diffStats->Rp = psFitsTableGetF32 (&status, table, i, "DIFF_R_P"); + @>PS1_DV1@ source->diffStats->SNp = psFitsTableGetF32 (&status, table, i, "DIFF_SN_P"); + @>PS1_DV1@ source->diffStats->Rm = psFitsTableGetF32 (&status, table, i, "DIFF_R_M"); + @>PS1_DV1@ source->diffStats->SNm = psFitsTableGetF32 (&status, table, i, "DIFF_SN_M"); + @PS1_DV?@ } + + @ALL@ source->mode = psFitsTableGetU32 (&status, table, i, "FLAGS"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->mode2 = psFitsTableGetU32 (&status, table, i, "FLAGS2"); + @ALL@ source->nFrames = psFitsTableGetU16 (&status, table, i, "N_FRAMES"); + assert (status); + + sources->data[i] = source; + } + psFree (table); + return sources; +} + +// read in a readout from the fits file +psArray *pmSourcesRead_CMF_@CMFMODE@_Old (psFits *fits, psMetadata *header) +{ + // fprintf (stderr, "reading with %s\n", __func__); + + 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 + int defaultModelType = pmModelClassGetType ("PS_MODEL_GAUSS"); + int modelType = -1; + + // if header does not define the model, default to a gaussian + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL"); + if (PSF_NAME != NULL) { + modelType = pmModelClassGetType (PSF_NAME); + } + // work around bug in psphotFullForce + if (modelType < 0) { + modelType = defaultModelType; + } + // assert (modelType > -1); + + // do we expect to find lensing parameters? + bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ"); + bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF"); + + // 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 entries for (int i = 0; i < numSources; i++) { psMetadata *row = psFitsReadTableRow(fits, i); // Table row @@ -382,8 +1041,8 @@ @ALL@ axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); @ALL@ axes.theta = axes.theta * PS_RAD_DEG; - - @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { - @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ } + + @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) { + @>PS1_V4,>PS1_SV2,>PS1_DV3@ PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ } @ALL@ PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); @@ -396,12 +1055,12 @@ @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); - @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX"); + @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG"); // XXX use these to determine PAR[PM_PAR_I0] if they exist? - // XXX add these to PS1_SV1? - @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -424,9 +1083,9 @@ @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); + @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); @ALL@ source->apRadius = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS"); - @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); + @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -443,6 +1102,6 @@ @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, - // we are storing enough information so the output will be consistent with the input + // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data, + // we are storing enough information so the output will be consistent with the input @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.0; @@ -456,44 +1115,44 @@ @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0; - // Lensing parameters (on read if PS1_V5+) - if (haveLensOBJ) { - source->lensingOBJ = pmSourceLensingAlloc (); - source->lensingOBJ->smear = pmLensingParsAlloc(); - source->lensingOBJ->shear = pmLensingParsAlloc(); - - @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); - @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); - } - - @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM"); - @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X"); - @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y"); - - if (haveLensPSF) { - source->lensingPSF = pmSourceLensingAlloc (); - source->lensingPSF->smear = pmLensingParsAlloc(); - source->lensingPSF->shear = pmLensingParsAlloc(); - - @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); - @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); - @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); - @>PS1_V4@ source->lensingPSF->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_PSF"); - @>PS1_V4@ source->lensingPSF->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_PSF"); - } + // Lensing parameters (on read if PS1_V5+) + if (haveLensOBJ) { + source->lensingOBJ = pmSourceLensingAlloc (); + source->lensingOBJ->smear = pmLensingParsAlloc(); + source->lensingOBJ->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ"); + @>PS1_V4@ source->lensingOBJ->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ"); + } + + @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM"); + @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X"); + @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y"); + + if (haveLensPSF) { + source->lensingPSF = pmSourceLensingAlloc (); + source->lensingPSF->smear = pmLensingParsAlloc(); + source->lensingPSF->shear = pmLensingParsAlloc(); + + @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF"); + @>PS1_V4@ source->lensingPSF->smear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF"); + @>PS1_V4@ source->lensingPSF->shear->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF"); + @>PS1_V4@ source->lensingPSF->e1 = psMetadataLookupF32 (&status, row, "LENS_E1_PSF"); + @>PS1_V4@ source->lensingPSF->e2 = psMetadataLookupF32 (&status, row, "LENS_E2_PSF"); + } @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); @@ -508,18 +1167,18 @@ @>PS1_V3@ source->skySlope = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE"); - @PS1_DV?@ int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); - @PS1_DV?@ if (nPos) { - @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); - @PS1_DV?@ source->diffStats->nGood = nPos; - @PS1_DV?@ source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); - @PS1_DV?@ source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); - @PS1_DV?@ source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); - @PS1_DV?@ source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); - - @>PS1_DV1@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); - @>PS1_DV1@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); - @>PS1_DV1@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); - @>PS1_DV1@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); - @PS1_DV?@ } + @PS1_DV?@ int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); + @PS1_DV?@ if (nPos) { + @PS1_DV?@ source->diffStats = pmSourceDiffStatsAlloc(); + @PS1_DV?@ source->diffStats->nGood = nPos; + @PS1_DV?@ source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); + @PS1_DV?@ source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); + @PS1_DV?@ source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); + @PS1_DV?@ source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); + + @>PS1_DV1@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); + @>PS1_DV1@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); + @>PS1_DV1@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); + @>PS1_DV1@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); + @PS1_DV?@ } @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); @@ -533,4 +1192,10 @@ return sources; +} + +psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) { + // psArray *array = pmSourcesRead_CMF_@CMFMODE@_Old (fits, header); + psArray *array = pmSourcesRead_CMF_@CMFMODE@_New (fits, header); + return array; } @@ -599,17 +1264,17 @@ // write the radial profile apertures to header for (int i = 0; i < radMax->n; i++) { - sprintf (keyword1, "RMIN_%02d", i); - sprintf (keyword2, "RMAX_%02d", i); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); - psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]); + sprintf (keyword1, "RMIN_%02d", i); + sprintf (keyword2, "RMAX_%02d", i); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]); + psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]); } // we write out all sources, regardless of quality. the source flags tell us the state for (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // this is the source associated with this image pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without measurements @@ -639,11 +1304,11 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - float AxialRatio = NAN; - float AxialTheta = NAN; - pmSourceExtendedPars *extpars = source->extpars; - if (extpars) { - AxialRatio = extpars->axes.minor / extpars->axes.major; - AxialTheta = extpars->axes.theta; - } + float AxialRatio = NAN; + float AxialTheta = NAN; + pmSourceExtendedPars *extpars = source->extpars; + if (extpars) { + AxialRatio = extpars->axes.minor / extpars->axes.major; + AxialTheta = extpars->axes.theta; + } psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO", PS_DATA_F32, "Axial Ratio of radial profile", AxialRatio); psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA", PS_DATA_F32, "Angle of radial profile ellipse", AxialTheta); @@ -651,22 +1316,22 @@ // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured + // XXX check flags to see if Pet was measured if (doPetrosian) { - pmSourceExtendedPars *extpars = source->extpars; + pmSourceExtendedPars *extpars = source->extpars; if (extpars) { - // XXX note that this mag is either calibrated or instrumental depending on existence of zero point - float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point - // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - if (repairMagErrors) { - // I need to add the kron error in quadrature becasue pet_error ignores the object flux - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - if (isfinite (Krf) && isfinite (dKrf)) { - magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf)); - } else { - magErr = 1.0 / magErr; - } - } + // XXX note that this mag is either calibrated or instrumental depending on existence of zero point + float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point + // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired + float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point + if (repairMagErrors) { + // I need to add the kron error in quadrature becasue pet_error ignores the object flux + float Krf = source->moments ? source->moments->KronFlux : NAN; + float dKrf = source->moments ? source->moments->KronFluxErr : NAN; + if (isfinite (Krf) && isfinite (dKrf)) { + magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf)); + } else { + magErr = 1.0 / magErr; + } + } psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", mag); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", magErr); @@ -705,61 +1370,61 @@ // Flux Annuli (if we have extended source measurements, we have these. only optionally save them) if (doAnnuli) { - psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); - psVectorInit (radSB, NAN); - psVectorInit (radFlux, NAN); - psVectorInit (radFill, NAN); - if (!source->extpars) goto empty_annuli; - if (!source->extpars->radProfile) goto empty_annuli; - if (!source->extpars->radProfile->binSB) goto empty_annuli; - psAssert (source->extpars->radProfile->binSum, "programming error"); - psAssert (source->extpars->radProfile->binFill, "programming error"); - psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { - radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; - radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; - radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; - } - - empty_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radSB); - psFree (radFlux); - psFree (radFill); - } - if (nRow < 0) { - nRow = row->list->n; - } else { - psAssert (nRow == row->list->n, "inconsistent row lengths"); - } - psArrayAdd (table, 100, row); - psFree (row); + psVector *radSB = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32); + psVectorInit (radSB, NAN); + psVectorInit (radFlux, NAN); + psVectorInit (radFill, NAN); + if (!source->extpars) goto empty_annuli; + if (!source->extpars->radProfile) goto empty_annuli; + if (!source->extpars->radProfile->binSB) goto empty_annuli; + psAssert (source->extpars->radProfile->binSum, "programming error"); + psAssert (source->extpars->radProfile->binFill, "programming error"); + psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) { + radSB->data.F32[j] = source->extpars->radProfile->binSB->data.F32[j]; + radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j]; + radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j]; + } + + empty_annuli: + psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radSB); + psFree (radFlux); + psFree (radFill); + } + if (nRow < 0) { + nRow = row->list->n; + } else { + psAssert (nRow == row->list->n, "inconsistent row lengths"); + } + psArrayAdd (table, 100, row); + psFree (row); } if (table->n == 0) { - if (!psFitsWriteBlank (fits, outhead, extname)) { - psError(psErrorCodeLast(), false, "Unable to write empty sources file."); - psFree(outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; + psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; } psFree (outhead); @@ -885,25 +1550,25 @@ extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); - // magErr may have been saved in inverted form + // magErr may have been saved in inverted form float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); if (isfinite(mag)) { - extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); - if (isfinite(magErr)) { - if (repairMagErrors) { - // I need to add the kron error in quadrature becasue pet_error ignores the object flux - float Krf = source->moments ? source->moments->KronFlux : NAN; - float dKrf = source->moments ? source->moments->KronFluxErr : NAN; - if (isfinite (Krf) && isfinite (dKrf)) { - magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf)); - } else { - magErr = 1.0 / magErr; - } - extpars->petrosianFluxErr = extpars->petrosianFlux * magErr; - } else { - extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; - } - } - } + extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); + if (isfinite(magErr)) { + if (repairMagErrors) { + // I need to add the kron error in quadrature becasue pet_error ignores the object flux + float Krf = source->moments ? source->moments->KronFlux : NAN; + float dKrf = source->moments ? source->moments->KronFluxErr : NAN; + if (isfinite (Krf) && isfinite (dKrf)) { + magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf)); + } else { + magErr = 1.0 / magErr; + } + extpars->petrosianFluxErr = extpars->petrosianFlux * magErr; + } else { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + } extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); @@ -973,9 +1638,9 @@ int nParamMax = 0; for (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // this is the source associated with this image pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; if (source->modelFits == NULL) continue; @@ -1004,6 +1669,6 @@ pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // XXX if no model fits are saved, write out modelEXT? @@ -1018,9 +1683,9 @@ // pmSourceExtFitPars *extPars = source->extFitPars->data[j]; - // assert (extPars); - - // skip models which were not actually fitted - // XXX - if (model->flags & badModel) continue; + // assert (extPars); + + // skip models which were not actually fitted + // XXX + if (model->flags & badModel) continue; PAR = model->params->data.F32; @@ -1029,37 +1694,37 @@ yPos = PAR[PM_PAR_YPOS]; - // for the extended source models, we do not always fit the centroid in the non-linear fitting process - // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, - // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) - // TRAIL : X,Y are fitted - // - - // XXX this should be based on what happened, not on the model type - if (model->type == modelTypeTrail) { - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - } else { - // this is definitely an underestimate since it does not - // account for the extent of the source - xErr = fwhmMajor * model->magErr / 2.35; - yErr = fwhmMinor * model->magErr / 2.35; - } - - @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; - @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; - @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; - @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); - @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; - @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; - - float kronFlux = source->moments ? source->moments->KronFlux : NAN; - float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; + // for the extended source models, we do not always fit the centroid in the non-linear fitting process + // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, + // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) + // TRAIL : X,Y are fitted + // + + // XXX this should be based on what happened, not on the model type + if (model->type == modelTypeTrail) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // this is definitely an underestimate since it does not + // account for the extent of the source + xErr = fwhmMajor * model->magErr / 2.35; + yErr = fwhmMinor * model->magErr / 2.35; + } + + @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; + @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0; + @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0; + @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos); + @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD; + @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0; + + float kronFlux = source->moments ? source->moments->KronFlux : NAN; + float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN; row = psMetadataAlloc (); - // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation - // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" 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); @@ -1070,21 +1735,21 @@ @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT", 0, "EXT model ra coordinate", raPos); @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT", 0, "EXT model dec coordinate", decPos); - @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); + @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN; + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX", 0, "EXT fit instrumental counts", instFlux); 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); - @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); - @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); - @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); - @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE", PS_DATA_S32, "type for chosen EXT_MODEL", source->modelEXT ? source->modelEXT->type : -1); - - // EAM : adding for PV2 outputs: - @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); + @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; + @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); + @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE", PS_DATA_S32, "type for chosen EXT_MODEL", source->modelEXT ? source->modelEXT->type : -1); + + // EAM : adding for PV2 outputs: + @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags); @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); - @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); + @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE", 0, "plate scale at source (arcsec/pixel)", pltScale); // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX", 0, "second moment in x", extPars->Mxx); @@ -1103,29 +1768,29 @@ // XXX these should be major and minor, not 'x' and 'y' - if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", PAR[PM_PAR_SIGMA]); // this is not fitted - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", NAN); // this is not fitted, so error is NAN - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); - } else { - if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); - } else { - psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_SXX]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); - } - } + if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", PAR[PM_PAR_SIGMA]); // this is not fitted + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", PAR[PM_PAR_THETA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", NAN); // this is not fitted, so error is NAN + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } else { + if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY]) || !isfinite(PAR[PM_PAR_SXY])) { + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]); + } else { + psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width (major axis), length for trail", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width (minor axis), sigma for trail", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)", dPAR[PM_PAR_SXX]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SYY]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_SXY]); + } + } // write out the other generic parameters @@ -1140,5 +1805,5 @@ snprintf (name, 64, "EXT_PAR_%02d", k); - + if (k < model->params->n) { psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); @@ -1148,15 +1813,15 @@ } - // optionally, write out the covariance matrix values - // XXX do I need to pad this to match the biggest covar matrix? - if (false && model->covar) { - for (int iy = 0; iy < model->covar->numCols; iy++) { - for (int ix = iy; ix < model->covar->numCols; ix++) { - snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); - psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); - - } - } - } + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (false && model->covar) { + for (int iy = 0; iy < model->covar->numCols; iy++) { + for (int ix = iy; ix < model->covar->numCols; ix++) { + snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]); + + } + } + } psArrayAdd (table, 100, row); psFree (row); @@ -1280,28 +1945,28 @@ PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); // XXX add an error: - // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); - } - - // NOTE: we no longer write out the covariance matrix - if (false) { - // read the covariance matrix - int nparams = model->params->n; - psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); - for (int y = 0; y < nparams; y++) { - for (int x = 0; x < nparams; x++) { - char name[64]; - snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); - covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); - } - } - model->covar = covar; - } - - // we are only saving the values stored in dPAR[SXX,etc] + // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); + } + + // NOTE: we no longer write out the covariance matrix + if (false) { + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + } + + // we are only saving the values stored in dPAR[SXX,etc] dPAR[PM_PAR_SXX] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ_ERR"); dPAR[PM_PAR_SYY] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN_ERR"); dPAR[PM_PAR_SXY] = psMetadataLookupF32(&status, row, "EXT_THETA_ERR"); - // other parameters that we need to read + // other parameters that we need to read PAR[PM_PAR_SKY] = psMetadataLookupF32(&status, row, "SKY_EXT"); @@ -1341,6 +2006,6 @@ // perform full non-linear fits / extended source analysis? if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { - psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n"); + return true; } @@ -1385,77 +2050,77 @@ for (int i = 0; i < sources->n; i++) { - // this is the source associated with this image + // this is the source associated with this image pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; // skip sources without radial aper measurements (or insufficient) - if (source->radialAper == NULL) continue; + if (source->radialAper == NULL) continue; // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set"); - for (int entry = 0; entry < source->radialAper->n; entry++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; - assert (radialAper); - - if (pmSourcePositionUseMoments(source)) { - xPos = source->moments->Mx; - yPos = source->moments->My; - } else { - xPos = source->peak->xf; - yPos = source->peak->yf; - } - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation - // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment" - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER", 0, "Center of aperture measurements", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); - if (fwhmValues) { - psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); - } - - // XXX if we have raw radial apertures, write them out here - psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32); - psVectorInit (radFlux, NAN); - psVectorInit (radFluxErr, NAN); - psVectorInit (radFill, NAN); - if (!radialAper->flux) goto write_annuli; - if (!radialAper->fill) goto write_annuli; - psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths"); - psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths"); - - // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux - for (int j = 0; j < radialAper->flux->n; j++) { - radFlux->data.F32[j] = radialAper->flux->data.F32[j]; - radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; - radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j]; - radFill->data.F32[j] = radialAper->fill->data.F32[j]; - } - - write_annuli: - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); - psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFluxStdev); - psFree (radFill); - - psArrayAdd (table, 100, row); - psFree (row); - } + for (int entry = 0; entry < source->radialAper->n; entry++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmSourceRadialApertures *radialAper = source->radialAper->data[entry]; + assert (radialAper); + + if (pmSourcePositionUseMoments(source)) { + xPos = source->moments->Mx; + yPos = source->moments->My; + } else { + xPos = source->peak->xf; + yPos = source->peak->yf; + } + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation + // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment" + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER", 0, "Center of aperture measurements", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER", 0, "Center of aperture measurements", yPos); + if (fwhmValues) { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "FWHM of matched PSF", fwhmValues->data.F32[entry]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM", 0, "image is not FWHM-matched", NAN); + } + + // XXX if we have raw radial apertures, write them out here + psVector *radFlux = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFill = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32); + psVectorInit (radFlux, NAN); + psVectorInit (radFluxErr, NAN); + psVectorInit (radFill, NAN); + if (!radialAper->flux) goto write_annuli; + if (!radialAper->fill) goto write_annuli; + psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths"); + psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths"); + + // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux + for (int j = 0; j < radialAper->flux->n; j++) { + radFlux->data.F32[j] = radialAper->flux->data.F32[j]; + radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j]; + radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j]; + radFill->data.F32[j] = radialAper->fill->data.F32[j]; + } + + write_annuli: + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFluxStdev); + psFree (radFill); + + psArrayAdd (table, 100, row); + psFree (row); + } } @@ -1581,6 +2246,6 @@ // perform full non-linear fits / extended source analysis? if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { - psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); + return true; } @@ -1609,11 +2274,11 @@ pmSource *thisSource = sources->data[i]; - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - - // if we did not fit the galaxy model, modelFits will be NULL + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // if we did not fit the galaxy model, modelFits will be NULL if (source->modelFits == NULL) continue; - // if we did not fit the galaxy model, galaxyFits will also be NULL + // if we did not fit the galaxy model, galaxyFits will also be NULL if (source->galaxyFits == NULL) continue; Index: /trunk/psModules/src/objects/pmSourceIO_Ghosts.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_Ghosts.c (revision 41891) +++ /trunk/psModules/src/objects/pmSourceIO_Ghosts.c (revision 41892) @@ -136,5 +136,5 @@ pmFPAview *view = pmFPAviewAlloc (0); - //We need to check whether we are dealing with an old style ghost_model, or a new style model. Check if the mirror_rad polynomial exists + // We need to check whether we are dealing with an old style ghost_model, or a new style model. Check if the mirror_rad polynomial exists float mirCheck = 0; md = psMetadataLookupMetadata (&status, ghostModel, "GHOST.MIRROR.RAD"); @@ -153,5 +153,4 @@ GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor); - // select the input astrometry data (also carries the refstars) pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); @@ -211,5 +210,6 @@ double theta0 = atan2(ref->FP->y,ref->FP->x); - if(mirCheck) { + // EAM: XXX we should just use the existence of mirrorRad (!= NULL) instead of carrying a different bool + if (mirCheck) { //TdB: first mirror the reference star positions (around the 0,0 pixel) using the radial offset coefficients and the ghost/star angle double ghost_offset_rad = psPolynomial1DEval (mirrorRad, rSrc); @@ -220,6 +220,5 @@ ghost->FP->x = ghost_x_fpa_mirror + psPolynomial2DEval(centerX, ghost_x_fpa_mirror, ghost_y_fpa_mirror); ghost->FP->y = ghost_y_fpa_mirror + psPolynomial2DEval(centerY, ghost_x_fpa_mirror, ghost_y_fpa_mirror); - } - if(!mirCheck) { + } else { //Use the old-style ghost position determination ghost->FP->x = -ref->FP->x + psPolynomial2DEval(centerX, -ref->FP->x, -ref->FP->y); @@ -303,4 +302,5 @@ psFree (outerMajor); psFree (outerMinor); + psFree (mirrorRad); psFree (ghostModel); psFree (view); @@ -327,4 +327,5 @@ psFree (outerMajor); psFree (outerMinor); + psFree (mirrorRad); psFree (ghostModel); psFree (view);