IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 19, 2011, 10:40:50 AM (15 years ago)
Author:
eugene
Message:

remove invalid references to isophot and kron from extended source analysis; use the measured sky stdev instead of the user guess; optionally allow stars as well as galaxies to get petrosians; fix target vs input seeing in psphotStack; allow radial apertures for single-image analysis as well as stack; calculate the radial aperture flux errors; calculate the petrosian parameter errors; turn on threads for ext source fits; save covar matrix for ext source fits; calculate petrosian fill factor

Location:
branches/eam_branches/ipp-20110213/psphot
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110213/psphot

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianStats.c

    r28013 r30707  
    66// generate the Petrosian radius and flux from the mean surface brightness (r_i)
    77
    8 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
     8float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
     9float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
     10float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1);
    911
    1012bool psphotPetrosianStats (pmSource *source) {
    … …  
    2527    psVector *binRad     = profile->radialBins;
    2628    psVector *area       = profile->area;
     29    psVector *binFill    = profile->binFill;
    2730
    2831    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    … …  
    3336    psVector *meanSB      = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3437    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     38    psVector *apixSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3539
    3640    float petRadius = NAN;
     41    float petRadiusErr = NAN;
    3742    float petFlux = NAN;
     43    float petFluxErr = NAN;
     44    float petArea = NAN;
     45    float petApix = NAN;
    3846
    3947    bool anyPetro = false;
    … …  
    4149    bool above = true;
    4250    float Asum = 0.0;
     51    float Psum = 0.0;
    4352    float Fsum = 0.0;
    4453    float dFsum2 = 0.0;
    … …  
    5665
    5766        float Area = area->data.F32[i];
    58         Asum += Area;
     67        Asum += Area;                   // Asum is the cumulative area interior to this bin
    5968        Fsum += binSB->data.F32[i] * Area;
    6069        dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);
     70        Psum += Area*binFill->data.F32[i]; // Psum is the cumulative number of pixels interior to this bin
    6171
    6272        float areaInner = 0.5 * Area;
    … …  
    8797
    8898        psVectorAppend(areaSum, Asum);
     99        psVectorAppend(apixSum, Psum);
    89100        psVectorAppend(fluxSum, Fsum);
    90101        psVectorAppend(fluxSumErr2, dFsum2);
    91102        psVectorAppend(refRadius, binRad->data.F32[i]);
    92103
    93         psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
     104        psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f  %5.1f\n",
    94105                 i, refRadius->data.F32[nOut],
    95106                 binSB->data.F32[i], binSBstdev->data.F32[i],
    96107                 meanSB->data.F32[nOut], meanSBerr,
    97108                 petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
    98                  fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
     109                 fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], apixSum->data.F32[nOut], areaInner);
    99110   
    100111        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    … …  
    104115            if (i == 0) {
    105116                // assume Fmax @ R = 0.0
    106                 petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
    107             } else {
    108                 petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     117                petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     118                petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
     119            } else {
     120                petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     121                petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
    109122            }
    110123            above = false;
    … …  
    128141    }
    129142
     143    // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
    130144    if (!anyPetro) {
    131145        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    132146        if (lowestSignificantRadius == 0) {
    133147            // assume Fmax @ R = 0.0
    134             petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     148            petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     149            petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
     150
    135151        } else {
    136             petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     152            int n0 = lowestSignificantRadius-1;
     153            int n1 = lowestSignificantRadius;
     154            petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
     155            petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]);
    137156        }
    138157    }
    … …  
    147166                continue;
    148167            } else {
    149                 petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     168                petFlux    = InterpolateValues     (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     169                petFluxErr = InterpolateValuesErrY (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
     170                petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
     171                petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
    150172                break;
    151173            }
    … …  
    158180    float R50 = NAN;
    159181    float R90 = NAN;
     182    float R50err = NAN;
     183    float R90err = NAN;
    160184    bool found50 = false;
    161185    bool found90 = false;
    … …  
    167191                continue;
    168192            } else {
    169                 R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     193                R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     194                R50err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
    170195                found50 = true;
    171196            }
    … …  
    176201                continue;
    177202            } else {
    178                 R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     203                R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     204                R90err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i]));
    179205                found90 = true;
    180206            }
    … …  
    188214    source->extpars->petrosianR50    = R50;
    189215    source->extpars->petrosianR90    = R90;
     216    source->extpars->petrosianFill      = petApix / petArea;
    190217   
    191218    // XXX add the errors
    192     source->extpars->petrosianRadiusErr = NAN;
    193     source->extpars->petrosianFluxErr   = NAN;
    194     source->extpars->petrosianR50Err    = NAN;
    195     source->extpars->petrosianR90Err    = NAN;
     219    source->extpars->petrosianRadiusErr = petRadiusErr;
     220    source->extpars->petrosianFluxErr   = petFluxErr;
     221    source->extpars->petrosianR50Err    = R50err;
     222    source->extpars->petrosianR90Err    = R90err;
    196223
    197224    // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
    … …  
    205232    psFree(meanSB);
    206233    psFree(areaSum);
     234    psFree(apixSum);
    207235
    208236    return true;
    … …  
    210238
    211239float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
    212     float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0);
     240    float dydx = (Y1 - Y0) / (X1 - X0);
     241    float Y = Y0 + dydx * (X - X0);
    213242    return Y;
    214243}
    215244
     245float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1) {
     246
     247    float dydx = (Y1 - Y0) / (X1 - X0);
     248    float dxdx = (X  - X0) / (X1 - X0);
     249   
     250    float dY = sqrt(PS_SQR(dX1*dydx*dxdx) + PS_SQR(dX0*dydx*(dxdx - 1.0)));
     251    return dY;
     252}
     253
     254float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1) {
     255
     256    float dxdx = (X  - X0) / (X1 - X0);
     257   
     258    float dY = sqrt(PS_SQR(dY1*dxdx) + PS_SQR(dY0*(1.0 - dxdx)));
     259    return dY;
     260}
     261
Note: See TracChangeset for help on using the changeset viewer.