IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 14, 2008, 9:07:09 AM (18 years ago)
Author:
eugene
Message:

return null for no stars

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroLuminosityFunction.c

    r17556 r17677  
    55static void pmLumFuncFree (pmLumFunc *func) {
    66
    7   if (func == NULL) return;
     7    if (func == NULL) return;
    88
    9   return;
     9    return;
    1010}
    1111
    1212pmLumFunc *pmLumFuncAlloc () {
    1313
    14   pmLumFunc *func = (pmLumFunc *) psAlloc(sizeof(pmLumFunc));
    15   psMemSetDeallocator(func, (psFreeFunc) pmLumFuncFree);
     14    pmLumFunc *func = (pmLumFunc *) psAlloc(sizeof(pmLumFunc));
     15    psMemSetDeallocator(func, (psFreeFunc) pmLumFuncFree);
    1616
    17   func->mMin = 0.0;
    18   func->mMax = 0.0;
    19   func->slope = 0.0;
    20   func->offset = 0.0;
     17    func->mMin = 0.0;
     18    func->mMax = 0.0;
     19    func->slope = 0.0;
     20    func->offset = 0.0;
    2121
    22   func->mPeak = 0.0;
    23   func->nPeak = 0;
    24   func->sPeak = 0;
     22    func->mPeak = 0.0;
     23    func->nPeak = 0;
     24    func->sPeak = 0;
    2525
    26   return func;
     26    return func;
    2727}
    2828
    2929pmLumFunc *psastroLuminosityFunction (psArray *stars) {
    3030
    31   // determine the max and min magnitude for the array of stars
    32   pmAstromObj *star = stars->data[0];
    33   double mMin = star->Mag;
    34   double mMax = star->Mag;
    35   for (int i = 0; i < stars->n; i++) {
    36     star = stars->data[i];
    37     mMin = PS_MIN (star->Mag, mMin);
    38     mMax = PS_MAX (star->Mag, mMax);
    39   }
     31    if (stars->n == 0) return NULL;
    4032
    41   int nBin = 1 + (mMax - mMin) / dMag;
    42   if (nBin <= 1) {
    43     psLogMsg ("psastro", 4, "insufficient valid stars for this readout\n");
    44     return NULL;
    45   }
     33    // determine the max and min magnitude for the array of stars
     34    pmAstromObj *star = stars->data[0];
     35    double mMin = star->Mag;
     36    double mMax = star->Mag;
     37    for (int i = 0; i < stars->n; i++) {
     38        star = stars->data[i];
     39        mMin = PS_MIN (star->Mag, mMin);
     40        mMax = PS_MAX (star->Mag, mMax);
     41    }
    4642
    47   // create a histogram of the magnitudes
    48   // bin[0] = mMin : mMin + dMag
    49   // bin[i] = mMin + i*dMag : mMin + (i+1)*dMag
    50   psVector *nMags = psVectorAlloc (nBin, PS_TYPE_F32);
    51   psVectorInit (nMags, 0);
     43    int nBin = 1 + (mMax - mMin) / dMag;
     44    if (nBin <= 1) {
     45        psLogMsg ("psastro", 4, "insufficient valid stars for this readout\n");
     46        return NULL;
     47    }
    5248
    53   for (int i = 0; i < stars->n; i++) {
    54     star = stars->data[i];
    55     if (!isfinite(star->Mag)) continue;
    56     int bin = (star->Mag - mMin) / dMag;
    57     if (bin < 0) psAbort("bin cannot be negative!");
    58     if (bin >= nBin) psAbort("bin cannot be > %d!", nBin);
    59     nMags->data.F32[bin] += 1.0;
    60   }
     49    // create a histogram of the magnitudes
     50    // bin[0] = mMin : mMin + dMag
     51    // bin[i] = mMin + i*dMag : mMin + (i+1)*dMag
     52    psVector *nMags = psVectorAlloc (nBin, PS_TYPE_F32);
     53    psVectorInit (nMags, 0);
    6154
    62   // find the peak and position
    63   int iPeak = 0;
    64   int nPeak = 0;
    65   int sPeak = 0;
    66   for (int i = 0; i < nMags->n; i++) {
    67     if (nMags->data.F32[i] < nPeak) continue;
    68     iPeak = i;
    69     nPeak = nMags->data.F32[i];
    70     sPeak += nPeak;
    71   }
     55    for (int i = 0; i < stars->n; i++) {
     56        star = stars->data[i];
     57        if (!isfinite(star->Mag)) continue;
     58        int bin = (star->Mag - mMin) / dMag;
     59        if (bin < 0) psAbort("bin cannot be negative!");
     60        if (bin >= nBin) psAbort("bin cannot be > %d!", nBin);
     61        nMags->data.F32[bin] += 1.0;
     62    }
    7263
    73   psVector *lnMag = psVectorAllocEmpty (nBin, PS_TYPE_F32);
    74   psVector *Mag = psVectorAllocEmpty (nBin, PS_TYPE_F32);
     64    // find the peak and position
     65    int iPeak = 0;
     66    int nPeak = 0;
     67    int sPeak = 0;
     68    for (int i = 0; i < nMags->n; i++) {
     69        if (nMags->data.F32[i] < nPeak) continue;
     70        iPeak = i;
     71        nPeak = nMags->data.F32[i];
     72        sPeak += nPeak;
     73    }
    7574
    76   // create 2 vectors represnting the dlogN/dlogS line
    77   // exclude bins with no stars
    78   // exclude points after the peak with N < 0.8*peak
    79   int n = 0;
    80   for (int i = 0; i < nMags->n; i++) {
    81     if (nMags->data.F32[i] < 1) continue;
    82     if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue;
    83     lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
    84     Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
    85     n++;
    86   }
    87   lnMag->n = n;
    88   Mag->n = n;
    89   psLogMsg ("psastro", 4, "fitting %d points to luminosity function\n", n);
     75    psVector *lnMag = psVectorAllocEmpty (nBin, PS_TYPE_F32);
     76    psVector *Mag = psVectorAllocEmpty (nBin, PS_TYPE_F32);
    9077
    91   psVector *mask = psVectorAlloc (Mag->n, PS_TYPE_MASK);
    92   psVectorInit (mask, 0);
     78    // create 2 vectors represnting the dlogN/dlogS line
     79    // exclude bins with no stars
     80    // exclude points after the peak with N < 0.8*peak
     81    int n = 0;
     82    for (int i = 0; i < nMags->n; i++) {
     83        if (nMags->data.F32[i] < 1) continue;
     84        if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue;
     85        lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
     86        Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
     87        n++;
     88    }
     89    lnMag->n = n;
     90    Mag->n = n;
     91    psLogMsg ("psastro", 4, "fitting %d points to luminosity function\n", n);
    9392
    94   psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    95   psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    96   stats->clipSigma = 3.0;
    97   stats->clipIter = 3;
     93    psVector *mask = psVectorAlloc (Mag->n, PS_TYPE_MASK);
     94    psVectorInit (mask, 0);
    9895
    99   if (!psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag)) {
    100       psLogMsg ("psastro", 4, "Failed to fit the luminosity function\n");
    101       return(NULL);
    102   }
     96    psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     97    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     98    stats->clipSigma = 3.0;
     99    stats->clipIter = 3;
    103100
    104   // find min and max unmasked magnitudes
    105   double mMinValid = NAN;
    106   double mMaxValid = NAN;
    107   for (int i = 0; i < Mag->n; i++) {
    108     if (mask->data.U8[i]) continue;
    109     if (isnan(mMinValid) || (Mag->data.F32[i] < mMinValid)) {
    110       mMinValid = Mag->data.F32[i];
     101    if (!psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag)) {
     102        psLogMsg ("psastro", 4, "Failed to fit the luminosity function\n");
     103        return(NULL);
    111104    }
    112     if (isnan(mMaxValid) || (Mag->data.F32[i] > mMaxValid)) {
    113       mMaxValid = Mag->data.F32[i];
     105
     106    // find min and max unmasked magnitudes
     107    double mMinValid = NAN;
     108    double mMaxValid = NAN;
     109    for (int i = 0; i < Mag->n; i++) {
     110        if (mask->data.U8[i]) continue;
     111        if (isnan(mMinValid) || (Mag->data.F32[i] < mMinValid)) {
     112            mMinValid = Mag->data.F32[i];
     113        }
     114        if (isnan(mMaxValid) || (Mag->data.F32[i] > mMaxValid)) {
     115            mMaxValid = Mag->data.F32[i];
     116        }
    114117    }
    115   }
    116118
    117   psLogMsg ("psastro", 4, "logN = %f + %f mag\n", line->coeff[0], line->coeff[1]);
    118   psLogMsg ("psastro", 4, "logN vs mag residuals: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     119    psLogMsg ("psastro", 4, "logN = %f + %f mag\n", line->coeff[0], line->coeff[1]);
     120    psLogMsg ("psastro", 4, "logN vs mag residuals: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    119121
    120   pmLumFunc *lumFunc = pmLumFuncAlloc ();
    121   lumFunc->mMin = mMinValid;
    122   lumFunc->mMax = mMaxValid;
    123   lumFunc->offset = line->coeff[0];
    124   lumFunc->slope = line->coeff[1];
    125   lumFunc->mPeak = mMin + (iPeak + 0.5)*dMag;
    126   lumFunc->nPeak = nPeak;
    127   lumFunc->sPeak = sPeak;
     122    pmLumFunc *lumFunc = pmLumFuncAlloc ();
     123    lumFunc->mMin = mMinValid;
     124    lumFunc->mMax = mMaxValid;
     125    lumFunc->offset = line->coeff[0];
     126    lumFunc->slope = line->coeff[1];
     127    lumFunc->mPeak = mMin + (iPeak + 0.5)*dMag;
     128    lumFunc->nPeak = nPeak;
     129    lumFunc->sPeak = sPeak;
    128130
    129   psFree (lnMag);
    130   psFree (nMags);
    131   psFree (Mag);
    132   psFree (mask);
    133   psFree (line);
    134   psFree (stats);
     131    psFree (lnMag);
     132    psFree (nMags);
     133    psFree (Mag);
     134    psFree (mask);
     135    psFree (line);
     136    psFree (stats);
    135137
    136   return (lumFunc);
     138    return (lumFunc);
    137139}
    138140
Note: See TracChangeset for help on using the changeset viewer.