Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 19977) +++ trunk/psastro/src/psastro.h (revision 20036) @@ -17,11 +17,11 @@ // logN = offset + slope * logS typedef struct { - double mMin; // minimum magnitude bin with data - double mMax; // maximum magnitude bin with data - double offset; // fitted line offset - double slope; // fitted line slope - double mPeak; // mag of peak bin - int nPeak; // # of stars in peak bin - int sPeak; // sum of stars to peak bin + double mMin; // minimum magnitude bin with data + double mMax; // maximum magnitude bin with data + double offset; // fitted line offset + double slope; // fitted line slope + double mPeak; // mag of peak bin + int nPeak; // # of stars in peak bin + int sPeak; // sum of stars to peak bin } pmLumFunc; @@ -47,5 +47,6 @@ bool psastroChooseRefstars (pmConfig *config, psArray *refs); bool psastroRefstarSubset (pmReadout *readout); -pmLumFunc *psastroLuminosityFunction (psArray *stars); +pmLumFunc *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc); +bool psastroLuminosityFunctionPlot(psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc); psArray *psastroRemoveClumps (psArray *input, int scale); @@ -56,36 +57,36 @@ // mosaic fitting functions -bool psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe, int pass); +bool psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe, int pass); psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa); bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP); -bool psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe); +bool psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP); -bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe); -bool psastroMosaicAstrom (pmConfig *config); -bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration); -bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration); -bool psastroMosaicSetAstrom (pmFPA *fpa); -bool psastroMosaicHeaders (pmConfig *config); -bool psastroMosaicRescaleChips (pmFPA *fpa); -bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration); +bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe); +bool psastroMosaicAstrom (pmConfig *config); +bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration); +bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration); +bool psastroMosaicSetAstrom (pmFPA *fpa); +bool psastroMosaicHeaders (pmConfig *config); +bool psastroMosaicRescaleChips (pmFPA *fpa); +bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration); // Return version strings. -psString psastroVersion(void); -psString psastroVersionLong(void); +psString psastroVersion(void); +psString psastroVersionLong(void); // demo plots -bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); -bool psastroPlotRefstars (psArray *refstars, psMetadata *recipe); -bool psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); +bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); +bool psastroPlotRefstars (psArray *refstars, psMetadata *recipe); +bool psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); -bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip); -bool psastroDumpRefstars (psArray *refstars, char *filename); -bool psastroDumpMatches (pmFPA *fpa, char *filename); -bool psastroDumpStars (psArray *stars, char *filename); +bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip); +bool psastroDumpRefstars (psArray *refstars, char *filename); +bool psastroDumpMatches (pmFPA *fpa, char *filename); +bool psastroDumpStars (psArray *stars, char *filename); bool psastroDumpMatchedStars (char *filename, psArray *rawstars, psArray *refstars, psArray *match); bool psastroDumpGradients (psArray *gradients, char *filename); -bool psastroMosaicSetAstrom_tmp (pmFPA *fpa); -int psastroSortByMag (const void *a, const void *b); +bool psastroMosaicSetAstrom_tmp (pmFPA *fpa); +int psastroSortByMag (const void *a, const void *b); bool psastroFixChips (pmConfig *config, psMetadata *recipe); Index: trunk/psastro/src/psastroLuminosityFunction.c =================================================================== --- trunk/psastro/src/psastroLuminosityFunction.c (revision 19977) +++ trunk/psastro/src/psastroLuminosityFunction.c (revision 20036) @@ -27,5 +27,5 @@ } -pmLumFunc *psastroLuminosityFunction (psArray *stars) { +pmLumFunc *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc) { if (stars->n == 0) return NULL; @@ -36,13 +36,13 @@ 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,10 +54,10 @@ 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,8 +67,8 @@ 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,9 +81,9 @@ 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,6 +100,6 @@ 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,11 +108,11 @@ 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,4 +129,6 @@ lumFunc->sPeak = sPeak; + psastroLuminosityFunctionPlot(lnMag, Mag, lumFunc, rawFunc); + psFree (lnMag); psFree (nMags); Index: trunk/psastro/src/psastroLuminosityFunctionPlot.c =================================================================== --- trunk/psastro/src/psastroLuminosityFunctionPlot.c (revision 20036) +++ trunk/psastro/src/psastroLuminosityFunctionPlot.c (revision 20036) @@ -0,0 +1,100 @@ +#include "psastroInternal.h" + +bool psastroLuminosityFunctionPlot (psVector *lnMag, psVector *Mag, + pmLumFunc *lumFunc, pmLumFunc *rawFunc) { + //Set up graph window + Graphdata graphdata; + int kapa = pmKapaOpen(true); + KapaSection section1 = {"s1", .05, .05, .9, .4}; + KapaSection section2 = {"s2", .05, .55, .9, .4}; + KapaSection section; + if(rawFunc == NULL) + section = section1; + else + section = section2; + int DX = 1000, DY = 800; + KapaResize(kapa, DX, DY); + if (rawFunc == NULL) + KapaClearPlots(kapa); + KapaInitGraph(&graphdata); + + //Determine Plot Limits + graphdata.xmin = 50; + graphdata.ymin = 50; + graphdata.xmax = -50; + graphdata.ymax = -50; + + for(int i = 0; i < lnMag->n; i++) { + graphdata.xmin = PS_MIN(graphdata.xmin, Mag->data.F32[i]); + graphdata.ymin = PS_MIN(graphdata.ymin, lnMag->data.F32[i]); + graphdata.xmax = PS_MAX(graphdata.xmax, Mag->data.F32[i]); + graphdata.ymax = PS_MAX(graphdata.ymax, lnMag->data.F32[i]); + } + + float range = graphdata.xmax - graphdata.xmin; + graphdata.xmin -= .05 * range; + graphdata.xmax += .05 * range; + + range = graphdata.ymax - graphdata.ymin; + graphdata.ymin -= .05 * range; + graphdata.ymax += .05 * range; + + //Make a line for the fit + float x[2] = {graphdata.xmin, graphdata.xmax}; + float y[2] = {lumFunc->offset + x[0] * lumFunc->slope, + lumFunc->offset + x[1] * lumFunc->slope}; + + //Plot Data + KapaSetSection(kapa, §ion); + KapaSetLimits(kapa, &graphdata); + + KapaSetFont (kapa, "helvetica", 14); + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM); + KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM); + if (rawFunc == NULL) + KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP); + else + KapaSendLabel (kapa, + "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff", + KAPA_LABEL_XP); + graphdata.color=KapaColorByName("black"); + graphdata.style = 1; + KapaPrepPlot (kapa, lnMag->n, &graphdata); + KapaPlotVector(kapa, lnMag->n, Mag->data.F32, "x"); + KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y"); + + //Overplot fit + graphdata.style=0; + KapaPrepPlot(kapa,2,&graphdata); + KapaPlotVector(kapa, 2, x, "x"); + KapaPlotVector(kapa, 2, y, "y"); + + //If rawFunc was supplied, calculate the cutoff magnitude + if( rawFunc != NULL) { + double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax); + double logRho = mRef * lumFunc->slope + lumFunc->offset; + double mRaw = (logRho - rawFunc->offset) / rawFunc->slope; + double deltaM = mRef - mRaw; + double mRefMax = rawFunc->mMax + deltaM; + + float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM}; + float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin, + rawFunc->offset + (rawFunc->slope) * rawFunc->mMax}; + float x[2] = {mRefMax, mRefMax}; + float y[2] = {graphdata.ymin, graphdata.ymax}; + graphdata.color= KapaColorByName("red"); + KapaPrepPlot(kapa, 2, &graphdata); + KapaPlotVector(kapa, 2, x, "x"); + KapaPlotVector(kapa, 2, y, "y"); + KapaPrepPlot (kapa, 2, &graphdata); + KapaPlotVector (kapa, 2, xraw, "x"); + KapaPlotVector (kapa, 2, yraw, "y"); + + KapaPNG(kapa,"LumPlot.png"); + //char c; + //fscanf(stdin, "%c", &c); + // pmKapaClose(kapa); + } + return true; +}