Index: /trunk/psModules/src/objects/Makefile.am =================================================================== --- /trunk/psModules/src/objects/Makefile.am (revision 34258) +++ /trunk/psModules/src/objects/Makefile.am (revision 34259) @@ -80,5 +80,6 @@ models/pmModel_SERSIC.c \ models/pmModel_EXP.c \ - models/pmModel_DEV.c + models/pmModel_DEV.c \ + models/pmModel_TRAIL.c pkginclude_HEADERS = \ @@ -126,10 +127,11 @@ models/pmModel_SERSIC.h \ models/pmModel_EXP.h \ - models/pmModel_DEV.h + models/pmModel_DEV.h \ + models/pmModel_TRAIL.h CLEANFILES = *~ # pmSourceID_CMF_* functions use a common framework -BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c pmSourceIO_CMF_PS1_V4.c +BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c pmSourceIO_CMF_PS1_DV1.c pmSourceIO_CMF_PS1_DV2.c pmSourceIO_CMF_PS1_SV1.c pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl @@ -145,3 +147,12 @@ mksource.pl pmSourceIO_CMF.c.in PS1_V4 pmSourceIO_CMF_PS1_V4.c +pmSourceIO_CMF_PS1_DV1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_DV1 pmSourceIO_CMF_PS1_DV1.c + +pmSourceIO_CMF_PS1_DV2.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_DV2 pmSourceIO_CMF_PS1_DV2.c + +pmSourceIO_CMF_PS1_SV1.c : pmSourceIO_CMF.c.in mksource.pl + mksource.pl pmSourceIO_CMF.c.in PS1_SV1 pmSourceIO_CMF_PS1_SV1.c + # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in Index: /trunk/psModules/src/objects/mksource.pl =================================================================== --- /trunk/psModules/src/objects/mksource.pl (revision 34258) +++ /trunk/psModules/src/objects/mksource.pl (revision 34259) @@ -1,3 +1,4 @@ #!/usr/bin/env perl +$DEBUG = 0; # this program takes the pmSourceIO_CMF.in.c template file and generates the .c version based on the given I/O format made @@ -17,8 +18,6 @@ "PS1_V2", 2, "PS1_V3", 3, - "PS1_V4", 4); - -print "1: $cmfmodes{1}\n"; -print "PS1_V1: $cmfmodes{'PS1_V1'}\n"; + "PS1_V4", 4, + ); open (FILE, "$template") || die "failed to open template $template\n"; @@ -35,4 +34,15 @@ # @ MODE +# XXX need to add features: split @foo,bar,baz@ by commas +# treat each chunk as a rule +# add the following options: +# !MODE -- exclude the given entry (defaults to all, or is ALL required?) +# * and ? regexp + +# some examples: +# @ALL,!PS1_V1@ +# @PS1_DV?@ +# @PS1_V?,!PS1_V1@ + foreach $line (@list) { @@ -40,36 +50,123 @@ $line =~ s|\@CMFMODE\@|$cmfmode|g; - if ($line =~ m|\@ALL\@|) { - $line =~ s|\@ALL\@\s*||; + # print and continue if we do not match @RULES@ + unless ($line =~ m|\@.*\@|) { + print "no rule\n" if $DEBUG; + print FILE $line; + next; } - ($isMode) = ($line =~ m|\@=(\S*)\@|); - ($gtMode) = ($line =~ m|\@>(\S*)\@|); - ($ltMode) = ($line =~ m|\@<(\S*)\@|); + # grab the rules and the rest of the line + ($prefix,$rules,$content) = ($line =~ m|(.*)\@(.*)\@\s*(.*)|); - if ($isMode) { - if ($isMode ne $cmfmode) { next; } - $line =~ s|\@=\S*\@\s*||; + # split the rules into separate items + @rules = split (",", $rules); + + $keepLine = 0; + # does $cmfmode match any of the rules? + foreach $rule (@rules) { + print "rule: $rule\n" if $DEBUG; + + # special rule "ALL" + if ($rule eq "ALL") { + print "ALL match\n" if $DEBUG; + $keepLine = 1; + next; + } # look for other rules (esp !foo) + + # pure match + if ($cmfmode eq $rule) { + print "simple match\n" if $DEBUG; + $keepLine = 1; + next; + } # skip to end? + + # NOT match + if ($rule =~ m|^!|) { + print "NOT rule: $rule\n" if $DEBUG; + ($realrule) = ($rule =~ m|^!(.*)|); + if ($cmfmode eq $realrule) { $keepLine = 0; } # skip to end? + next; + } + + # simple regexp: foo* + if ($rule =~ m|\*$|) { + print "regexp * rule: $rule\n" if $DEBUG; + ($realrule) = ($rule =~ m|(.*)\*$|); + if ($cmfmode =~ m|$realrule|) { $keepLine = 1; } # skip to end? + next; + } + + # simple regexp: foo? + if ($rule =~ m|\?$|) { + print "regexp ? rule: $rule\n" if $DEBUG; + ($realrule) = ($rule =~ m|(.*)\?$|); + if (substr($cmfmode,0,-1) eq $realrule) { $keepLine = 1; } # skip to end? + next; + } + + # rule: =FOO + if ($rule =~ m|^=|) { + print "= rule: $rule\n" if $DEBUG; + $realrule = substr($cmfmode,1); + if ($cmfmode eq $realrule) { $keepLine = 1; } # skip to end? + next; + } + + # only apply the < > <= >= rules if cmfmode is one of cmfmodes + # rule: >=FOO + if ($rule =~ m|^>=|) { + print ">= rule: $rule\n" if $DEBUG; + if ($cmfmodes{$cmfmode} == 0) { next; } + $realrule = substr($rule,2); + $thisLevel = $cmfmodes{$realrule}; + $myLevel = $cmfmodes{$cmfmode}; + print "levels: $cmfmode, $realrule, $myLevel, $thisLevel\n" if $DEBUG; + if ($myLevel >= $thisLevel) { $keepLine = 1; } + next; + } + + # rule: >FOO + if ($rule =~ m|^>|) { + print "> rule: $rule\n" if $DEBUG; + if ($cmfmodes{$cmfmode} == 0) { next; } + $realrule = substr($rule,1); + $thisLevel = $cmfmodes{$realrule}; + $myLevel = $cmfmodes{$cmfmode}; + print "levels: $cmfmode, $realrule, $myLevel, $thisLevel\n" if $DEBUG; + if ($myLevel > $thisLevel) { $keepLine = 1; } + next; + } + + # rule: <=FOO + if ($rule =~ m|^<=|) { + print "<= rule: $rule\n" if $DEBUG; + if ($cmfmodes{$cmfmode} == 0) { next; } + $realrule = substr($rule,2); + $thisLevel = $cmfmodes{$realrule}; + $myLevel = $cmfmodes{$cmfmode}; + print "levels: $cmfmode, $realrule, $myLevel, $thisLevel\n" if $DEBUG; + if ($myLevel <= $thisLevel) { $keepLine = 1; } + next; + } + + # rule: \S*\@\s*||; + if ($keepLine) { + print FILE "$prefix $content\n"; } - - if ($ltMode) { - # print "ltMode : $line\n"; - $thisLevel = $cmfmodes{$ltMode}; - $myLevel = $cmfmodes{$cmfmode}; - print "ltMode : $ltMode vs $cmfmode, $thisLevel, $myLevel\n"; - if ($myLevel >= $thisLevel) { next; } - $line =~ s|\@<\S*\@\s*||; - } - - print FILE $line; } Index: /trunk/psModules/src/objects/models/pmModel_TRAIL.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 34259) +++ /trunk/psModules/src/objects/models/pmModel_TRAIL.c (revision 34259) @@ -0,0 +1,358 @@ +/****************************************************************************** + * this file defines the TRAIL source shape model. This represents an infinitely thin + * line of length convolved with a Gaussian PSF. The models use a psVector to represent + * the set of parameters, with the sequence used to specify the meaning of the parameter. + * The meaning of the parameters may thus vary depending on the specifics of the model. + * All models which are used as a PSF representations share a few parameters, for which # + * define names are listed in pmModel.h: + + pure Gaussian: + exp(-z) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - flux normalization + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_LENGTH 4 - trail length + * PM_PAR_THETA 5 - position angle + * PM_PAR_SIGMA 6 - PSF Gaussian sigma (not fitted?) + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmModelClass.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + +#include "pmModel_TRAIL.h" + +# define PM_MODEL_NPARAM 7 +# define PM_MODEL_FUNC pmModelFunc_TRAIL +# define PM_MODEL_FLUX pmModelFlux_TRAIL +# define PM_MODEL_GUESS pmModelGuess_TRAIL +# define PM_MODEL_LIMITS pmModelLimits_TRAIL +# define PM_MODEL_RADIUS pmModelRadius_TRAIL +# define PM_MODEL_FROM_PSF pmModelFromPSF_TRAIL +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_TRAIL +# define PM_MODEL_FIT_STATUS pmModelFitStatus_TRAIL +# define PM_MODEL_SET_LIMITS pmModelSetLimits_TRAIL + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2, 0.5, -3.3, -0.5 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e+8, +1.0e4, +1.0e4, 150.0, +3.3 , 5.0 }; + +// Moderate parameter limits +static float *paramsMinModerate = paramsMinLax; +static float *paramsMaxModerate = paramsMaxLax; + +// Strict parameter limits +static float *paramsMinStrict = paramsMinLax; +static float *paramsMaxStrict = paramsMaxLax; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 2.0, 0.1, 0.1 }; + +static bool limitsApply = true; // Apply limits? + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords +psF32 PM_MODEL_FUNC(psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 ST = sin(PAR[PM_PAR_THETA]); + psF32 CT = cos(PAR[PM_PAR_THETA]); + + psF32 S2 = 2.0 * PS_SQR(PAR[PM_PAR_SIGMA]); + + psF32 Zp = (X*CT + Y*ST + 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + psF32 Zm = (X*CT + Y*ST - 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + + // psF32 Zp = (X*CT + Y*ST + 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + // psF32 Zm = (X*CT + Y*ST - 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + + psF32 Ep = erf(Zp); + psF32 Em = erf(Zm); + + psF32 Rxy = Y*CT - X*ST; + psF32 Gxy = exp(-Rxy*Rxy/S2); + + psF32 Pxy = Gxy * (Ep - Em); + psF32 f = Pxy * PAR[PM_PAR_I0] + PAR[PM_PAR_SKY]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + dPAR[PM_PAR_SKY] = 1.0; + dPAR[PM_PAR_I0] = Pxy; + + float dGdR = -2.0 * Rxy * Gxy / S2; // -R Gxy / (2 Sigma^2) + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dRdX = +ST; + float dRdY = -CT; + float dRdT = -Y*ST - X*CT; + + float dGdX = dGdR * dRdX; + float dGdY = dGdR * dRdY; + float dGdT = dGdR * dRdT; + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dZpdX = -CT / sqrt(S2); + float dZmdX = dZpdX; // float dZmdX = -CT / sqrt(S2); dZmdX = dZpdX + + float dZpdY = -ST / sqrt(S2); // float dZmdY = -ST / sqrt(S2); dZmdY = dZpdY + float dZmdY = dZpdY; + + float dZpdL = +0.5 / sqrt(S2); + float dZmdL = -0.5 / sqrt(S2); + + float dZpdT = (-X*ST + Y*CT) / sqrt(S2); + float dZmdT = dZpdT; // dZpdT = dZmdT + + float dEdZp = exp (-Zp*Zp) * M_2_SQRTPI; + float dEdZm = exp (-Zm*Zm) * M_2_SQRTPI; + + float dEpdX = dEdZp * dZpdX; + float dEmdX = dEdZm * dZmdX; + + float dEpdY = dEdZp * dZpdY; + float dEmdY = dEdZm * dZmdY; + + float dEpdL = dEdZp * dZpdL; + float dEmdL = dEdZm * dZmdL; + + float dEpdT = dEdZp * dZpdT; + float dEmdT = dEdZm * dZmdT; + + float dPdX = dGdX * (Ep - Em) + Gxy * (dEpdX - dEmdX); + float dPdY = dGdY * (Ep - Em) + Gxy * (dEpdY - dEmdY); + + // dGdL is 0.0 because dRdL is 0.0 + float dPdL = Gxy * (dEpdL - dEmdL); + + float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT); + + dPAR[PM_PAR_XPOS] = PAR[PM_PAR_I0] * dPdX; + dPAR[PM_PAR_YPOS] = PAR[PM_PAR_I0] * dPdY; + + dPAR[PM_PAR_LENGTH] = PAR[PM_PAR_I0] * dPdL; + dPAR[PM_PAR_THETA] = PAR[PM_PAR_I0] * dPdT; + dPAR[PM_PAR_SIGMA] = 0; // we don't actually allow this to vary, so we do not need to calculate it + } + return(f); +} + +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 + +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + +// make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters +bool PM_MODEL_GUESS (pmModel *model, pmSource *source) +{ + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 + PAR[PM_PAR_SKY] = 0.0; + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape psfShape; + psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); + + psEllipseMoments emoments; + emoments.x2 = source->moments->Mxx; + emoments.xy = source->moments->Mxy; + emoments.y2 = source->moments->Myy; + + // force the axis ratio to be < 20.0 + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); + + if (!isfinite(axes.major)) return false; + if (!isfinite(axes.minor)) return false; + if (!isfinite(axes.theta)) return false; + + float size = (axes.major > source->moments->Mrf) ? axes.major : source->moments->Mrf; + + // axes.major is a sigma in the major direction; scale to + PAR[PM_PAR_LENGTH] = 1.5*2.35*size; // a tophat of length L has L = 1.5 * 2.35 * sigma + PAR[PM_PAR_THETA] = axes.theta; // theta in radians + PAR[PM_PAR_SIGMA] = psfAxes.major; // psf major axes (sigma of the psf) + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + return(true); +} + +psF64 PM_MODEL_FLUX (const psVector *params) +{ + psF32 *PAR = params->data.F32; + psF64 Flux = PAR[PM_PAR_I0] * PAR[PM_PAR_LENGTH] * PAR[PM_PAR_SIGMA] * 2.0 * sqrt(2.0 * M_PI); + return(Flux); +} + +// return the radius which yields the requested flux +// this function is never allowed to return <= 0 +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF32 *PAR = params->data.F32; + + // PAR_LENGTH is the unconvolved length. add a bit for safety + return (0.5*PAR[PM_PAR_LENGTH] + 2); +} + +// construct the PSF model from the FLT model and the psf +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) +{ + psWarning ("do you really want to use a trail as a PSF model??"); + return false; +} + +// generate a model based on a the psf model +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psWarning ("do you really want to use a trail as a PSF model??"); + return false; +} + +// check the status of the fitted model +// this test is invalid if the parameters are derived +// from the PSF model +// XXX how is this used? it prevents forced photometry from ever being 'successful' +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: /trunk/psModules/src/objects/models/pmModel_TRAIL.h =================================================================== --- /trunk/psModules/src/objects/models/pmModel_TRAIL.h (revision 34259) +++ /trunk/psModules/src/objects/models/pmModel_TRAIL.h (revision 34259) @@ -0,0 +1,15 @@ +#ifndef PM_MODEL_TRAIL_H + +#include "pmModel.h" + +psF32 pmModelFunc_TRAIL(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_TRAIL(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_TRAIL(pmModel *model, pmSource *source); +psF64 pmModelFlux_TRAIL(const psVector *params); +psF64 pmModelRadius_TRAIL(const psVector *params, psF64 flux); +bool pmModelFromPSF_TRAIL(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_TRAIL(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_TRAIL(pmModel *model); +void pmModelSetLimits_TRAIL(pmModelLimitsType type); + +#endif Index: /trunk/psModules/src/objects/pmModelClass.c =================================================================== --- /trunk/psModules/src/objects/pmModelClass.c (revision 34258) +++ /trunk/psModules/src/objects/pmModelClass.c (revision 34259) @@ -51,4 +51,5 @@ # include "models/pmModel_EXP.h" # include "models/pmModel_DEV.h" +# include "models/pmModel_TRAIL.h" static pmModelClass defaultModels[] = { @@ -59,6 +60,7 @@ {"PS_MODEL_RGAUSS", 8, (pmModelFunc)pmModelFunc_RGAUSS, (pmModelFlux)pmModelFlux_RGAUSS, (pmModelRadius)pmModelRadius_RGAUSS, (pmModelLimits)pmModelLimits_RGAUSS, (pmModelGuessFunc)pmModelGuess_RGAUSS, (pmModelFromPSFFunc)pmModelFromPSF_RGAUSS, (pmModelParamsFromPSF)pmModelParamsFromPSF_RGAUSS, (pmModelFitStatusFunc)pmModelFitStatus_RGAUSS, (pmModelSetLimitsFunc)pmModelSetLimits_RGAUSS }, {"PS_MODEL_SERSIC", 8, (pmModelFunc)pmModelFunc_SERSIC, (pmModelFlux)pmModelFlux_SERSIC, (pmModelRadius)pmModelRadius_SERSIC, (pmModelLimits)pmModelLimits_SERSIC, (pmModelGuessFunc)pmModelGuess_SERSIC, (pmModelFromPSFFunc)pmModelFromPSF_SERSIC, (pmModelParamsFromPSF)pmModelParamsFromPSF_SERSIC, (pmModelFitStatusFunc)pmModelFitStatus_SERSIC, (pmModelSetLimitsFunc)pmModelSetLimits_SERSIC }, - {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, - {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, + {"PS_MODEL_EXP", 7, (pmModelFunc)pmModelFunc_EXP, (pmModelFlux)pmModelFlux_EXP, (pmModelRadius)pmModelRadius_EXP, (pmModelLimits)pmModelLimits_EXP, (pmModelGuessFunc)pmModelGuess_EXP, (pmModelFromPSFFunc)pmModelFromPSF_EXP, (pmModelParamsFromPSF)pmModelParamsFromPSF_EXP, (pmModelFitStatusFunc)pmModelFitStatus_EXP, (pmModelSetLimitsFunc)pmModelSetLimits_EXP }, + {"PS_MODEL_DEV", 7, (pmModelFunc)pmModelFunc_DEV, (pmModelFlux)pmModelFlux_DEV, (pmModelRadius)pmModelRadius_DEV, (pmModelLimits)pmModelLimits_DEV, (pmModelGuessFunc)pmModelGuess_DEV, (pmModelFromPSFFunc)pmModelFromPSF_DEV, (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV, (pmModelFitStatusFunc)pmModelFitStatus_DEV, (pmModelSetLimitsFunc)pmModelSetLimits_DEV }, + {"PS_MODEL_TRAIL", 7, (pmModelFunc)pmModelFunc_TRAIL, (pmModelFlux)pmModelFlux_TRAIL, (pmModelRadius)pmModelRadius_TRAIL, (pmModelLimits)pmModelLimits_TRAIL, (pmModelGuessFunc)pmModelGuess_TRAIL, (pmModelFromPSFFunc)pmModelFromPSF_TRAIL, (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL, (pmModelFitStatusFunc)pmModelFitStatus_TRAIL, (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL }, }; Index: /trunk/psModules/src/objects/pmModelFuncs.h =================================================================== --- /trunk/psModules/src/objects/pmModelFuncs.h (revision 34258) +++ /trunk/psModules/src/objects/pmModelFuncs.h (revision 34259) @@ -81,4 +81,9 @@ #define PM_PAR_8 8 ///< Model-dependent parameter +// these are used by pmModel_TRAIL, with refers to L and Theta explicitly +#define PM_PAR_LENGTH 4 ///< trail length +#define PM_PAR_THETA 5 ///< position angle +#define PM_PAR_SIGMA 6 ///< position angle + /*** these prototype classes are used to define elements of the pmModelClass structure below ***/ Index: /trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /trunk/psModules/src/objects/pmSourceFitModel.c (revision 34258) +++ /trunk/psModules/src/objects/pmSourceFitModel.c (revision 34259) @@ -182,4 +182,11 @@ constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 1; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + break; + case PM_SOURCE_FIT_TRAIL: + // special mode for pmModel_TRAIL: Io, Xo, Yo, Length, and Theta (not Io or Sigma) + nParams = params->n - 3; + psVectorInit (constraint->paramMask, 0); + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; + constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SIGMA] = 1; break; case PM_SOURCE_FIT_INDEX: Index: /trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /trunk/psModules/src/objects/pmSourceFitModel.h (revision 34258) +++ /trunk/psModules/src/objects/pmSourceFitModel.h (revision 34259) @@ -22,4 +22,5 @@ PM_SOURCE_FIT_INDEX, PM_SOURCE_FIT_NO_INDEX, + PM_SOURCE_FIT_TRAIL, } pmSourceFitMode; Index: /trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 34258) +++ /trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 34259) @@ -89,5 +89,9 @@ // 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++) { - pmSource *source = (pmSource *) sources->data[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 @@ -106,78 +110,107 @@ 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 (); - @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); + @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); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - @>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@ 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", outputs.apRadius); - @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_SV1,PS1_DV2@ 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", outputs.apRadius); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); + + @PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - - @>PS1_V2@ 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@ 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@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - @>PS1_V2@ 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@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 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); - - @ALL@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); - - // XXX not sure how to get this : need to load Nimages with weight? - @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); + @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); + @ALL,!PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); + + @>=PS1_V3@ 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@ 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@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); + @>PS1_V2,PS1_SV1,PS1_DV2@ 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_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); + @>PS1_V2,PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); + + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); + @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); + @>PS1_V2,PS1_SV1,PS1_DV2@ 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_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_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); + @PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); + @PS1_DV2@ 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_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); + @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2", PS_DATA_S32, "more padding", 0); + @PS1_SV1@ + + // 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); @@ -282,9 +315,10 @@ @ALL@ source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); @ALL@ source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - @>PS1_V2@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->apMagRaw = psMetadataLookupF32 (&status, row, "AP_MAG_RAW"); // XXX use these to determine PAR[PM_PAR_I0] if they exist? - @>PS1_V2@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - @>PS1_V2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); + // XXX add these to PS1_SV1? + @>PS1_V2,PS1_SV1,PS1_DV?@ source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); + @>PS1_V2,PS1_SV1,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); // XXX this scaling is incorrect: does not include the 2 \pi AREA factor @@ -307,5 +341,5 @@ @ALL@ source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - @>PS1_V2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); @ALL@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); @ALL@ source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); @@ -325,11 +359,23 @@ @ALL@ source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - @>PS1_V2@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); - @>PS1_V2@ source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); - @>PS1_V2@ source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); - @>PS1_V2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); - - @>PS1_V2@ source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); - @>PS1_V2@ source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); + // 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_SV1@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); + @>PS1_V2,PS1_SV1@ source->moments->Mxxy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxyy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); + + @>PS1_V2,PS1_SV1@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); + @>PS1_V2,PS1_SV1@ source->moments->Mxxxy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxxyy = 0.0; + @>PS1_V2,PS1_SV1@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); + @>PS1_V2,PS1_SV1@ source->moments->Myyyy = 0.0; + + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); @>PS1_V3@ source->skyRadius = psMetadataLookupF32 (&status, row, "SKY_LIMIT_RAD"); @@ -337,19 +383,21 @@ @>PS1_V3@ source->skySlope = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE"); - // 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@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); - @>PS1_V2@ source->moments->Mxxy = 0.0; - @>PS1_V2@ source->moments->Mxyy = 0.0; - @>PS1_V2@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); - - @>PS1_V2@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); - @>PS1_V2@ source->moments->Mxxxy = 0.0; - @>PS1_V2@ source->moments->Mxxyy = 0.0; - @>PS1_V2@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); - @>PS1_V2@ source->moments->Myyyy = 0.0; - - @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - @>PS1_V2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); + @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_DV2@ source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); + @PS1_DV2@ source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); + @PS1_DV2@ source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); + @PS1_DV2@ source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); + @PS1_DV?@ } + + @ALL@ source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); + @>PS1_V2,PS1_SV1,PS1_DV2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); assert (status); @@ -409,13 +457,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, "min 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, "min 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++) { - pmSource *source = sources->data[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; // skip sources without measurements @@ -470,4 +522,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); } else { psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); @@ -479,4 +532,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); } } @@ -539,10 +593,104 @@ psError(psErrorCodeLast(), false, "writing ext data %s\n", extname); psFree (outhead); - psFree(table); - return false; + psFree(table); + return false; } psFree (outhead); psFree (table); + return true; +} + +bool pmSourcesRead_CMF_@CMFMODE@_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); + return false; + } + + // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point + float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); + float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); + float magOffset = zeropt + 2.5*log10(exptime); + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + pmSource *source = NULL; +#ifndef ASSUME_SORTED + long j = seq < sources->n ? seq : sources->n - 1; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } +#else + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + source = sources->data[j]; +#endif + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + + if (!source->extpars) { + source->extpars = pmSourceExtendedParsAlloc (); + } + pmSourceExtendedPars *extpars = source->extpars; + + // Assume that X_EXT Y_EXT and sigmas match the psf src so skip + + // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by + // psf size or something? + extpars->axes.major = 1.0; + extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); + extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); + + 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)) { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + + extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); + extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); + extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); + extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); + extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); + extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); + extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); + + psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); + psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); + psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); + + if (radSB && radSB->n > 0) { + extpars->radProfile = pmSourceRadialProfileAlloc(); + extpars->radProfile->binSB = psMemIncrRefCounter(radSB); + extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); + extpars->radProfile->binFill = psMemIncrRefCounter(radFill); + } + + psFree(row); + } + return true; } @@ -572,5 +720,10 @@ int nParamMax = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[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->modelFits == NULL) continue; for (int j = 0; j < source->modelFits->n; j++) { @@ -586,5 +739,8 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + pmSource *thisSource = sources->data[i]; + + // 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? @@ -641,7 +797,19 @@ // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + 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]); + 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)", dPAR[PM_PAR_SIGMA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } else { + 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_LENGTH]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)", dPAR[PM_PAR_SIGMA]); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR", 0, "EXT orientation angle (error)", dPAR[PM_PAR_THETA]); + } // write out the other generic parameters @@ -658,14 +826,21 @@ if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN); } } - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - + // optionally, write out the covariance matrix values + // XXX do I need to pad this to match the biggest covar matrix? + if (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); @@ -697,6 +872,338 @@ } +bool pmSourcesRead_CMF_@CMFMODE@_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); + return false; + } + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (!source->modelFits) { + // XXX: where to find the number of models to expect? + source->modelFits = psArrayAllocEmpty(5); + } + psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); + if (!modelName) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); + psFree(row); + return false; + } + pmModelType modelType = pmModelClassGetType(modelName); + if (modelType < 0) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); + psFree(row); + return false; + } + pmModel *model = pmModelAlloc(modelType); + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); + + model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); + model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + + psEllipseAxes axes; + axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); + axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); + axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); + if (!pmPSF_AxesToModel(PAR, axes, modelType)) { + // Do we need to fail here or can this happen? + psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); + psFree(model); + psFree(row); + return false; + } + // XXX: clean this up + if (model->params->n > 7) { + PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); + } + // 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; + + psArrayAdd(source->modelFits, 1, model); + psFree(model); + + psFree(row); + } + + return true; +} + +// **** write out the radial flux values for the sources for a given matched-PSF image +// **** how do we distinguish the matched-PSF images from the non-matched version bool pmSourcesWrite_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) { + bool status = false; + psArray *table; + psMetadata *row; + psF32 xPos, yPos; + char keyword1[80], keyword2[80]; + + // 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; + } + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname); + + // we use this just to define the output vectors (which must be present for all objects) + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + psAssert (radMax, "this must have been defined and tested earlier!"); + psAssert (radMax->n, "this must have been defined and tested earlier!"); + psAssert (radMin->n == radMax->n, "inconsistent annular bins"); + + // 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, "min radius for SB profile", radMax->data.F32[i]); + } + + // the FWHM values are available if we measured a psf-matched convolved set + psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + + table = psArrayAllocEmpty (sources->n); + + // we write out all sources, regardless of quality. the source flags tell us the state + for (int i = 0; i < sources->n; i++) { + + // 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; + + // skip sources without radial aper measurements (or insufficient) + 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) + 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: + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); + psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psFree (radFlux); + psFree (radFluxErr); + psFree (radFluxStdev); + psFree (radFill); + + 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; + } + + 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; + } + psFree (outhead); + psFree (table); return true; } + +bool pmSourcesRead_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); + return false; + } + + long seq_first = -1; + long seq_last = -1; + psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); + long max_entries = -1; + long num_entries = -1; + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (seq_first == -1) { + seq_first = seq; + } + if (seq == seq_first) { + psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM"); + psVectorAppend(fwhmValues, value); + } + if (seq == seq_last) { + num_entries++; + } else { + num_entries = 1; + seq_last = seq; + } + if (num_entries > max_entries) { + max_entries = num_entries; + } + + if (!source->radialAper) { + // XXX: where to find the number of models to expect? + source->radialAper = psArrayAllocEmpty(5); + } + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); + + radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); + radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); + radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); + radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); + + psArrayAdd(source->radialAper, 1, radialAper); + + psFree(radialAper); + psFree(row); + } + + // check for consistency between the length of fwhmValues and the maximum number of entries for each row + if (fwhmValues->n != max_entries) { + psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n", + fwhmValues->n, max_entries); + psAssert(0, "fixme"); + } + + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); + psFree(fwhmValues); + + return true; +} Index: unk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c (revision 34258) +++ (revision ) @@ -1,618 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -// This version has elements intended for difference images & forced photometry: -// diffStats entries (good for dipoles); flux + flux error (for insignificant detections) -bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - 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); - } - } - - float magOffset; - float zeroptErr; - float fwhmMajor; - float fwhmMinor; - pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); - - table = psArrayAllocEmpty (sources->n); - - // we write out PSF-fits for all sources, regardless of quality. the source flags tell us the state - // by the time we call this function, all values should be assigned. let's use asserts to be sure in some cases. - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - pmSourceDiffStats diffStats; - pmSourceDiffStatsInit(&diffStats); - if (source->diffStats) { - diffStats = *source->diffStats; - } - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - psFitsWriteBlank(fits, header, extname); - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - assert (status); - - int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); - if (nPos) { - source->diffStats = pmSourceDiffStatsAlloc(); - source->diffStats->nGood = nPos; - source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); - source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); - source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); - source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); - } - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - - psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - assert (radialBinsLower->n == radialBinsUpper->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - -# if (0) - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - if (doPetrosian) { - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - } - } - - // Kron measurements - if (doKron) { - pmSourceKronValues *kron = source->extpars->kron; - if (kron) { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); - } - } - - // Isophot measurements - // XXX insert header data: isophotal level - if (doIsophotal) { - pmSourceIsophotalValues *isophot = source->extpars->isophot; - if (isophot) { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); - } - } - - // Flux Annuli - if (doAnnuli) { - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } - } else { - for (int j = 0; j < radialBinsLower->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // XXX if no model fits are saved, write out modelEXT? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_DV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: unk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c (revision 34258) +++ (revision ) @@ -1,644 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -// This version has elements intended for difference images & forced photometry: -// diffStats entries (good for dipoles); flux + flux error (for insignificant detections) -bool pmSourcesWrite_CMF_PS1_DV2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - 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++) { - pmSource *source = (pmSource *) sources->data[i]; - - // If source->seq is -1, source was generated in this analysis. If source->seq is - // not -1, source was read from elsewhere: in the latter case, preserve the source - // ID. source.seq is used instead of source.id since the latter is a const - // generated on Alloc, and would thus be wrong for read in sources. - if (source->seq == -1) { - source->seq = i; - } - - // set the 'best' values for various output fields: - pmSourceOutputs outputs; - pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset); - - pmSourceOutputsMoments moments; - pmSourceOutputsSetMoments (&moments, source); - - pmSourceDiffStats diffStats; - pmSourceDiffStatsInit(&diffStats); - if (source->diffStats) { - diffStats = *source->diffStats; - } - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); // XXX this is only measured for non-linear fits - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfFluxErr); - - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in real aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX", PS_DATA_F32, "instrumental flux in standard aperture", source->apFlux); - psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG", PS_DATA_F32, "aperture flux error", source->apFluxErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 1.0 R1)", moments.Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 4.0 R1)", moments.Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS", PS_DATA_S32, "nPos (n pix > 3 sigma)", diffStats.nGood); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO", PS_DATA_F32, "fPos / (fPos + fNeg)", diffStats.fRatio); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD", PS_DATA_F32, "nPos / (nPos + nNeg)", diffStats.nRatioBad); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)", diffStats.nRatioMask); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL", PS_DATA_F32, "nPos / (nGood + nMask + nBad)", diffStats.nRatioAll); - - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P", PS_DATA_F32, "distance to positive match source", diffStats.Rp); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P", PS_DATA_F32, "signal-to-noise of pos match src", diffStats.SNp); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M", PS_DATA_F32, "distance to negative match source", diffStats.Rm); - psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M", PS_DATA_F32, "signal-to-noise of neg match src", diffStats.SNm); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags (group 2)", source->mode2); - - // XXX not sure how to get this : need to load Nimages with weight? - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - psMetadata *header = psMetadataCopy(NULL, tableHeader); - pmSourceMasksHeader(header); - - if (table->n == 0) { - psFitsWriteBlank(fits, header, extname); - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(fits, header, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - - return true; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfFlux = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX"); - source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG"); - - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); - assert (status); - - int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS"); - if (nPos) { - source->diffStats = pmSourceDiffStatsAlloc(); - source->diffStats->nGood = nPos; - source->diffStats->fRatio = psMetadataLookupF32 (&status, row, "DIFF_FRATIO"); - source->diffStats->nRatioBad = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD"); - source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK"); - source->diffStats->nRatioAll = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL"); - - source->diffStats->Rp = psMetadataLookupF32 (&status, row, "DIFF_R_P"); - source->diffStats->SNp = psMetadataLookupF32 (&status, row, "DIFF_SN_P"); - source->diffStats->Rm = psMetadataLookupF32 (&status, row, "DIFF_R_M"); - source->diffStats->SNm = psMetadataLookupF32 (&status, row, "DIFF_SN_M"); - } - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - - psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - assert (radialBinsLower->n == radialBinsUpper->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? - - pmSource *source = sources->data[i]; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - -# if (0) - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - if (doPetrosian) { - pmSourcePetrosianValues *petrosian = source->extpars->petrosian; - if (petrosian) { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", petrosian->mag); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - } - } - - // Kron measurements - if (doKron) { - pmSourceKronValues *kron = source->extpars->kron; - if (kron) { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", kron->mag); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", kron->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); - } - } - - // Isophot measurements - // XXX insert header data: isophotal level - if (doIsophotal) { - pmSourceIsophotalValues *isophot = source->extpars->isophot; - if (isophot) { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", isophot->mag); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); - } - } - - // Flux Annuli - if (doAnnuli) { - pmSourceAnnuli *annuli = source->extpars->annuli; - if (annuli) { - psVector *fluxVal = annuli->flux; - psVector *fluxErr = annuli->fluxErr; - psVector *fluxVar = annuli->fluxVar; - - for (int j = 0; j < fluxVal->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } - } else { - for (int j = 0; j < radialBinsLower->n; j++) { - char name[32]; - sprintf (name, "FLUX_VAL_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); - sprintf (name, "FLUX_ERR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); - sprintf (name, "FLUX_VAR_R_%02d", j); - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); - } - } - } - -# endif - psArrayAdd (table, 100, row); - psFree (row); - } - - if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_DV2_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // XXX if no model fits are saved, write out modelEXT? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); - } - } - - // XXX other parameters which may be set. - // XXX flag / value to define the model - // XXX write out the model type, fit status flags - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (outhead); - psFree (table); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree (outhead); - psFree(table); - return false; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesWrite_CMF_PS1_DV2_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - return true; -} Index: unk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 34258) +++ (revision ) @@ -1,1116 +1,0 @@ -/** @file pmSourceIO.c - * - * @author EAM, IfA - * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 02:44:19 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "pmConfig.h" -#include "pmErrorCodes.h" -#include "pmDetrendDB.h" - -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" - -#include "pmTrend2D.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmSpan.h" -#include "pmFootprintSpans.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmModelFuncs.h" -#include "pmModel.h" -#include "pmModelUtils.h" -#include "pmModelClass.h" -#include "pmSourceMasks.h" -#include "pmSourceExtendedPars.h" -#include "pmSourceDiffStats.h" -#include "pmSource.h" -#include "pmSourceFitModel.h" -#include "pmPSF.h" -#include "pmPSFtry.h" - -#include "pmSourceIO.h" -#include "pmSourceOutputs.h" - -// panstarrs-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data - -// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes -// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes) - -bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe) -{ - 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. - // Careful: if we are working with child sources, then we need to sort by the parent info, - // not our info - if (sources->n > 0) { - pmSource *source = sources->data[0]; - if (source->seq == -1) { - 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 is unique and generated in this analysis. If - // source->seq is not -1, source was read from elsewhere or tied to other source (eg - // from another image): 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); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF", PS_DATA_F32, "PSF x coordinate", outputs.xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF", PS_DATA_F32, "PSF y coordinate", outputs.yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG", PS_DATA_F32, "Sigma in PSF x coordinate", outputs.xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG", PS_DATA_F32, "Sigma in PSF y coordinate", outputs.yErr); - psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE", PS_DATA_F32, "position angle at source (degrees)", outputs.posAngle); - psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE", PS_DATA_F32, "plate scale at source (arcsec/pixel)", outputs.pltScale); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG", PS_DATA_F32, "PSF fit instrumental magnitude", source->psfMag); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->psfMagErr); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX", PS_DATA_F32, "PSF fit instrumental flux (counts)", source->psfFlux); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux", source->psfFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG", PS_DATA_F32, "magnitude in standard aperture", source->apMag); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW", PS_DATA_F32, "magnitude in reported aperture", source->apMagRaw); - psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS", PS_DATA_F32, "radius used for aperture mags", outputs.apRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude", outputs.peakMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG", PS_DATA_F32, "PSF Magnitude using supplied calibration", outputs.calMag); - psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG", PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr); - psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF", PS_DATA_F64, "PSF RA coordinate (degrees)", outputs.ra); - psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF", PS_DATA_F64, "PSF DEC coordinate (degrees)", outputs.dec); - psMetadataAdd (row, PS_LIST_TAIL, "SKY", PS_DATA_F32, "Sky level", source->sky); - psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA", PS_DATA_F32, "Sigma of sky level", source->skyErr); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "Chisq of PSF-fit", outputs.chisq); - psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to CF", source->crNsigma); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA", PS_DATA_F32, "Nsigma deviations from PSF to EXT", source->extNsigma); - - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR", PS_DATA_F32, "PSF width (major axis)", outputs.psfMajor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR", PS_DATA_F32, "PSF width (minor axis)", outputs.psfMinor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA", PS_DATA_F32, "PSF orientation angle", outputs.psfTheta); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF", PS_DATA_F32, "PSF coverage/quality factor (bad)", source->pixWeightNotBad); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT", PS_DATA_F32, "PSF coverage/quality factor (poor)", source->pixWeightNotPoor); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF", PS_DATA_S32, "degrees of freedom", outputs.nDOF); - psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX", PS_DATA_S32, "number of pixels in fit", outputs.nPix); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX", PS_DATA_F32, "second moments (X^2)", moments.Mxx); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY", PS_DATA_F32, "second moments (X*Y)", moments.Mxy); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY", PS_DATA_F32, "second moments (Y*Y)", moments.Myy); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C", PS_DATA_F32, "third momemt cos theta", moments.M_c3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S", PS_DATA_F32, "third momemt sin theta", moments.M_s3); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C", PS_DATA_F32, "fourth momemt cos theta", moments.M_c4); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S", PS_DATA_F32, "fourth momemt sin theta", moments.M_s4); - - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1", PS_DATA_F32, "first radial moment", moments.Mrf); - psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH", PS_DATA_F32, "half radial moment", moments.Mrh); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Krf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR", PS_DATA_F32, "Kron Flux Error", moments.dKrf); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kinner); - psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER", PS_DATA_F32, "Kron Flux (in 2.5 R1)", moments.Kouter); - - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS", PS_DATA_U32, "psphot analysis flags", source->mode); - psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2", PS_DATA_U32, "psphot analysis flags", source->mode2); - - psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES", PS_DATA_U16, "Number of frames overlapping source center", source->nFrames); - psMetadataAdd (row, PS_LIST_TAIL, "PADDING", PS_DATA_S16, "padding", 0); - - psArrayAdd (table, 100, row); - psFree (row); - - // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not - // subtracted - - // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image - // edge; 3) any pixels in the 3x3 peak region are masked; - } - - // 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); - - if (table->n == 0) { - if (!psFitsWriteBlank(fits, header, extname)) { - psError(psErrorCodeLast(), false, "Unable to write blank sources file."); - psFree(table); - psFree(header); - return false; - } - psFree(table); - psFree(header); - return true; - } - - psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); - if (!psFitsWriteTable(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; -} - -// read in a readout from the fits file -psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(header, false); - - bool status; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - - // define PSF model type - // XXX need to carry the extra model parameters - int modelType = pmModelClassGetType ("PS_MODEL_GAUSS"); - - char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); - if (PSF_NAME != NULL) { - modelType = pmModelClassGetType (PSF_NAME); - } - assert (modelType > -1); - - // We get the size of the table, and allocate the array of sources first because the table - // is large and ephemeral --- when the table gets blown away, whatever is allocated after - // the table is read blocks the free. In fact, it's better to read the table row by row. - long numSources = psFitsTableSize(fits); // Number of sources in table - psArray *sources = psArrayAlloc(numSources); // Array of sources, to return - - // convert the table to the pmSource entriesa - for (int i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i); - psFree(sources); - return NULL; - } - - pmSource *source = pmSourceAlloc (); - pmModel *model = pmModelAlloc (modelType); - source->modelPSF = model; - source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags - - // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models. - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - source->seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - PAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG"); - axes.major = psMetadataLookupF32 (&status, row, "PSF_MAJOR"); - axes.minor = psMetadataLookupF32 (&status, row, "PSF_MINOR"); - axes.theta = psMetadataLookupF32 (&status, row, "PSF_THETA"); - - PAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY"); - dPAR[PM_PAR_SKY] = psMetadataLookupF32 (&status, row, "SKY_SIGMA"); - source->sky = PAR[PM_PAR_SKY]; - source->skyErr = dPAR[PM_PAR_SKY]; - - // XXX use these to determine PAR[PM_PAR_I0]? - source->psfMag = psMetadataLookupF32 (&status, row, "PSF_INST_MAG"); - source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG"); - source->apMag = psMetadataLookupF32 (&status, row, "AP_MAG"); - - // XXX this scaling is incorrect: does not include the 2 \pi AREA factor - PAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN; - dPAR[PM_PAR_I0] = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN; - - pmPSF_AxesToModel (PAR, axes, modelType); - - float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG"); - float peakFlux = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN; - - // recreate the peak to match (xPos, yPos) +/- (xErr, yErr) - source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->rawFlux = peakFlux; - source->peak->smoothFlux = peakFlux; - source->peak->xf = PAR[PM_PAR_XPOS]; // more accurate position - source->peak->yf = PAR[PM_PAR_YPOS]; // more accurate position - source->peak->dx = dPAR[PM_PAR_XPOS]; - source->peak->dy = dPAR[PM_PAR_YPOS]; - - source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF"); - source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT"); - source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); - source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); - source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); - - // note that some older versions used PSF_PROBABILITY: this was not well defined. - model->chisq = psMetadataLookupF32 (&status, row, "PSF_CHISQ"); - model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); - model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - - source->moments = pmMomentsAlloc (); - source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf - source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf - - source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX"); - source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY"); - source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY"); - - source->moments->Mrf = psMetadataLookupF32 (&status, row, "MOMENTS_R1"); - source->moments->Mrh = psMetadataLookupF32 (&status, row, "MOMENTS_RH"); - source->moments->KronFlux = psMetadataLookupF32 (&status, row, "KRON_FLUX"); - source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR"); - - source->moments->KronFinner = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER"); - source->moments->KronFouter = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER"); - - // 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 - source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C"); - source->moments->Mxxy = 0.0; - source->moments->Mxyy = 0.0; - source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S"); - - source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C"); - source->moments->Mxxxy = 0.0; - source->moments->Mxxyy = 0.0; - source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S"); - source->moments->Myyyy = 0.0; - - source->mode = psMetadataLookupU32 (&status, row, "FLAGS"); - source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2"); - assert (status); - - sources->data[i] = source; - psFree(row); - } - - return sources; -} - -bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - - bool status; - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psF32 xPos, yPos; - psF32 xErr, yErr; - int nRow = -1; - char keyword1[80], keyword2[80]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - pmChip *chip = readout->parent->parent; - pmFPA *fpa = chip->parent; - - // zero point corrections - bool status1 = false; - bool status2 = false; - float magOffset = 0.0; - float exptime = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE"); - float zeropt = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP"); - if (!isfinite(zeropt)) { - zeropt = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS"); - } - if (status1 && status2 && (exptime > 0.0)) { - magOffset = zeropt + 2.5*log10(exptime); - } - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // which extended source analyses should we perform? - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // 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, "min 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 - pmSource *thisSource = sources->data[i]; - - // this is the "real" version of this source - pmSource *source = thisSource->parent ? thisSource->parent : thisSource; - - // skip sources without measurements - if (source->extpars == NULL) continue; - - // we require a PSF model fit (ignore the real crud) - pmModel *model = source->modelPSF; - if (model == NULL) continue; - - // XXX I need to split the extended models from the extended aperture measurements - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - - 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); - - // Petrosian measurements - // XXX insert header data: petrosian ref radius, flux ratio - // XXX check flags to see if Pet was measured - if (doPetrosian) { - 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 - float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point - 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); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill); - } else { - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50", PS_DATA_F32, "Petrosian R50 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90", PS_DATA_F32, "Petrosian R90 (pix)", NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); - psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL", PS_DATA_F32, "Petrosian Fill Factor", NAN); - } - } - - // 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); - } - - 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; - } - - 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; - } - psFree (outhead); - psFree (table); - - return true; -} - -bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); - return false; - } - - // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point - float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); - float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); - float magOffset = zeropt + 2.5*log10(exptime); - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - pmSource *source = NULL; -#ifndef ASSUME_SORTED - long j = seq < sources->n ? seq : sources->n - 1; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } -#else - long j = sourceIndex[seq]; - psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); - source = sources->data[j]; -#endif - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - - if (!source->extpars) { - source->extpars = pmSourceExtendedParsAlloc (); - } - pmSourceExtendedPars *extpars = source->extpars; - - // Assume that X_EXT Y_EXT and sigmas match the psf src so skip - - // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by - // psf size or something? - extpars->axes.major = 1.0; - extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); - extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); - - 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)) { - extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; - } - } - - extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); - extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); - extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); - extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); - extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); - extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); - extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); - - psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); - psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); - psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); - - if (radSB && radSB->n > 0) { - extpars->radProfile = pmSourceRadialProfileAlloc(); - extpars->radProfile->binSB = psMemIncrRefCounter(radSB); - extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); - extpars->radProfile->binFill = psMemIncrRefCounter(radFill); - } - - psFree(row); - } - - return true; -} - -// XXX this layout is still the same as PS1_DEV_1 -bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) -{ - - psArray *table; - psMetadata *row; - psF32 *PAR, *dPAR; - psEllipseAxes axes; - psF32 xPos, yPos; - psF32 xErr, yErr; - char name[64]; - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams - int nParamMax = 0; - for (int i = 0; i < sources->n; i++) { - // 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->modelFits == NULL) continue; - for (int j = 0; j < source->modelFits->n; j++) { - pmModel *model = source->modelFits->data[j]; - assert (model); - nParamMax = PS_MAX (nParamMax, model->params->n); - } - } - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - pmSource *thisSource = sources->data[i]; - - // 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? - if (source->modelFits == NULL) continue; - - // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) - for (int j = 0; j < source->modelFits->n; j++) { - - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelFits->data[j]; - assert (model); - - // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0, model->type); - - row = psMetadataAlloc (); - - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); - psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); - - psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); - psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); - psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); - - // write out the other generic parameters - for (int k = 0; k < nParamMax; k++) { - if (k == PM_PAR_I0) continue; - if (k == PM_PAR_SKY) continue; - if (k == PM_PAR_XPOS) continue; - if (k == PM_PAR_YPOS) continue; - if (k == PM_PAR_SXX) continue; - if (k == PM_PAR_SXY) continue; - if (k == PM_PAR_SYY) continue; - - snprintf (name, 64, "EXT_PAR_%02d", k); - - if (k < model->params->n) { - psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); - } else { - psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN); - } - } - - // optionally, write out the covariance matrix values - // XXX do I need to pad this to match the biggest covar matrix? - if (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); - } - } - - 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; - } - - 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; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); - return false; - } - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq. - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - long j = seq < sources->n ? seq : sources->n - 1; - pmSource *source = NULL; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - if (!source->modelFits) { - // XXX: where to find the number of models to expect? - source->modelFits = psArrayAllocEmpty(5); - } - psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); - if (!modelName) { - psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); - psFree(row); - return false; - } - pmModelType modelType = pmModelClassGetType(modelName); - if (modelType < 0) { - psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); - psFree(row); - return false; - } - pmModel *model = pmModelAlloc(modelType); - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); - PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); - dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); - dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); - - model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); - model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); - - psEllipseAxes axes; - axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); - axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); - axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); - if (!pmPSF_AxesToModel(PAR, axes, modelType)) { - // Do we need to fail here or can this happen? - psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); - psFree(model); - psFree(row); - return false; - } - // XXX: clean this up - if (model->params->n > 7) { - PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); - } - // 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; - - psArrayAdd(source->modelFits, 1, model); - psFree(model); - - psFree(row); - } - - return true; -} - -// **** write out the radial flux values for the sources for a given matched-PSF image -// **** how do we distinguish the matched-PSF images from the non-matched version -bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe) -{ - bool status = false; - psArray *table; - psMetadata *row; - psF32 xPos, yPos; - char keyword1[80], keyword2[80]; - - // 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; - } - - // create a header to hold the output data - psMetadata *outhead = psMetadataAlloc (); - - // write the links to the image header - psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname); - - // we use this just to define the output vectors (which must be present for all objects) - psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); - psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); - psAssert (radMax, "this must have been defined and tested earlier!"); - psAssert (radMax->n, "this must have been defined and tested earlier!"); - psAssert (radMin->n == radMax->n, "inconsistent annular bins"); - - // 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, "min radius for SB profile", radMax->data.F32[i]); - } - - // the FWHM values are available if we measured a psf-matched convolved set - psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); - - // let's write these out in S/N order - sources = psArraySort (sources, pmSourceSortByFlux); - - table = psArrayAllocEmpty (sources->n); - - // we write out all sources, regardless of quality. the source flags tell us the state - for (int i = 0; i < sources->n; i++) { - - // 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; - - // skip sources without radial aper measurements (or insufficient) - 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) - 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: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); - psFree (radFlux); - psFree (radFluxErr); - psFree (radFluxStdev); - psFree (radFill); - - 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; - } - - 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; - } - psFree (outhead); - psFree (table); - return true; -} - -bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) -{ - PS_ASSERT_PTR_NON_NULL(fits, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - - bool status; - long numSources = psFitsTableSize(fits); // Number of sources in table - if (numSources == 0) { - psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); - return false; - } - - long seq_first = -1; - long seq_last = -1; - psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); - long max_entries = -1; - long num_entries = -1; - - for (long i = 0; i < numSources; i++) { - psMetadata *row = psFitsReadTableRow(fits, i); // Table row - if (!row) { - psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); - psFree(row); - return false; - } - // Find the source with this sequence number. - // XXX: I am assuming that sources is sorted in order of seq. - long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); - long j = seq < sources->n ? seq : sources->n - 1; - pmSource *source = NULL; - for (; j >= 0; j--) { - source = sources->data[j]; - if (source->seq == seq) { - break; - } - } - if (!source) { - psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); - psFree(row); - return false; - } - if (seq_first == -1) { - seq_first = seq; - } - if (seq == seq_first) { - psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM"); - psVectorAppend(fwhmValues, value); - } - if (seq == seq_last) { - num_entries++; - } else { - num_entries = 1; - seq_last = seq; - } - if (num_entries > max_entries) { - max_entries = num_entries; - } - - if (!source->radialAper) { - // XXX: where to find the number of models to expect? - source->radialAper = psArrayAllocEmpty(5); - } - pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); - - radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); - radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); - radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); - radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); - - psArrayAdd(source->radialAper, 1, radialAper); - - psFree(radialAper); - psFree(row); - } - - // check for consistency between the length of fwhmValues and the maximum number of entries for each row - if (fwhmValues->n != max_entries) { - psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n", - fwhmValues->n, max_entries); - psAssert(0, "fixme"); - } - - if (!readout->analysis) { - readout->analysis = psMetadataAlloc(); - } - - psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); - psFree(fwhmValues); - - return true; -}