Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5072) +++ /trunk/psphot/src/psphotOutput.c (revision 5073) @@ -73,5 +73,5 @@ pmModel *model; psF32 *PAR, *dPAR; - float skyFlux, dmag, apMag, fitMag; + float sky, dmag, apMag, fitMag; psLine *line = psLineAlloc (104); // 104 is dophot-defined line length @@ -117,6 +117,14 @@ dPAR = model->dparams->data.F32; + 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, RADIUS, "OR", PSPHOT_MASK_MARKED); + // save local sky for later - skyFlux = model->params->data.F32[0]; + sky = model->params->data.F32[0]; // replace source flux @@ -124,8 +132,12 @@ model->params->data.F32[0] = 0; pmSourceAddModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = skyFlux; + model->params->data.F32[0] = sky; // measure object photometry status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + rflux = pow (10.0, 0.4*fitMag); + apMag -= rflux * psf->skyBias * (M_PI * PS_SQR(RADIUS)); fitMag += psf->ApResid; @@ -134,5 +146,8 @@ model->params->data.F32[0] = 0; pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = skyFlux; + model->params->data.F32[0] = sky; + + // unmask aperture + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); if (status == FALSE) continue; @@ -178,4 +193,6 @@ // aperture radius for ap magnitude + // XXX EAM I probably need to measure the fit and ap magnitudes at one location elsewhere + // this would have to happen *after* all objects are detected and subtracted float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); @@ -255,5 +272,5 @@ if (status == FALSE) continue; - // XXX EAM : fix this to match sextractor types? + // XXX EAM : fix this to return sextractor type and flags // type = pmSourceSextractType (source); @@ -300,5 +317,5 @@ fseek (f, 0, SEEK_END); - psLine *line = psLineAlloc (104); // 104 is dophot-defined line length + psLine *line = psLineAlloc (66); // 66 is imclean-defined line length // write sources with models first @@ -336,9 +353,14 @@ dPAR = model->dparams->data.F32; + 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, RADIUS, "OR", PSPHOT_MASK_MARKED); + // save local sky for later sky = model->params->data.F32[0]; - - // XXX EAM : get Theta from PAR[4,5,6] - Theta = 0.0; // replace source flux @@ -350,4 +372,8 @@ // measure object photometry status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + rflux = pow (10.0, 0.4*fitMag); + apMag -= rflux * psf->skyBias * (M_PI * PS_SQR(RADIUS)); fitMag += psf->ApResid; @@ -357,4 +383,7 @@ pmSourceSubModel (source->pixels, source->mask, model, false); model->params->data.F32[0] = sky; + + // unmask aperture + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); if (status == FALSE) continue; @@ -372,5 +401,5 @@ psLineAdd (line, "%3.1f ", lsky); psLineAdd (line, "%6.3f ", 32.0); // should be 'Mgal - psLineAdd (line, "%6.3f ", apMag + 25.0); + psLineAdd (line, "%6.3f ", apMag + 25.0); // 25.0 should come from config psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' @@ -383,9 +412,117 @@ // elixir-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data bool pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { - // write imdata->header to file - - // write FITS table to file (use autocode tools) + psArray *table; + psMetadata *row; + psMetadata *theader; + + // set NAXIS to 0 + psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0); + + table = psArrayAlloc (sources->n); + sources->n = 0; + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // use the correct model + switch (source->type) { + + // use PSF model with stars + case PM_SOURCE_PSFSTAR: + case PM_SOURCE_SATSTAR: + case PM_SOURCE_GOODSTAR: + 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 + case PM_SOURCE_DEFECT: + case PM_SOURCE_SATURATED: + case PM_SOURCE_FAINTSTAR: + case PM_SOURCE_POOR_FIT_PSF: + case PM_SOURCE_POOR_FIT_GAL: + case PM_SOURCE_FAIL_FIT_GAL: + default: + model = NULL; + } + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + 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, RADIUS, "OR", PSPHOT_MASK_MARKED); + + // save local sky for later + sky = model->params->data.F32[0]; + + // replace source flux + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceAddModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // measure object photometry + status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS + rflux = pow (10.0, 0.4*fitMag); + apMag -= rflux * psf->skyBias * (M_PI * PS_SQR(RADIUS)); + fitMag += psf->ApResid; + + // subtract object, leave local sky + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceSubModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // unmask aperture + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + if (status == FALSE) continue; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + table = psMetdataAlloc (); + psMetadataAdd (table, PS_LIST_HEAD, "X_PIX", PS_META_F32, "", PAR[2]); + psMetadataAdd (table, PS_LIST_HEAD, "Y_PIX", PS_META_F32, "", PAR[3]); + psMetadataAdd (table, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0); + psMetadataAdd (table, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag)); + psMetadataAdd (table, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type); + psMetadataAdd (table, PS_LIST_HEAD, "MAG_AP", PS_META_F32, "", lsky); + psMetadataAdd (table, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0); // should be 'Mgal' + psMetadataAdd (table, PS_LIST_HEAD, "FWHM_X", PS_META_F32, "", apMag + 25.0); // 25.0 should come from config + psMetadataAdd (table, PS_LIST_HEAD, "FWHM_Y", PS_META_F32, "", PAR[4]); // should be 'FHWM x' + psMetadataAdd (table, PS_LIST_HEAD, "THETA", PS_META_F32, "", PAR[5]); // should be 'FHWM y' + psMetadataAdd (table, PS_LIST_HEAD, "DOPHOT", PS_META_S8, "", Theta); // should be theta + // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY", PS_META_STR, "", NULL); + + psArrayAdd (table, 100, row); + psFree (row); + } + + image = psImageAlloc (0, 0, PS_TYPE_F32); + theader = psFitsHeaderFromTable (NULL, table, "SMPFILE"); + + unlink (filename); + psFits *fits = psFitsAlloc (filename); + psFitsWriteImage (fits, header, image, 0, NULL); + psFitsWriteTable (fits, theader, table); + psFree (fits); return true;