Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5836) +++ /trunk/psphot/Makefile (revision 5837) @@ -97,5 +97,5 @@ $(MODELTEST): $(SRC)/psphot.h -test.install: psphotTest.install +test-install: psphotTest.install test: $(BIN)/psphotTest.$(ARCH) $(BIN)/psphotTest.$(ARCH) : $(TEST) Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 5836) +++ /trunk/psphot/src/psphotApResid.c (revision 5837) @@ -5,6 +5,6 @@ int Npsf; bool status; - float x, y; - float fitMag, obsMag; + pmModel *model; + pmSource *source; // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS @@ -26,8 +26,6 @@ // select the NNN brightest, non-saturated sources, or just select PSFSTARs? for (int i = 0; (i < sources->n) && (Npsf < 300); i++) { - pmSource *source = sources->data[i]; - pmModel *model = source->modelPSF; + source = sources->data[i]; - if (model == NULL) continue; if (source->type != PM_SOURCE_STAR) continue; if (source->mode & PM_SOURCE_SATSTAR) continue; @@ -36,22 +34,13 @@ if (source->mode & PM_SOURCE_POOR) continue; - // set shape for this model based on PSF - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - - pmSourceAddModel (source->pixels, source->mask, model, false, false); - - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - pmSourceSubModel (source->pixels, source->mask, model, false, false); - if (!status) continue; + // get magnitudes, uncorrected for (x, y, rflux) + model = pmSourceMagnitudes (source, NULL, RADIUS); + if (model == NULL) continue; mask->data.U8[Npsf] = 0; - xPos->data.F64[Npsf] = x; - yPos->data.F64[Npsf] = y; - rflux->data.F64[Npsf] = pow(10.0, 0.4*fitMag); - apResid->data.F64[Npsf] = obsMag - fitMag; + xPos->data.F64[Npsf] = model->params->data.F32[2]; + yPos->data.F64[Npsf] = model->params->data.F32[3]; + rflux->data.F64[Npsf] = pow(10.0, 0.4*source->fitMag); + apResid->data.F64[Npsf] = source->apMag - source->fitMag; psVectorExtend (mask, 100, 1); Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 5836) +++ /trunk/psphot/src/psphotArguments.c (revision 5837) @@ -61,5 +61,6 @@ bool ModelTest = false; - float ModelTest_X, ModelTest_Y; + float ModelTest_X = 0; + float ModelTest_Y = 0; char *model = NULL; char *fitset = NULL; Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 5836) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5837) @@ -74,5 +74,5 @@ bool UseAnalysisRegion = false; - psRegion AnalysisRegion; + psRegion AnalysisRegion = {0, 0, 0, 0}; char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION"); if (status) { @@ -236,5 +236,5 @@ // XXX EAM : need to free up many things here - psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f\n", psTimerMark ("psphot")); + psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot")); return true; } Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 5836) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 5837) @@ -1,5 +1,5 @@ # include "psphot.h" -// XXX EAM : the apMag should only be calculated for the brighter sources +// XXX EAM : the apMag should only be calculated for the brighter sources? // XXX EAM : SN limit set by user? // XXX EAM : masked region should be (optionally) elliptical @@ -7,26 +7,37 @@ int status; + bool isPSF; float x, y; float rflux; + float radius; pmModel *model; - // use the correct model (PSF vs FLT) - if (psf != NULL) { - model = source->modelPSF; - } else { - model = source->modelFLT; + switch (source->type) { + case PM_SOURCE_STAR: + model = source->modelPSF; + if (model == NULL) return NULL; + radius = apRadius; + isPSF = true; + break; + + case PM_SOURCE_GALAXY: + model = source->modelFLT; + if (model == NULL) return NULL; + radius = model->radius; + isPSF = false; + break; + + default: + return NULL; } - if (model == NULL) return NULL; - x = source->peak->x; - y = source->peak->y; - - // we have already (psphotApplyPSF) defined pixels at least OUTER_RADIUS from source - // we need to mask pixels to measure the aperture magnitude - // set aperture mask circle of PSF_FIT_RADIUS - psImageKeepCircle (source->mask, x, y, apRadius, "OR", PSPHOT_MASK_MARKED); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; // replace source flux pmSourceAddModel (source->pixels, source->mask, model, false, false); + + // set aperture mask circle of PSF_FIT_RADIUS + psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED); // measure object photometry @@ -34,5 +45,5 @@ // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS - if (psf != NULL) { + if (isPSF && (psf != NULL)) { rflux = pow (10.0, 0.4*source->fitMag); source->apMag -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius)); @@ -44,29 +55,8 @@ // unmask aperture - psImageKeepCircle (source->mask, x, y, apRadius, "AND", ~PSPHOT_MASK_MARKED); + psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); + + if (!status) return NULL; return model; -} - -pmModel *pmSourceSelectModel (pmSource *source) { - - pmModel *model; - - switch (source->type) { - - // use PSF model with stars - case PM_SOURCE_STAR: - model = source->modelPSF; - break; - - // use FLT model with galaxies - case PM_SOURCE_GALAXY: - model = source->modelFLT; - break; - - // skip defects and poorly fitted stars or galaxies - default: - model = NULL; - } - return (model); } Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 5836) +++ /trunk/psphot/src/psphotModelTest.c (revision 5837) @@ -8,5 +8,5 @@ float obsMag, fitMag, value; char name[64]; - pmPSF *psf; + pmPSF *psf = NULL; psMetadataItem *item = NULL; Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5836) +++ /trunk/psphot/src/psphotOutput.c (revision 5837) @@ -24,4 +24,20 @@ return; } + if (!strcasecmp (outputMode, "OBJ")) { + pmSourcesWriteOBJ (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "SX")) { + pmSourcesWriteSX (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "CMP")) { + pmSourcesWriteCMP (imdata, config, outputFile, sources, psf, sky); + return; + } + if (!strcasecmp (outputMode, "CMF")) { + pmSourcesWriteCMF (imdata, config, outputFile, sources, psf, sky); + return; + } if (!strcasecmp (outputMode, "TEXT")) { pmSourcesWriteText (imdata, config, outputFile, sources, psf, sky); @@ -29,45 +45,5 @@ } - if (!strcasecmp (outputMode, "OBJ")) { - pmSourcesWriteOBJ (imdata, config, outputFile, sources, psf, sky); - return; - } - - if (!strcasecmp (outputMode, "SX")) { - pmSourcesWriteSX (imdata, config, outputFile, sources, psf, sky); - return; - } - - if (!strcasecmp (outputMode, "CMP")) { - pmSourcesWriteCMP (imdata, config, outputFile, sources, psf, sky); - return; - } - - if (!strcasecmp (outputMode, "CMF")) { - pmSourcesWriteCMF (imdata, config, outputFile, sources, psf, sky); - return; - } - psAbort ("psphot", "unknown output mode %s", outputMode); -} - -bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky) { - - char *name = (char *) psAlloc (strlen(filename) + 10); - - sprintf (name, "%s.psf.dat", filename); - pmModelWritePSFs (sources, config, name, psf); - - sprintf (name, "%s.flt.dat", filename); - pmModelWriteFLTs (sources, name); - - sprintf (name, "%s.nul.dat", filename); - pmModelWriteNULLs (sources, name); - - sprintf (name, "%s.mnt.dat", filename); - pmMomentsWriteText (sources, name); - - psFree (name); - return true; } @@ -115,5 +91,5 @@ psLineAdd (line, "%9.3f", PAR[5]); psLineAdd (line, "%7.2f", PAR[6]); - psLineAdd (line, "%8.3f", 32.0); + psLineAdd (line, "%8.3f", 99.999); psLineAdd (line, "%8.3f", source->apMag); psLineAdd (line, "%8.2f\n", apResid); @@ -166,5 +142,5 @@ psLineAdd (line, "%9.2f", 0.0); // should be FWHMy psLineAdd (line, "%6.1f", 0.0); // should be Theta - psLineAdd (line, "%9.4f", 32.0); // should be MAG_ISO + psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO psLineAdd (line, "%9.4f", source->apMag); psLineAdd (line, "%4d\n", 0); // should be flags @@ -182,5 +158,5 @@ psMetadataItem *mdi; psF32 *PAR, *dPAR; - float dmag, apResid, lsky; + float dmag, lsky; bool status; @@ -188,4 +164,5 @@ float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); + if (!status) ZERO_POINT = 25.0; char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); @@ -237,5 +214,4 @@ dmag = dPAR[1] / PAR[1]; type = pmSourceDophotType (source); - apResid = source->apMag - source->fitMag; lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); @@ -243,10 +219,10 @@ psLineAdd (line, "%6.1f ", PAR[2]); psLineAdd (line, "%6.1f ", PAR[3]); - psLineAdd (line, "%6.3f ", PS_MIN (32.0, source->fitMag + ZERO_POINT)); - psLineAdd (line, "%03d ", (int)(1000*dmag)); + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); psLineAdd (line, "%2d ", type); psLineAdd (line, "%3.1f ", lsky); - psLineAdd (line, "%6.3f ", 32.0); // should be 'Mgal - psLineAdd (line, "%6.3f ", PS_MIN (32.0, source->apMag + ZERO_POINT)); + psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal + psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' @@ -297,5 +273,5 @@ table = psArrayAlloc (sources->n); - sources->n = 0; + table->n = 0; for (i = 0; i < sources->n; i++) { @@ -313,15 +289,15 @@ row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_HEAD, "X_PIX", PS_DATA_F32, "", PAR[2]); - psMetadataAdd (row, PS_LIST_HEAD, "Y_PIX", PS_DATA_F32, "", PAR[3]); - psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_HEAD, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); - psMetadataAdd (row, PS_LIST_HEAD, "MAG_GAL", PS_DATA_F32, "", type); - psMetadataAdd (row, PS_LIST_HEAD, "MAG_AP", PS_DATA_F32, "", lsky); - psMetadataAdd (row, PS_LIST_HEAD, "LOG_SKY", PS_DATA_F32, "", 32.0); // should be 'Mgal' - psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X", PS_DATA_F32, "", source->apMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_HEAD, "FWHM_Y", PS_DATA_F32, "", PAR[4]); // should be 'FHWM x' - psMetadataAdd (row, PS_LIST_HEAD, "THETA", PS_DATA_F32, "", PAR[5]); // should be 'FHWM y' - psMetadataAdd (row, PS_LIST_HEAD, "DOPHOT", PS_TYPE_S8, "", 0); // should be theta + psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", type); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0); // should be 'Mgal' + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", source->apMag + ZERO_POINT); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); // should be 'FHWM x' + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); // should be 'FHWM y' + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_TYPE_S32, "", 0); // should be theta // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY", PS_DATA_STRING, "", NULL); @@ -339,7 +315,19 @@ psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); + psImage *image; + + image = psImageAlloc (512, 512, PS_DATA_F32); + for (int j = 0; j < 512; j++) { + for (int i = 0; i < 512; i++) { + image->data.F32[j][i] = i; + } + } + psFits *fits = psFitsOpen (filename, "w"); - psFitsWriteHeader (imdata->header, fits); - psFitsWriteTable (fits, theader, table); + // psFitsWriteHeader (imdata->header, fits); + // psFitsWriteTable (fits, theader, table); + + psFitsWriteImage (fits, NULL, image, 0); + psFitsWriteTable (fits, NULL, table); psFitsClose (fits); @@ -349,23 +337,21 @@ /***** Text Output Methods *****/ -// write the peaks to an output file -bool pmPeaksWriteText (psArray *peaks, char *filename) { - - int i; - FILE *f; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); - return false; - } - - for (i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - if (peak == NULL) continue; - fprintf (f, "%5d %5d %7.1f\n", - peak->x, peak->y, peak->counts); - } - fclose (f); +bool pmSourcesWriteText (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky) { + + char *name = (char *) psAlloc (strlen(filename) + 10); + + sprintf (name, "%s.psf.dat", filename); + pmModelWritePSFs (sources, config, name, psf); + + sprintf (name, "%s.flt.dat", filename); + pmModelWriteFLTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + pmModelWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + pmMomentsWriteText (sources, name); + + psFree (name); return true; } @@ -381,4 +367,6 @@ bool status; + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + f = fopen (filename, "w"); if (f == NULL) { @@ -387,6 +375,4 @@ } - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); - // write sources with models first for (i = 0; i < sources->n; i++) { @@ -394,7 +380,4 @@ if (source->type != PM_SOURCE_STAR) continue; - // if (source->mode & PM_SOURCE_FAIL) continue; - // if (source->mode & PM_SOURCE_POOR) continue; - model = pmSourceMagnitudes (source, psf, RADIUS); if (model == NULL) continue; @@ -407,5 +390,5 @@ dMag = dPAR[1] / PAR[1]; - fprintf (f, "%7.1f %7.1f %5.1f %8.4f %7.4f %7.4f ", + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); @@ -449,7 +432,6 @@ if (source->type != PM_SOURCE_GALAXY) continue; - - if (source->modelFLT == NULL) continue; - model = pmSourceMagnitudes (source, NULL, source->modelFLT->radius); + model = pmSourceMagnitudes (source, NULL, 0.0); + if (model == NULL) continue; PAR = model->params->data.F32; @@ -458,12 +440,8 @@ // dPos is shape error // XXX these are hardwired for SGAUSS - dPos = 0; - dPos += PS_SQR(dPAR[4] / PAR[4]); - dPos += PS_SQR(dPAR[5] / PAR[5]); - // dPos += PS_SQR(dPAR[7] / PAR[7]); - dPos = sqrt (dPos); + dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); dMag = dPAR[1] / PAR[1]; - fprintf (f, "%7.1f %7.1f %7.1f %7.3f %7.4f %7.4f ", + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); @@ -560,4 +538,26 @@ } +// write the peaks to an output file +bool pmPeaksWriteText (psArray *peaks, char *filename) { + + int i; + FILE *f; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + fprintf (f, "%5d %5d %7.1f\n", + peak->x, peak->y, peak->counts); + } + fclose (f); + return true; +} + // translations between psphot object types and dophot object types int pmSourceDophotType (pmSource *source) { Index: /trunk/psphot/src/psphotTest.c =================================================================== --- /trunk/psphot/src/psphotTest.c (revision 5836) +++ /trunk/psphot/src/psphotTest.c (revision 5837) @@ -3,32 +3,34 @@ int main (int argc, char **argv) { - bool status; + char line[80]; + psImage *image; - psMetadata *config = psphotTestArguments (&argc, argv); - eamReadout *imdata = psphotSetup (config); - - int x = psMetadataLookupF32 (&status, config, "TEST_FIT_X"); - int y = psMetadataLookupF32 (&status, config, "TEST_FIT_Y"); - float fraction = psMetadataLookupF32 (&status, config, "TEST_FRACTION"); - float threshold = psMetadataLookupF32 (&status, config, "TEST_THRESHOLD"); - - if (threshold == -1) { - threshold = fraction * imdata->image->data.F32[y][x]; + image = psImageAlloc (512, 512, PS_DATA_F32); + for (int j = 0; j < 512; j++) { + for (int i = 0; i < 512; i++) { + image->data.F32[j][i] = i; + } } - psArray *contour = pmSourceContour_EAM (imdata->image, x, y, threshold); - if (contour == NULL) { - fprintf (stderr, "error: invalid contour\n"); - exit (1); + psMetadata *row = NULL; + + psArray *table = psArrayAlloc (100); + table->n = 0; + + for (int i = 0; i < 100; i++) { + row = psMetadataAlloc (); + sprintf (line, "line %03d", i); + psMetadataAdd (row, PS_LIST_HEAD, "FLOAT_DATA", PS_DATA_F32, "", 0.2*i); + psMetadataAdd (row, PS_LIST_HEAD, "INT_DATA", PS_DATA_S32, "", i); + psMetadataAdd (row, PS_LIST_HEAD, "CHAR_DATA", PS_DATA_STRING, "", line); + psArrayAdd (table, 100, row); + // psFree (row); } - psVector *xV = contour->data[0]; - psVector *yV = contour->data[1]; + psFits *fits = psFitsOpen ("test.fits", "w"); + psFitsWriteImage (fits, NULL, image, 0); + psFitsWriteTable (fits, NULL, table); + psFitsClose (fits); - for (int i = 0; i < xV->n; i++) { - fprintf (stdout, "%d %f %f\n", i, xV->data.F32[i], yV->data.F32[i]); - } - - // psSparseMatrixTest (); exit (0); }