Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 14652) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 14937) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-24 00:11:02 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-09-21 00:06:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,4 +24,5 @@ #include "pmResiduals.h" #include "pmGrowthCurve.h" +#include "pmTrend2D.h" #include "pmPSF.h" #include "pmModel.h" @@ -70,7 +71,7 @@ test->psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); - test->metric = psVectorAlloc (sources->n, PS_TYPE_F64); - test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F64); - test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F64); + test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); + test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); + test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F32); test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); @@ -79,7 +80,7 @@ test->sources->data[i] = pmSourceCopy (sources->data[i]); test->mask->data.U8[i] = 0; - test->metric->data.F64[i] = 0; - test->metricErr->data.F64[i] = 0; - test->fitMag->data.F64[i] = 0; + test->metric->data.F32[i] = 0; + test->metricErr->data.F32[i] = 0; + test->fitMag->data.F32[i] = 0; } @@ -183,7 +184,7 @@ } - psfTry->fitMag->data.F64[i] = source->psfMag; - psfTry->metric->data.F64[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F64[i] = source->errMag; + psfTry->fitMag->data.F32[i] = source->psfMag; + psfTry->metric->data.F32[i] = source->apMag - source->psfMag; + psfTry->metricErr->data.F32[i] = source->errMag; Npsf ++; } @@ -195,6 +196,6 @@ // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) // this should be linear for Poisson errors and quadratic for constant sky errors - psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK); @@ -203,10 +204,10 @@ pmSource *source = psfTry->sources->data[i]; if (source->modelPSF == NULL) { - flux->data.F64[i] = 0.0; - chisq->data.F64[i] = 0.0; + flux->data.F32[i] = 0.0; + chisq->data.F32[i] = 0.0; mask->data.U8[i] = 0xff; } else { - flux->data.F64[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; - chisq->data.F64[i] = source->modelPSF->chisq / source->modelPSF->nDOF; + flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; + chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF; mask->data.U8[i] = 0; } @@ -259,14 +260,14 @@ // r2rflux = radius^2 * ten(0.4*fitMag); - psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); for (int i = 0; i < psfTry->sources->n; i++) { if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); - } - - // XXX this analysis of the apResid statistics is only approximate. The fitted magnitudes - // measure at this point (in the PSF fit) use Poisson errors, and are thus biased as a + r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]); + } + + // This analysis of the apResid statistics is only approximate. The fitted magnitudes + // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a // function of magnitude. We re-measure the apResid statistics later in psphot using the // linear, constant-error fitting. Do not reject outliers with excessive vigor here. @@ -277,9 +278,4 @@ stats->clipIter = 3; - // XXX we used to include an asymmetric clipping in order to toss out contaminated stars. - // test this on the very crowded field data. - // stats->min = 1.0; - // stats->max = 3.0; - // fit ApTrend only to r2rflux, ignore x,y,flux variations for now // linear clipped fit of ApResid to r2rflux @@ -287,4 +283,5 @@ poly->mask[1] = 1; // fit only a constant offset (no SKYBIAS) + // XXX replace this with a psVectorStats call? since we are not fitting the trend bool result = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, psfTry->metricErr, r2rflux); if (!result) { @@ -312,9 +309,9 @@ fprintf (f, "%d %d, %f %f %f %f %f %f %f\n", i, keep, x, y, - psfTry->fitMag->data.F64[i], - r2rflux->data.F64[i], - psfTry->metric->data.F64[i], - psfTry->metricErr->data.F64[i], - apfit->data.F64[i]); + psfTry->fitMag->data.F32[i], + r2rflux->data.F32[i], + psfTry->metric->data.F32[i], + psfTry->metricErr->data.F32[i], + apfit->data.F32[i]); } fclose (f); @@ -322,8 +319,5 @@ } - pmPSFMaskApTrend (psfTry->psf->ApTrend, PM_PSF_APTREND_SKYBIAS); - psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0]; - psfTry->psf->ApTrend->coeff[0][0][1][0] = 0; - + // XXX drop the skyBias value, or include above?? psfTry->psf->skyBias = poly->coeff[1]; psfTry->psf->ApResid = poly->coeff[0]; @@ -356,11 +350,11 @@ // construct the fit vectors from the collection of objects - psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); psVector *dz = NULL; if (applyWeight) { - dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); } @@ -371,11 +365,10 @@ continue; - // use F64 for polynomial fitting - x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; - y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; + x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; // weight by the error on the source flux if (dz) { - dz->data.F64[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0]; + dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0]; } } @@ -397,7 +390,7 @@ // convert the measured source shape paramters to polarization terms - psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - psVector *e1 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); - psVector *e2 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *e1 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); + psVector *e2 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); for (int i = 0; i < psfTry->sources->n; i++) { pmSource *source = psfTry->sources->data[i]; @@ -406,7 +399,7 @@ psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); - e0->data.F64[i] = pol.e0; - e1->data.F64[i] = pol.e1; - e2->data.F64[i] = pol.e2; + e0->data.F32[i] = pol.e0; + e1->data.F32[i] = pol.e1; + e2->data.F32[i] = pol.e2; } @@ -427,6 +420,6 @@ for (int i = 0; i < e0->n; i++) { fprintf (f, "%f %f : %f %f %f : %d\n", - x->data.F64[i], y->data.F64[i], - e0->data.F64[i], e1->data.F64[i], e2->data.F64[i], psfTry->mask->data.U8[i]); + x->data.F32[i], y->data.F32[i], + e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]); } fclose (f); @@ -457,5 +450,5 @@ if (source->modelEXT == NULL) continue; - z->data.F64[j] = source->modelEXT->params->data.F32[i]; + z->data.F32[j] = source->modelEXT->params->data.F32[i]; }