IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 7, 2007, 10:33:31 AM (19 years ago)
Author:
magnier
Message:

convert to pmTrend2D for ApTrend; convert exiting vectors to psF32 (the mix of psF64 entries and psF32 entries led to confusion)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20070830/psModules/src/objects/pmPSFtry.c

    r14652 r14783  
    55 *  @author EAM, IfA
    66 *
    7  *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-08-24 00:11:02 $
     7 *  @version $Revision: 1.44.2.1 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2007-09-07 20:33:31 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "pmResiduals.h"
    2525#include "pmGrowthCurve.h"
     26#include "pmTrend2D.h"
    2627#include "pmPSF.h"
    2728#include "pmModel.h"
    … …  
    7071
    7172    test->psf       = pmPSFAlloc (type, poissonErrors, psfTrendMask);
    72     test->metric    = psVectorAlloc (sources->n, PS_TYPE_F64);
    73     test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F64);
    74     test->fitMag    = psVectorAlloc (sources->n, PS_TYPE_F64);
     73    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
     74    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
     75    test->fitMag    = psVectorAlloc (sources->n, PS_TYPE_F32);
    7576    test->mask      = psVectorAlloc (sources->n, PS_TYPE_U8);
    7677
    … …  
    7980        test->sources->data[i] = pmSourceCopy (sources->data[i]);
    8081        test->mask->data.U8[i]  = 0;
    81         test->metric->data.F64[i] = 0;
    82         test->metricErr->data.F64[i] = 0;
    83         test->fitMag->data.F64[i] = 0;
     82        test->metric->data.F32[i] = 0;
     83        test->metricErr->data.F32[i] = 0;
     84        test->fitMag->data.F32[i] = 0;
    8485    }
    8586
    … …  
    183184        }
    184185
    185         psfTry->fitMag->data.F64[i] = source->psfMag;
    186         psfTry->metric->data.F64[i] = source->apMag - source->psfMag;
    187         psfTry->metricErr->data.F64[i] = source->errMag;
     186        psfTry->fitMag->data.F32[i] = source->psfMag;
     187        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
     188        psfTry->metricErr->data.F32[i] = source->errMag;
    188189        Npsf ++;
    189190    }
    … …  
    195196    // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    196197    // this should be linear for Poisson errors and quadratic for constant sky errors
    197     psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    198     psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     198    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     199    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    199200    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK);
    200201
    … …  
    203204        pmSource *source = psfTry->sources->data[i];
    204205        if (source->modelPSF == NULL) {
    205             flux->data.F64[i] = 0.0;
    206             chisq->data.F64[i] = 0.0;
     206            flux->data.F32[i] = 0.0;
     207            chisq->data.F32[i] = 0.0;
    207208            mask->data.U8[i] = 0xff;
    208209        } else {
    209             flux->data.F64[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
    210             chisq->data.F64[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
     210            flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
     211            chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
    211212            mask->data.U8[i] = 0;
    212213        }
    … …  
    259260
    260261    // r2rflux = radius^2 * ten(0.4*fitMag);
    261     psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     262    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    262263
    263264    for (int i = 0; i < psfTry->sources->n; i++) {
    264265        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
    265266            continue;
    266         r2rflux->data.F64[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
    267     }
    268 
    269     // XXX this analysis of the apResid statistics is only approximate.  The fitted magnitudes
    270     // measure at this point (in the PSF fit) use Poisson errors, and are thus biased as a
     267        r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
     268    }
     269
     270    // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
     271    // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
    271272    // function of magnitude.  We re-measure the apResid statistics later in psphot using the
    272273    // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
    … …  
    277278    stats->clipIter = 3;
    278279
    279     // XXX we used to include an asymmetric clipping in order to toss out contaminated stars.
    280     // test this on the very crowded field data.
    281     // stats->min = 1.0;
    282     // stats->max = 3.0;
    283 
    284280    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    285281    // linear clipped fit of ApResid to r2rflux
    … …  
    287283    poly->mask[1] = 1; // fit only a constant offset (no SKYBIAS)
    288284
     285    // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    289286    bool result = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, psfTry->metricErr, r2rflux);
    290287    if (!result) {
    … …  
    312309            fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
    313310                     i, keep, x, y,
    314                      psfTry->fitMag->data.F64[i],
    315                      r2rflux->data.F64[i],
    316                      psfTry->metric->data.F64[i],
    317                      psfTry->metricErr->data.F64[i],
    318                      apfit->data.F64[i]);
     311                     psfTry->fitMag->data.F32[i],
     312                     r2rflux->data.F32[i],
     313                     psfTry->metric->data.F32[i],
     314                     psfTry->metricErr->data.F32[i],
     315                     apfit->data.F32[i]);
    319316        }
    320317        fclose (f);
    … …  
    322319    }
    323320
    324     pmPSFMaskApTrend (psfTry->psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    325     psfTry->psf->ApTrend->coeff[0][0][0][0] = poly->coeff[0];
    326     psfTry->psf->ApTrend->coeff[0][0][1][0] = 0;
    327 
     321    // XXX drop the skyBias value, or include above??
    328322    psfTry->psf->skyBias  = poly->coeff[1];
    329323    psfTry->psf->ApResid  = poly->coeff[0];
    … …  
    356350
    357351    // construct the fit vectors from the collection of objects
    358     psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    359     psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    360     psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     352    psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     353    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     354    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    361355
    362356    psVector *dz = NULL;
    363357    if (applyWeight) {
    364         dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     358        dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    365359    }
    366360
    … …  
    371365            continue;
    372366
    373         // use F64 for polynomial fitting
    374         x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    375         y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
     367        x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
     368        y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
    376369
    377370        // weight by the error on the source flux
    378371        if (dz) {
    379             dz->data.F64[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0];
     372            dz->data.F32[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0];
    380373        }
    381374    }
    … …  
    397390
    398391    // convert the measured source shape paramters to polarization terms
    399     psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    400     psVector *e1   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
    401     psVector *e2   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
     392    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     393    psVector *e1   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     394    psVector *e2   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    402395    for (int i = 0; i < psfTry->sources->n; i++) {
    403396        pmSource *source = psfTry->sources->data[i];
    … …  
    406399        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
    407400
    408         e0->data.F64[i] = pol.e0;
    409         e1->data.F64[i] = pol.e1;
    410         e2->data.F64[i] = pol.e2;
     401        e0->data.F32[i] = pol.e0;
     402        e1->data.F32[i] = pol.e1;
     403        e2->data.F32[i] = pol.e2;
    411404    }
    412405
    … …  
    427420        for (int i = 0; i < e0->n; i++) {
    428421            fprintf (f, "%f %f  :  %f %f %f  : %d\n",
    429                      x->data.F64[i], y->data.F64[i],
    430                      e0->data.F64[i], e1->data.F64[i], e2->data.F64[i], psfTry->mask->data.U8[i]);
     422                     x->data.F32[i], y->data.F32[i],
     423                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]);
    431424        }
    432425        fclose (f);
    … …  
    457450            if (source->modelEXT == NULL)
    458451                continue;
    459             z->data.F64[j] = source->modelEXT->params->data.F32[i];
     452            z->data.F32[j] = source->modelEXT->params->data.F32[i];
    460453        }
    461454
Note: See TracChangeset for help on using the changeset viewer.