IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 24, 2006, 2:07:03 PM (20 years ago)
Author:
Paul Price
Message:

Changed definition of psVectorAlloc and psArrayAlloc

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psphot/src/psphotApResid.c

    r9568 r9734  
    44
    55// measure the aperture residual statistics
    6 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 
     6bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
    77
    88    int Nfail = 0;
     
    1818    float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");
    1919
    20     // set limits on the aperture magnitudes 
     20    // set limits on the aperture magnitudes
    2121    pmSourceMagnitudesInit (recipe);
    2222
     
    2525    psf->growth = pmGrowthCurveAlloc (3.0, 100.0, REF_RADIUS);
    2626    psphotGrowthCurve (readout, psf);
    27    
    28     psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
    29     psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
    30     psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
    31     psVector *flux    = psVectorAlloc (300, PS_TYPE_F64);
    32     psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64);
    33     psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
    34     psVector *dMag    = psVectorAlloc (300, PS_TYPE_F64);
    35     mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0;
     27
     28    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
     29    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
     30    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
     31    psVector *flux    = psVectorAllocEmpty (300, PS_TYPE_F64);
     32    psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64);
     33    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
     34    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
    3635    Npsf = 0;
    3736
    3837    // select all good PM_SOURCE_TYPE_STAR entries
    3938    for (int i = 0; i < sources->n; i++) {
    40         source = sources->data[i];
    41         model = source->modelPSF;
    42 
    43         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    44         if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
    45         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    46         if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
    47         if (source->mode &  PM_SOURCE_MODE_POOR) continue;
    48 
    49         // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
    50         // will fail if below S/N threshold or model is missing
    51         if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {
    52             Nskip ++;
    53             continue;
    54         }
    55 
    56         apResid->data.F64[Npsf] = source->apMag - source->psfMag;
    57 
    58         xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    59         yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    60 
    61         flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
    62         r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];
    63 
    64         mask->data.U8[Npsf] = 0;
    65 
    66         // XXX sanity clip?
    67         // XXX need to see if all data were tossed?
    68         // XXX need to subtract median?
    69         // XXX need to put this in the config data...
    70         // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
    71         if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
    72             Nfail ++;
    73             continue;
    74         }
    75 
    76         dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    77 
    78         psVectorExtend (mask,    100, 1);
    79         psVectorExtend (xPos,    100, 1);
    80         psVectorExtend (yPos,    100, 1);
    81         psVectorExtend (flux,    100, 1);
    82         psVectorExtend (r2rflux, 100, 1);
    83         psVectorExtend (dMag,    100, 1);
    84         psVectorExtend (apResid, 100, 1);
    85         Npsf ++;
    86     }
    87     psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", 
    88               psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
     39        source = sources->data[i];
     40        model = source->modelPSF;
     41
     42        if (source->type != PM_SOURCE_TYPE_STAR) continue;
     43        if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
     44        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     45        if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
     46        if (source->mode &  PM_SOURCE_MODE_POOR) continue;
     47
     48        // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
     49        // will fail if below S/N threshold or model is missing
     50        if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {
     51            Nskip ++;
     52            continue;
     53        }
     54
     55        apResid->data.F64[Npsf] = source->apMag - source->psfMag;
     56
     57        xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
     58        yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
     59
     60        flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
     61        r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];
     62
     63        mask->data.U8[Npsf] = 0;
     64
     65        // XXX sanity clip?
     66        // XXX need to see if all data were tossed?
     67        // XXX need to subtract median?
     68        // XXX need to put this in the config data...
     69        // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
     70        if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
     71            Nfail ++;
     72            continue;
     73        }
     74
     75        dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
     76
     77        psVectorExtend (mask,    100, 1);
     78        psVectorExtend (xPos,    100, 1);
     79        psVectorExtend (yPos,    100, 1);
     80        psVectorExtend (flux,    100, 1);
     81        psVectorExtend (r2rflux, 100, 1);
     82        psVectorExtend (dMag,    100, 1);
     83        psVectorExtend (apResid, 100, 1);
     84        Npsf ++;
     85    }
     86    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n",
     87              psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
    8988
    9089    // XXX choose a better value here?
    9190    if (Npsf < 5) {
    92         psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up");
    93         return false;
    94     }
    95 
    96     // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 
     91        psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up");
     92        return false;
     93    }
     94
     95    // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
    9796    // APTREND options are used in the switch block below
    9897    pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
     
    10099    if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
    101100    if (ApTrendOption == PM_PSF_APTREND_ERROR) {
    102         psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
    103         return false;
     101        psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
     102        return false;
    104103    }
    105104
     
    112111    switch (ApTrendOption) {
    113112      case PM_PSF_APTREND_NONE:
    114         // remove ApTrend fit from pmPSFtry
    115         psf->ApTrend->coeff[0][0][0][0] = 0;
    116         break;
     113        // remove ApTrend fit from pmPSFtry
     114        psf->ApTrend->coeff[0][0][0][0] = 0;
     115        break;
    117116      case PM_PSF_APTREND_CONSTANT:
    118         stats->clipIter = 2;
    119         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    120         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    121         if (psf->ApTrend == NULL) {
    122             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    123             return false;
    124         }
    125         // apply the fit
    126         stats->clipIter = 3;
    127         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    128         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    129         if (psf->ApTrend == NULL) {
    130             psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
    131             return false;
    132         }
    133         break;
     117        stats->clipIter = 2;
     118        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     119        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     120        if (psf->ApTrend == NULL) {
     121            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     122            return false;
     123        }
     124        // apply the fit
     125        stats->clipIter = 3;
     126        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     127        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     128        if (psf->ApTrend == NULL) {
     129            psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
     130            return false;
     131        }
     132        break;
    134133      case PM_PSF_APTREND_SKYBIAS:
    135         // first clip out objects which are too far from the median
    136         stats->clipIter = 2;
    137         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    138         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    139         if (psf->ApTrend == NULL) {
    140             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    141             return false;
    142         }
    143         // apply the fit
    144         stats->clipIter = 3;
    145         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    146         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    147         if (psf->ApTrend == NULL) {
    148             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    149             return false;
    150         }
    151         break;
     134        // first clip out objects which are too far from the median
     135        stats->clipIter = 2;
     136        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     137        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     138        if (psf->ApTrend == NULL) {
     139            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     140            return false;
     141        }
     142        // apply the fit
     143        stats->clipIter = 3;
     144        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     145        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     146        if (psf->ApTrend == NULL) {
     147            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     148            return false;
     149        }
     150        break;
    152151      case PM_PSF_APTREND_SKYSAT:
    153         // first clip out objects which are too far from the median
    154         stats->clipIter = 2;
    155         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    156         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    157         if (psf->ApTrend == NULL) {
    158             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    159             return false;
    160         }
    161         // apply the fit
    162         stats->clipIter = 2;
    163         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    164         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    165         if (psf->ApTrend == NULL) {
    166             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    167             return false;
    168         }
    169         // apply the fit
    170         stats->clipIter = 3;
    171         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
    172         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    173         if (psf->ApTrend == NULL) {
    174             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
    175             return false;
    176         }
    177         break;
     152        // first clip out objects which are too far from the median
     153        stats->clipIter = 2;
     154        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     155        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     156        if (psf->ApTrend == NULL) {
     157            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     158            return false;
     159        }
     160        // apply the fit
     161        stats->clipIter = 2;
     162        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     163        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     164        if (psf->ApTrend == NULL) {
     165            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     166            return false;
     167        }
     168        // apply the fit
     169        stats->clipIter = 3;
     170        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
     171        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     172        if (psf->ApTrend == NULL) {
     173            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
     174            return false;
     175        }
     176        break;
    178177      case PM_PSF_APTREND_XY_LIN:
    179         // first clip out objects which are too far from the median
    180         stats->clipIter = 2;
    181         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    182         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    183         if (psf->ApTrend == NULL) {
    184             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    185             return false;
    186         }
    187         // apply the fit
    188         stats->clipIter = 3;
    189         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
    190         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    191         if (psf->ApTrend == NULL) {
    192             psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
    193             return false;
    194         }
    195         break;
     178        // first clip out objects which are too far from the median
     179        stats->clipIter = 2;
     180        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     181        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     182        if (psf->ApTrend == NULL) {
     183            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     184            return false;
     185        }
     186        // apply the fit
     187        stats->clipIter = 3;
     188        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
     189        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     190        if (psf->ApTrend == NULL) {
     191            psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
     192            return false;
     193        }
     194        break;
    196195      case PM_PSF_APTREND_XY_QUAD:
    197         // first clip out objects which are too far from the median
    198         stats->clipIter = 2;
    199         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    200         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    201         if (psf->ApTrend == NULL) {
    202             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    203             return false;
    204         }
    205         // apply the fit
    206         stats->clipIter = 3;
    207         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
    208         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    209         if (psf->ApTrend == NULL) {
    210             psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
    211             return false;
    212         }
    213         break;
     196        // first clip out objects which are too far from the median
     197        stats->clipIter = 2;
     198        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     199        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     200        if (psf->ApTrend == NULL) {
     201            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     202            return false;
     203        }
     204        // apply the fit
     205        stats->clipIter = 3;
     206        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
     207        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     208        if (psf->ApTrend == NULL) {
     209            psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
     210            return false;
     211        }
     212        break;
    214213      case PM_PSF_APTREND_SKY_XY_LIN:
    215         // first clip out objects which are too far from the median
    216         stats->clipIter = 2;
    217         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    218         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    219         if (psf->ApTrend == NULL) {
    220             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    221             return false;
    222         }
    223         // apply the fit
    224         stats->clipIter = 3;
    225         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
    226         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    227         if (psf->ApTrend == NULL) {
    228             psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
    229             return false;
    230         }
    231         break;
     214        // first clip out objects which are too far from the median
     215        stats->clipIter = 2;
     216        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     217        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     218        if (psf->ApTrend == NULL) {
     219            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     220            return false;
     221        }
     222        // apply the fit
     223        stats->clipIter = 3;
     224        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
     225        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     226        if (psf->ApTrend == NULL) {
     227            psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
     228            return false;
     229        }
     230        break;
    232231      case PM_PSF_APTREND_SKYSAT_XY_LIN:
    233         // first clip out objects which are too far from the median
    234         stats->clipIter = 2;
    235         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    236         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    237         if (psf->ApTrend == NULL) {
    238             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    239             return false;
    240         }
    241         // apply the fit
    242         stats->clipIter = 3;
    243         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    244         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    245         if (psf->ApTrend == NULL) {
    246             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    247             return false;
    248         }
    249         // apply the fit
    250         stats->clipIter = 3;
    251         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
    252         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    253         if (psf->ApTrend == NULL) {
    254             psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
    255             return false;
    256         }
    257         break;
     232        // first clip out objects which are too far from the median
     233        stats->clipIter = 2;
     234        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     235        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     236        if (psf->ApTrend == NULL) {
     237            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     238            return false;
     239        }
     240        // apply the fit
     241        stats->clipIter = 3;
     242        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     243        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     244        if (psf->ApTrend == NULL) {
     245            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     246            return false;
     247        }
     248        // apply the fit
     249        stats->clipIter = 3;
     250        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
     251        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     252        if (psf->ApTrend == NULL) {
     253            psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
     254            return false;
     255        }
     256        break;
    258257      case PM_PSF_APTREND_ALL:
    259         // first clip out objects which are too far from the median
    260         stats->clipIter = 2;
    261         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    262         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    263         if (psf->ApTrend == NULL) {
    264             psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
    265             return false;
    266         }
    267         // fit just SkyBias and clip out objects which are too far from the median
    268         stats->clipIter = 2;
    269         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    270         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    271         if (psf->ApTrend == NULL) {
    272             psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
    273             return false;
    274         }
    275         // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    276         stats->clipIter = 3;
    277         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
    278         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    279         if (psf->ApTrend == NULL) {
    280             psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
    281             return false;
    282         }
    283         break;
     258        // first clip out objects which are too far from the median
     259        stats->clipIter = 2;
     260        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     261        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     262        if (psf->ApTrend == NULL) {
     263            psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
     264            return false;
     265        }
     266        // fit just SkyBias and clip out objects which are too far from the median
     267        stats->clipIter = 2;
     268        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     269        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     270        if (psf->ApTrend == NULL) {
     271            psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
     272            return false;
     273        }
     274        // finally, fit x, y, SkyBias and clip out objects which are too far from the median
     275        stats->clipIter = 3;
     276        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
     277        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     278        if (psf->ApTrend == NULL) {
     279            psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
     280            return false;
     281        }
     282        break;
    284283      default:
    285         psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
    286         return false;
     284        psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
     285        return false;
    287286    }
    288287
     
    295294    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
    296295    for (int i = 0; i < dMag->n; i++) {
    297         if (dMag->data.F64[i] > 0.01) {
    298             mask->data.U8[i] |= 0x02;
    299         }
    300         if (! mask->data.U8[i]) Nkeep ++;
     296        if (dMag->data.F64[i] > 0.01) {
     297            mask->data.U8[i] |= 0x02;
     298        }
     299        if (! mask->data.U8[i]) Nkeep ++;
    301300    }
    302301    residStats  = psVectorStats (residStats, resid, NULL, mask, 0x03);
     
    320319
    321320    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
    322     psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", 
    323               psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);
    324     psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", 
    325               1e3*psf->ApTrend->coeff[1][0][0][0],
    326               1e6*psf->ApTrend->coeff[2][0][0][0],
    327               1e6*psf->ApTrend->coeff[1][1][0][0],
    328               1e3*psf->ApTrend->coeff[0][1][0][0],
    329               1e6*psf->ApTrend->coeff[0][2][0][0]);
     321    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n",
     322              psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);
     323    psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n",
     324              1e3*psf->ApTrend->coeff[1][0][0][0],
     325              1e6*psf->ApTrend->coeff[2][0][0][0],
     326              1e6*psf->ApTrend->coeff[1][1][0][0],
     327              1e3*psf->ApTrend->coeff[0][1][0][0],
     328              1e6*psf->ApTrend->coeff[0][2][0][0]);
    330329
    331330    psFree (mask);
Note: See TracChangeset for help on using the changeset viewer.