Changeset 20036
- Timestamp:
- Oct 9, 2008, 3:07:38 PM (18 years ago)
- Location:
- trunk
- Files:
-
- 1 added
- 3 edited
-
psModules/src/astrom/pmAstrometryObjects.c (modified) (9 diffs)
-
psastro/src/psastro.h (modified) (3 diffs)
-
psastro/src/psastroLuminosityFunction.c (modified) (8 diffs)
-
psastro/src/psastroLuminosityFunctionPlot.c (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/astrom/pmAstrometryObjects.c
r17437 r20036 8 8 * @author EAM, IfA 9 9 * 10 * @version $Revision: 1. 39$ $Name: not supported by cvs2svn $11 * @date $Date: 2008- 04-11 07:42:05$10 * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2008-10-10 01:07:38 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 17 17 #include <config.h> 18 18 #endif 19 19 20 20 21 /******************************************************************************/ … … 28 29 #include <unistd.h> // for unlink 29 30 #include <pslib.h> 31 #include <kapa.h> // cnb: do I need an IFDEF? 30 32 31 33 #include "pmHDU.h" 32 34 #include "pmFPA.h" 33 35 #include "pmAstrometryObjects.h" 36 #include "pmKapaPlots.h" 34 37 35 38 #define PM_ASTROMETRYOBJECTS_DEBUG 1 39 36 40 /****************************************************************************** 37 41 pmAstromObjSortByMag(**a, **b): sort by mag (descending) … … 243 247 // need to use the stats lookups functions to get the width and center 244 248 for (int i = 0; i < nIter; i++) { 245 249 if (!psVectorClipFitPolynomial2D (map->x, results->xStats, mask, 0xff, x, wt, X, Y)) { 246 250 psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for x\n"); 247 248 249 250 251 252 253 254 255 251 psFree (x); 252 psFree (y); 253 psFree (X); 254 psFree (Y); 255 psFree (wt); 256 psFree (mask); 257 258 return results; 259 } 256 260 psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n); 257 261 258 262 if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) { 259 263 psError(PS_ERR_UNKNOWN, false, "failure in clip-fitting for y\n"); 260 261 262 263 264 265 266 267 268 264 psFree (x); 265 psFree (y); 266 psFree (X); 267 psFree (Y); 268 psFree (wt); 269 psFree (mask); 270 271 return results; 272 } 269 273 psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n); 270 274 } … … 514 518 const psMetadata *config) 515 519 { 520 516 521 PS_ASSERT_PTR_NON_NULL(raw, NULL); 517 522 PS_ASSERT_PTR_NON_NULL(ref, NULL); … … 559 564 psF32 **D2 = gridD2->data.F32; 560 565 566 // vectors to hold dX and dY 567 int nplot = raw->n * ref->n; 568 float dXplot[nplot]; 569 float dYplot[nplot]; 570 561 571 // accumulate grids for focal plane (L,M) matches 562 572 for (int i = 0; i < raw->n; i++) { … … 567 577 dY = ob1->FP->y - ob2->FP->y; 568 578 579 dXplot[(i * ref->n) + j] = dX; 580 dYplot[(i * ref->n) + j] = dY; 569 581 // fprintf (f, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP->x, ob2->FP->x, ob1->FP->y, ob2->FP->y); 570 582 // find bin coordinates for this delta-delta … … 673 685 // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]); 674 686 675 psFree (sort); 687 /*** plot DXs and DYs ***/ 688 KapaSection section = {"s1", 0.00, 0.00, 1.0, 1.0}; 689 //KapaSection left = {"s2", 0.05, 0.05, 0.35, 0.8}; 690 //KapaSection right = {"s3", 0.55, 0.05, 0.35, 0.8}; 691 //KapaSection yprof = {"s4", 0.45, 0.05, 0.10, 0.9}; 692 //KapaSection xprof = {"s5", .05, .90, .35, .1}; 693 694 Graphdata graphdata; 695 int kapa = pmKapaOpen(true); 696 KapaClearPlots(kapa); 697 KapaInitGraph(&graphdata); 698 KapaSetSection(kapa, §ion); 699 graphdata.xmin = stats->offset.x - 1.5 * maxOffpix; 700 graphdata.xmax = stats->offset.x + 1.5 * maxOffpix; 701 graphdata.ymin = stats->offset.y - 1.5 * maxOffpix; 702 graphdata.ymax = stats->offset.y + 1.5 * maxOffpix; 703 704 KapaSetLimits(kapa, &graphdata); 705 KapaSetFont(kapa, "helvetica", 14); 706 KapaBox(kapa, &graphdata); 707 KapaSendLabel (kapa, "X offset", KAPA_LABEL_XM); 708 KapaSendLabel (kapa, "Y offset", KAPA_LABEL_YM); 709 KapaSendLabel (kapa, "pmAstromGridAngle residuals. Big Box: Serach Region. Small Box: Correlation Peak.", 710 KAPA_LABEL_XP); 711 graphdata.style = 2; 712 graphdata.ptype = 0; 713 graphdata.size = 0.4; 714 graphdata.color = KapaColorByName ("black"); 715 KapaPrepPlot(kapa, nplot, &graphdata); 716 KapaPlotVector (kapa, nplot, dXplot, "x"); 717 KapaPlotVector (kapa, nplot, dYplot, "y"); 718 719 //Overplot bounding box, peak of distribution 720 float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix}; 721 float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix}; 722 float xbin[5] = {stats->offset.x - 0.5 * Scale, stats->offset.x + 0.5 * Scale, stats->offset.x + 0.5 * Scale, 723 stats->offset.x - 0.5 * Scale, stats->offset.x - 0.5 * Scale}; 724 float ybin[5] = {stats->offset.y - 0.5 * Scale, stats->offset.y - 0.5 * Scale, stats->offset.y + 0.5 * Scale, 725 stats->offset.y + 0.5 * Scale, stats->offset.y - 0.5 * Scale}; 726 graphdata.color = KapaColorByName("red"); 727 graphdata.style = 0; 728 graphdata.size = 1.0; 729 KapaPrepPlot(kapa, 5, &graphdata); 730 KapaPlotVector (kapa, 5, xbound, "x"); 731 KapaPlotVector (kapa, 5, ybound, "y"); 732 KapaPrepPlot(kapa, 5, &graphdata); 733 KapaPlotVector (kapa, 5, xbin, "x"); 734 KapaPlotVector (kapa, 5, ybin, "y"); 735 736 KapaPNG(kapa, "dXdY.png"); 737 char c; 738 fscanf(stdin, "%c", &c); 739 /*** end ***/ 740 741 psFree (sort); 676 742 psFree (listNP); 677 743 psFree (gridNP); … … 730 796 minStat->minMetric = 1e10; 731 797 for (double scale = minScale; scale <= maxScale; scale += delScale) { 732 733 798 for (double angle = minAngle; angle <= maxAngle; angle += delAngle) { 799 rot = pmAstromRotateObj (raw, center, angle, scale); 734 800 735 801 # if 0 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 802 FILE *f1 = fopen ("raw.dat", "w"); 803 for (int i = 0; i < rot->n; i++) { 804 pmAstromObj *obj = rot->data[i]; 805 fprintf (f1, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); 806 } 807 fclose (f1); 808 FILE *f2 = fopen ("ref.dat", "w"); 809 for (int i = 0; i < ref->n; i++) { 810 pmAstromObj *obj = ref->data[i]; 811 fprintf (f2, "%8.2f %8.2f %6.2f\n", obj->FP->x, obj->FP->y, obj->Mag); 812 } 813 fclose (f2); 814 fprintf (stderr, "type return"); 815 char c; 816 fscanf (stdin, "%c", &c); 751 817 # endif 752 753 newStat = pmAstromGridAngle (rot, ref, config); 754 newStat->angle = angle; 755 newStat->scale = scale; 756 newStat->center = center; 757 758 if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) { 759 *minStat = *newStat; 760 psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", 761 minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); 762 } else { 763 psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", 764 newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); 765 766 } 767 psFree (newStat); 768 psFree (rot); 769 } 818 newStat = pmAstromGridAngle (rot, ref, config); 819 newStat->angle = angle; 820 newStat->scale = scale; 821 newStat->center = center; 822 823 if (isfinite(newStat->minMetric) && (newStat->minMetric > 0.0) && (newStat->minMetric < minStat->minMetric)) { 824 *minStat = *newStat; 825 psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric) *", 826 minStat->offset.x, minStat->offset.y, PS_DEG_RAD*minStat->angle, minStat->scale, minStat->nMatch, minStat->nSigma, minStat->minVar, log10(minStat->minMetric)); 827 } else { 828 psLogMsg ("psModule.astrom", 4, "grid test - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", 829 newStat->offset.x, newStat->offset.y, PS_DEG_RAD*newStat->angle, newStat->scale, newStat->nMatch, newStat->nSigma, newStat->minVar, log10(newStat->minMetric)); 830 831 } 832 psFree (newStat); 833 psFree (rot); 834 } 770 835 } 771 836 psLogMsg ("psModule.astrom.grid.match", 4, "grid best - offset: %7.2f,%7.2f @ %6.1f deg x %7.3f (%4d pts, %5.1f sig, %5.1f var, %6.3f log metric)", -
trunk/psastro/src/psastro.h
r19512 r20036 17 17 // logN = offset + slope * logS 18 18 typedef struct { 19 double mMin; 20 double mMax; 21 double offset; 22 double slope; 23 double mPeak; // mag of peak bin24 int nPeak; 25 int sPeak; 19 double mMin; // minimum magnitude bin with data 20 double mMax; // maximum magnitude bin with data 21 double offset; // fitted line offset 22 double slope; // fitted line slope 23 double mPeak; // mag of peak bin 24 int nPeak; // # of stars in peak bin 25 int sPeak; // sum of stars to peak bin 26 26 } pmLumFunc; 27 27 … … 47 47 bool psastroChooseRefstars (pmConfig *config, psArray *refs); 48 48 bool psastroRefstarSubset (pmReadout *readout); 49 pmLumFunc *psastroLuminosityFunction (psArray *stars); 49 pmLumFunc *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc); 50 bool psastroLuminosityFunctionPlot(psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc); 50 51 psArray *psastroRemoveClumps (psArray *input, int scale); 51 52 … … 56 57 57 58 // mosaic fitting functions 58 bool 59 bool psastroMosaicDistortion (pmFPA *fpa, psMetadata *recipe, int pass); 59 60 psPlaneTransform *psastroMosaicFitRotAndScale (pmFPA *fpa); 60 61 bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP); 61 bool 62 bool psastroMosaicDistortionFromGradients (pmFPA *fpa, psMetadata *recipe); 62 63 bool psastroMosaicCorrectDistortion (pmFPA *fpa, psPlaneTransform *TPtoFP); 63 bool 64 bool 65 bool 66 bool 67 bool 68 bool 69 bool 70 bool 64 bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe); 65 bool psastroMosaicAstrom (pmConfig *config); 66 bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration); 67 bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration); 68 bool psastroMosaicSetAstrom (pmFPA *fpa); 69 bool psastroMosaicHeaders (pmConfig *config); 70 bool psastroMosaicRescaleChips (pmFPA *fpa); 71 bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration); 71 72 72 73 // Return version strings. 73 psString 74 psString 74 psString psastroVersion(void); 75 psString psastroVersionLong(void); 75 76 76 77 // demo plots 77 bool 78 bool 79 bool 78 bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe); 79 bool psastroPlotRefstars (psArray *refstars, psMetadata *recipe); 80 bool psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe); 80 81 81 bool 82 bool 83 bool 84 bool 82 bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip); 83 bool psastroDumpRefstars (psArray *refstars, char *filename); 84 bool psastroDumpMatches (pmFPA *fpa, char *filename); 85 bool psastroDumpStars (psArray *stars, char *filename); 85 86 bool psastroDumpMatchedStars (char *filename, psArray *rawstars, psArray *refstars, psArray *match); 86 87 bool psastroDumpGradients (psArray *gradients, char *filename); 87 88 88 bool 89 int 89 bool psastroMosaicSetAstrom_tmp (pmFPA *fpa); 90 int psastroSortByMag (const void *a, const void *b); 90 91 91 92 bool psastroFixChips (pmConfig *config, psMetadata *recipe); -
trunk/psastro/src/psastroLuminosityFunction.c
r17677 r20036 27 27 } 28 28 29 pmLumFunc *psastroLuminosityFunction (psArray *stars ) {29 pmLumFunc *psastroLuminosityFunction (psArray *stars, pmLumFunc *rawFunc) { 30 30 31 31 if (stars->n == 0) return NULL; … … 36 36 double mMax = star->Mag; 37 37 for (int i = 0; i < stars->n; i++) { 38 39 40 38 star = stars->data[i]; 39 mMin = PS_MIN (star->Mag, mMin); 40 mMax = PS_MAX (star->Mag, mMax); 41 41 } 42 42 43 43 int nBin = 1 + (mMax - mMin) / dMag; 44 44 if (nBin <= 1) { 45 46 45 psLogMsg ("psastro", 4, "insufficient valid stars for this readout\n"); 46 return NULL; 47 47 } 48 48 … … 54 54 55 55 for (int i = 0; i < stars->n; i++) { 56 57 58 59 60 61 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 62 } 63 63 … … 67 67 int sPeak = 0; 68 68 for (int i = 0; i < nMags->n; i++) { 69 70 71 72 69 if (nMags->data.F32[i] < nPeak) continue; 70 iPeak = i; 71 nPeak = nMags->data.F32[i]; 72 sPeak += nPeak; 73 73 } 74 74 … … 81 81 int n = 0; 82 82 for (int i = 0; i < nMags->n; i++) { 83 84 85 86 87 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 88 } 89 89 lnMag->n = n; … … 100 100 101 101 if (!psVectorClipFitPolynomial1D(line, stats, mask, 0xff, lnMag, NULL, Mag)) { 102 103 102 psLogMsg ("psastro", 4, "Failed to fit the luminosity function\n"); 103 return(NULL); 104 104 } 105 105 … … 108 108 double mMaxValid = NAN; 109 109 for (int i = 0; i < Mag->n; i++) { 110 111 112 113 114 115 116 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 } 117 117 } 118 118 … … 129 129 lumFunc->sPeak = sPeak; 130 130 131 psastroLuminosityFunctionPlot(lnMag, Mag, lumFunc, rawFunc); 132 131 133 psFree (lnMag); 132 134 psFree (nMags);
Note:
See TracChangeset
for help on using the changeset viewer.
