Index: trunk/psastro/src/psastroLuminosityFunction.c
===================================================================
--- trunk/psastro/src/psastroLuminosityFunction.c	(revision 9587)
+++ trunk/psastro/src/psastroLuminosityFunction.c	(revision 9627)
@@ -15,18 +15,18 @@
   psMemSetDeallocator(func, (psFreeFunc) pmLumFuncFree);
 
-  func->mMin = 0;
-  func->mMax = 0;
-  func->slope = 0;
-  func->offset = 0;
+  func->mMin = mMin;
+  func->mMax = mMax;
+  func->slope = slope;
+  func->offset = offset;
 
   return func;
 }
 
-pmLumFunc *psastroRefstarSubset (psArray *stars) {
+pmLumFunc *psastroLuminosityFunction (psArray *stars) {
 
   // determine the max and min magnitude for the array of stars
   pmAstromObj *star = stars->data[0];
-  mMin = star->Mag;
-  mMax = star->Mag;
+  double mMin = star->Mag;
+  double mMax = star->Mag;
   for (int i = 0; i < stars->n; i++) {
     star = stars->data[i];
@@ -35,5 +35,5 @@
   }
 
-  int nBin = 1 + (nMax - nMin) / dMag;
+  int nBin = 1 + (mMax - mMin) / dMag;
   if (nBin <= 1) {
     psError(PSASTRO_ERR_DATA, true, "magnitude range of 0.0\n");
@@ -45,8 +45,20 @@
   // bin[i] = mMin + i*dMag : mMin + (i+1)*dMag
   psVector *nMags = psVectorAlloc (nBin, PS_TYPE_F32);
+  nMags->n = nBin;
+  psVectorInit (nMags, 0);
+
   for (int i = 0; i < stars->n; i++) {
     star = stars->data[i];
     int bin = (star->Mag - mMin) / dMag;
     nMags->data.F32[bin] += 1.0;
+  }
+
+  // find the peak and position
+  int iPeak = 0;
+  int nPeak = 0;
+  for (int i = 0; i < nMags->n; i++) {
+    if (nMags->data.F32[i] < nPeak) continue;
+    iPeak = i;
+    nPeak = nMags->data.F32[i];
   }
 
@@ -56,13 +68,21 @@
   // create 2 vectors represnting the dlogN/dlogS line
   // exclude bins with no stars
+  // exclude points after the peak with N < 0.8*peak
   int n = 0;
-  for (int i = 0; i < stars->n; i++) {
+  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++;
   }
+  lnMag->n = n;
+  Mag->n = n;
+  psLogMsg ("psastro", 4, "fitting %d points to luminosity function\n", n);
 
   psVector *mask = psVectorAlloc (nBin, PS_TYPE_MASK);
+  mask->n = Mag->n;
+  psVectorInit (mask, 0);
+
   psPolynomial1D *line = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
   psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
@@ -70,9 +90,9 @@
   stats->clipIter = 3;
 
-  line = psVectorClipFitPolynomial1D(line, stats, mask, 1, lnMag, NULL, Mag);
+  line = psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag);
 
   // find min and max unmasked magnitudes
-  mMinValid = NAN;
-  mMaxValid = NAN;
+  double mMinValid = NAN;
+  double mMaxValid = NAN;
   for (int i = 0; i < Mag->n; i++) {
     if (mask->data.U8[i]) continue;
@@ -89,4 +109,12 @@
 
   pmLumFunc *lumFunc = pmLumFuncAlloc (mMinValid, mMaxValid, line->coeff[0], line->coeff[1]);
+
+  psFree (lnMag);
+  psFree (nMags);
+  psFree (Mag);
+  psFree (mask);
+  psFree (line);
+  psFree (stats);
+
   return (lumFunc);
 }