Changeset 20036 for trunk/psastro/src/psastroLuminosityFunction.c

Timestamp:

Oct 9, 2008, 3:07:38 PM (18 years ago)

Author:

Paul Price

Message:

Merging from Chris Beaumont branch (cns_branch_20081009) to get luminosity function plots.

File:

• trunk/psastro/src/psastroLuminosityFunction.c (modified) (8 diffs)

trunk/psastro/src/psastroLuminosityFunction.c

@@ -27 +27 @@
 }
 
-pmLumFunc *psastroLuminosityFunction (psArray *stars) {
+pmLumFunc *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc) {
 
     if (stars->n == 0) return NULL;
@@ -36 +36 @@
     double mMax = star->Mag;
     for (int i = 0; i < stars->n; i++) {
-        star = stars->data[i];
-        mMin = PS_MIN (star->Mag, mMin);
-        mMax = PS_MAX (star->Mag, mMax);
+        star = stars->data[i];
+        mMin = PS_MIN (star->Mag, mMin);
+        mMax = PS_MAX (star->Mag, mMax);
     }
 
     int nBin = 1 + (mMax - mMin) / dMag;
     if (nBin <= 1) {
-        psLogMsg ("psastro", 4, "insufficient valid stars for this readout\n");
-        return NULL;
+        psLogMsg ("psastro", 4, "insufficient valid stars for this readout\n");
+        return NULL;
     }
 
@@ -54 +54 @@
 
     for (int i = 0; i < stars->n; i++) {
-        star = stars->data[i];
-        if (!isfinite(star->Mag)) continue;
-        int bin = (star->Mag - mMin) / dMag;
-        if (bin < 0) psAbort("bin cannot be negative!");
-        if (bin >= nBin) psAbort("bin cannot be > %d!", nBin);
-        nMags->data.F32[bin] += 1.0;
+        star = stars->data[i];
+        if (!isfinite(star->Mag)) continue;
+        int bin = (star->Mag - mMin) / dMag;
+        if (bin < 0) psAbort("bin cannot be negative!");
+        if (bin >= nBin) psAbort("bin cannot be > %d!", nBin);
+        nMags->data.F32[bin] += 1.0;
     }
 
@@ -67 +67 @@
     int sPeak = 0;
     for (int i = 0; i < nMags->n; i++) {
-        if (nMags->data.F32[i] < nPeak) continue;
-        iPeak = i;
-        nPeak = nMags->data.F32[i];
-        sPeak += nPeak;
+        if (nMags->data.F32[i] < nPeak) continue;
+        iPeak = i;
+        nPeak = nMags->data.F32[i];
+        sPeak += nPeak;
     }
 
@@ -81 +81 @@
     int n = 0;
     for (int i = 0; i < nMags->n; i++) {
-        if (nMags->data.F32[i] < 1) continue;
-        if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue;
-        lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
-        Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
-        n++;
+        if (nMags->data.F32[i] < 1) continue;
+        if ((i > iPeak) && (nMags->data.F32[i] < 0.8*nPeak)) continue;
+        lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
+        Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
+        n++;
     }
     lnMag->n = n;
@@ -100 +100 @@
 
     if (!psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag)) {
-        psLogMsg ("psastro", 4, "Failed to fit the luminosity function\n");
-        return(NULL);
+        psLogMsg ("psastro", 4, "Failed to fit the luminosity function\n");
+        return(NULL);
     }
 
@@ -108 +108 @@
     double mMaxValid = NAN;
     for (int i = 0; i < Mag->n; i++) {
-        if (mask->data.U8[i]) continue;
-        if (isnan(mMinValid) || (Mag->data.F32[i] < mMinValid)) {
-            mMinValid = Mag->data.F32[i];
-        }
-        if (isnan(mMaxValid) || (Mag->data.F32[i] > mMaxValid)) {
-            mMaxValid = Mag->data.F32[i];
-        }
+        if (mask->data.U8[i]) continue;
+        if (isnan(mMinValid) || (Mag->data.F32[i] < mMinValid)) {
+            mMinValid = Mag->data.F32[i];
+        }
+        if (isnan(mMaxValid) || (Mag->data.F32[i] > mMaxValid)) {
+            mMaxValid = Mag->data.F32[i];
+        }
     }
 
@@ -129 +129 @@
     lumFunc->sPeak = sPeak;
 
+    psastroLuminosityFunctionPlot(lnMag, Mag, lumFunc, rawFunc);
+
     psFree (lnMag);
     psFree (nMags);