IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31438


Ignore:
Timestamp:
May 4, 2011, 4:44:43 PM (15 years ago)
Author:
eugene
Message:

explicitly include the reference aperture in the growth curve; add some debugging output to pmGrowthCurveGenerate.c; use the reference aperture magnitude when using sources for the growth curve; extend the model integration to 150 sigmas2 (not just 50); re-generate modelFlux pixels if at all deviant from source

Location:
branches/eam_branches/ipp-20110404/psModules/src/objects
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_DEV.c

    r31153 r31438  
    268268    float f0 = 1.0;
    269269    float f1, f2;
    270     for (z = DZ; z < 50; z += DZ) {
     270    for (z = DZ; z < 150; z += DZ) {
    271271        f1 = exp(-pow(z,ALPHA));
    272272        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_EXP.c

    r31153 r31438  
    255255    float f0 = 1.0;
    256256    float f1, f2;
    257     for (z = DZ; z < 50; z += DZ) {
     257    for (z = DZ; z < 150; z += DZ) {
    258258        f1 = exp(-sqrt(z));
    259259        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_PGAUSS.c

    r31153 r31438  
    239239    float f0 = 1.0;
    240240    float f1, f2;
    241     for (z = DZ; z < 50; z += DZ) {
     241    for (z = DZ; z < 150; z += DZ) {
    242242        f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
    243243        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_PS1_V1.c

    r31153 r31438  
    261261    float f0 = 1.0;
    262262    float f1, f2;
    263     for (z = DZ; z < 50; z += DZ) {
     263    for (z = DZ; z < 150; z += DZ) {
    264264        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    265265        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_QGAUSS.c

    r31153 r31438  
    262262    float f0 = 1.0;
    263263    float f1, f2;
    264     for (z = DZ; z < 50; z += DZ) {
     264    for (z = DZ; z < 150; z += DZ) {
    265265        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
    266266        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_RGAUSS.c

    r31153 r31438  
    251251    float f0 = 1.0;
    252252    float f1, f2;
    253     for (z = DZ; z < 50; z += DZ) {
     253    for (z = DZ; z < 150; z += DZ) {
    254254        f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
    255255        z += DZ;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/models/pmModel_SERSIC.c

    r31153 r31438  
    316316    float f0 = 1.0;
    317317    float f1, f2;
    318     for (z = DZ; z < 50; z += DZ) {
     318    for (z = DZ; z < 150; z += DZ) {
    319319        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
    320320        f1 = exp(-pow(z,PAR[PM_PAR_7]));

  • branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c

    r31384 r31438  
    7373    // Fractional pixel radii are not well defined; use integer pixel radii.  Use 1 pixel steps
    7474    // until the scaling factor steps in intervals larger than 1 pixel
    75     float Rlin = 1.0 / (fR - 1.0);
     75    // float Rlin = 1.0 / (fR - 1.0);
    7676
    7777    growth->radius = psVectorAllocEmpty (NPTS, PS_DATA_F32);
    … …  
    7979    // there will be NPTS radii + a few extras
    8080    float radius = minRadius;
    81     while (radius < Rlin) {
     81    while (radius < refRadius) {
    8282        // fprintf (stderr, "r: %f\n", radius);
    8383        psVectorAppend (growth->radius, radius - 0.001);
    8484        radius += 1.0;
    8585    }   
     86    growth->refBin = growth->radius->n - 1;
    8687    while (radius < maxRadius) {
    8788        // fprintf (stderr, "r: %f\n", radius);
    … …  
    9798    growth->maxRadius = maxRadius;
    9899    growth->fitMag = NAN;
     100    growth->refMag = NAN;
    99101    growth->apRef  = NAN;
    100102    growth->apLoss = NAN;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.h

    r15562 r31438  
    2222    psF32 maxRadius;
    2323    psF32 fitMag;
     24    psF32 refMag;
    2425    psF32 apRef;   // apMag[refRadius]
    2526    psF32 apLoss;  // fitMag - apRef
     27    int refBin;
    2628}
    2729pmGrowthCurve;

  • branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c

    r31361 r31438  
    226226}
    227227
     228# define DEBUG 0
     229# if (DEBUG)
     230static FILE *fgr = NULL;
     231# endif
     232
    228233// we generate the growth curve for the center of the image with the specified psf model
    229234bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal)
    … …  
    240245    psArray *growths = psArrayAllocEmpty (100);
    241246
     247# if (DEBUG)
     248    fgr = fopen ("growth.mags.dat", "w");
     249# endif
     250
    242251    for (int i = 0; i < sources->n; i++) {
    243252
    … …  
    254263    psAssert (growths->n, "cannot build growth curve (no valid PSF stars?)");
    255264
    256 # if (0)
    257     FILE *f = fopen ("growth.dat", "w");
    258     assert (f);
    259 
    260     for (int j = 0; j < psf->growth->radius->n; j++) {
    261 
    262         fprintf (f, "%f : ", psf->growth->radius->data.F32[j]);
    263 
    264         // XXX dump the growth curves:
    265         for (int i = 0; i < growths->n; i++) {
    266            
    267             pmGrowthCurve *growth = growths->data[i];
    268             if (!growth) continue;
    269 
    270             fprintf (f, "%8.4f ", growth->apMag->data.F32[j]);
    271         }
    272         fprintf (f, "\n");
    273     }
    274     fclose (f);
     265# if (DEBUG)
     266    fclose (fgr);
    275267# endif
    276268
    … …  
    290282            pmGrowthCurve *growth = growths->data[j];
    291283            if (!isfinite(growth->apMag->data.F32[i])) continue;
    292             psVectorAppend (values, growth->apMag->data.F32[i] - growth->fitMag);
     284            psVectorAppend (values, growth->apMag->data.F32[i] - growth->refMag);
    293285        }
    294286        if (values->n == 0) {
    … …  
    354346        growth->apMag->data.F32[i] = source->apMag;
    355347    }
     348    psAssert(growth->refBin >= 0, "invalid growth reference bin");
     349    psAssert(growth->refBin < growth->apMag->n, "invalid growth reference bin");
     350    growth->refMag = growth->apMag->data.F32[growth->refBin];
     351
     352    // Loop over the range of radii
     353# if (DEBUG)
     354    for (int i = 0; i < growth->radius->n; i++) {
     355        fprintf (fgr, "%f %f  %f %f %f %f\n", xc, yc, growth->radius->data.F32[i], growth->apMag->data.F32[i], growth->fitMag, growth->refMag);
     356    }
     357# endif
     358
    356359    return growth;
    357360}
    358 
    359 # if (0)
    360     // solve for the best growth and the offsets per star to minimize the scatter
    361 
    362     // generate the per-star and per-radius sums
    363     int Nradius = psf->growth->radius->n;
    364     int Nstar = growths->n;
    365 
    366     psVector *Mstar   = psVectorAlloc (Nstar, PS_DATA_F32);
    367     psVector *Mradius = psVectorAlloc (Nradius, PS_DATA_F32);
    368 
    369     psVectorInit (Mstar, 0.0);
    370     psVectorInit (Mradius, 0.0);
    371 
    372     for (int i = 0; i < Nstar; i++) {
    373         pmGrowthCurve *growth = growths->data[i];
    374         for (int j = 0; j < Nradius; j++) {
    375             Mstar->data.F32[i] += growth->apMag->data.F32[j];
    376             Mradius->data.F32[j] += growth->apMag->data.F32[j];
    377         }
    378     }
    379 
    380     psImage *A = psImageAlloc(Nstar + Nradius, Nstar + Nradius, PS_DATA_F32);
    381     psVector *B = psVectorAlloc(Nstar + Nradius, PS_DATA_F32);
    382 
    383     psImageInit (A, 0.0);
    384     psVectorInit (B, 0.0);
    385 
    386     for (int i = 0; i < Nstar; i++) {
    387         A->data.F32[i][i] = Nradius;
    388         B->data.F32[i] = Mstar->data.F32[i];
    389     }
    390     for (int i = 0; i < Nradius; i++) {
    391         int j = i + Nstar;
    392         A->data.F32[j][j] = Nstar;
    393         B->data.F32[j] = Mradius->data.F32[i];
    394     }
    395    
    396     if (!psMatrixGJSolve(A, B)) {
    397         psAbort("GJ failure");
    398     }
    399 
    400     for (int i = 0; i < Nradius; i++) {
    401         psf->growth->apMag->data.F32[i] = B->data.F32[Nstar+i];
    402     }
    403 
    404     // XXX this analysis produces a biased growth curve: points with small radius are more
    405     // likely to err on the side of being too faint (bacause of mis-aligned aperture), while on
    406     // the bright side the tend to be too bright (because of contaminating neighbors).  Bright
    407     // stars are less likely to be biased,
    408 
    409     // XXX at this point, we could iterate and exclude some of the stars with low data quality
    410     // XXX we could also weight the above by flux or something
    411 
    412 # endif

  • branches/eam_branches/ipp-20110404/psModules/src/objects/pmPCMdata.c

    r31153 r31438  
    263263bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model) {
    264264
    265     bool newWindow = (source->pixels->numRows != pcm->modelFlux->numRows) || (source->pixels->numCols != pcm->modelFlux->numCols);
     265    bool sameWindow = (source->pixels->numRows == pcm->modelFlux->numRows);
     266    sameWindow     &= (source->pixels->numCols == pcm->modelFlux->numCols);
     267    sameWindow     &= (source->pixels->col0    == pcm->modelFlux->col0);
     268    sameWindow     &= (source->pixels->row0    == pcm->modelFlux->row0);
    266269
    267270    // re-count the number of unmasked pixels:
    268     if (newWindow) {
     271    if (!sameWindow) {
    269272        for (psS32 i = 0; i < source->pixels->numRows; i++) {
    270273            for (psS32 j = 0; j < source->pixels->numCols; j++) {
    … …  
    346349
    347350    // has the source pixel window changed?
    348     if (newWindow) {
     351    if (!sameWindow) {
    349352
    350353        // adjust all supporting images:

  • branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.c

    r31384 r31438  
    285285    extend |= (mySource->maskObj == NULL);
    286286    extend |= (mySource->maskView == NULL);
     287
     288    // if ((fabs(x - 2020) < 5) && (fabs(y - 366) < 5)) {
     289    //  if (extend) {
     290    //      fprintf (stderr, "extend T, %f, %f : %f, %f vs %f, %f : %f, %f\n",
     291    //               newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1,
     292    //               mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1);
     293    //  } else {
     294    //      fprintf (stderr, "extend F, %f, %f : %f, %f vs %f, %f : %f, %f\n",
     295    //               newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1,
     296    //               mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1);
     297    //  }
     298    // }
    287299
    288300    if (extend) {
Note: See TracChangeset for help on using the changeset viewer.