Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 5084) +++ trunk/psphot/src/psphotOutput.c (revision 5085) @@ -70,10 +70,9 @@ bool pmSourcesWriteOBJ (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { - int i, type, status; + int type; pmModel *model; psF32 *PAR, *dPAR; - float sky, dmag, apMag, fitMag; - float x, y, rflux; - bool result; + float dmag, apResid; + bool status; psLine *line = psLineAlloc (104); // 104 is dophot-defined line length @@ -85,35 +84,11 @@ } - float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); - - // write sources with models first - for (i = 0; i < sources->n; i++) { + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + + // write sources with models + for (int 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; - } + model = pmSourceMagnitudes (source, psf, RADIUS); if (model == NULL) continue; @@ -121,47 +96,13 @@ 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); + apResid = source->apMag - source->fitMag; + psLineInit (line); - psLineAdd (line, "%3d", source->type); + psLineAdd (line, "%3d", type); psLineAdd (line, "%8.2f", PAR[2]); psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", fitMag); + psLineAdd (line, "%8.3f", source->fitMag); psLineAdd (line, "%6.3f", dmag); psLineAdd (line, "%9.2f", PAR[0]); @@ -170,6 +111,6 @@ psLineAdd (line, "%7.2f", PAR[6]); psLineAdd (line, "%8.3f", 32.0); - psLineAdd (line, "%8.3f", apMag); - psLineAdd (line, "%8.2f\n", apMag - fitMag); + psLineAdd (line, "%8.3f", source->apMag); + psLineAdd (line, "%8.2f\n", apResid); fwrite (line->line, 1, line->Nline, f); } @@ -181,10 +122,8 @@ bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) { - int i, status, flags; pmModel *model; psF32 *PAR, *dPAR; - float sky, dmag, apMag, fitMag; - float x, y, rflux; - bool result; + float dmag; + bool status; psLine *line = psLineAlloc (110); // 110 is sextractor line length @@ -192,47 +131,15 @@ FILE *f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - // 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 + psLogMsg ("WriteSourceSX", 3, "can't open output file for output %s\n", filename); + return false; + } + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); - // write sources with models first - for (i = 0; i < sources->n; i++) { + // write sources with models + for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - flags = 0; - - // use the correct model - switch (source->type) { - - // use PSF model with stars - case PM_SOURCE_SATSTAR: - flags = 4; - model = source->modelPSF; - break; - - case PM_SOURCE_PSFSTAR: - 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; - } + + model = pmSourceMagnitudes (source, psf, RADIUS); if (model == NULL) continue; @@ -240,49 +147,13 @@ 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; - - // XXX EAM : fix this to return sextractor type and flags - // type = pmSourceSextractType (source); + // pmSourceSextractType (source, &type, &flags); XXX EAM : implement me... dmag = dPAR[1] / PAR[1]; + psLineInit (line); psLineAdd (line, "%5.2f", 0.0); // should be type psLineAdd (line, "%11.3f", PAR[2]); psLineAdd (line, "%11.3f", PAR[3]); - psLineAdd (line, "%9.4f", fitMag); + psLineAdd (line, "%9.4f", source->fitMag); psLineAdd (line, "%9.4f", dmag); psLineAdd (line, "%13.4f", PAR[0]); @@ -291,5 +162,5 @@ psLineAdd (line, "%6.1f", 0.0); // should be Theta psLineAdd (line, "%9.4f", 32.0); // should be MAG_ISO - psLineAdd (line, "%9.4f", apMag); + psLineAdd (line, "%9.4f", source->apMag); psLineAdd (line, "%4d\n", 0); // should be flags fwrite (line->line, 1, line->Nline, f); @@ -310,4 +181,22 @@ bool result; + // find config information for output header + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&result, config, "ZERO_POINT"); + char *PHOTCODE = psMetadataLookupF32 (&result, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_META_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT", PS_META_F32 | PS_META_REPLACE, "zero point", ZERO_POINT); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT", PS_META_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_META_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS", PS_META_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_META_STR | PS_META_REPLACE, "photometry code", PHOTCODE); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X", PS_META_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_META_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_META_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_META_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_META_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + // create file, write-out header unlink (filename); @@ -333,44 +222,9 @@ psLine *line = psLineAlloc (67); // 66 is imclean-defined line length - float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); - - // XXX EAM : need to add these concepts to the header (from psf and config) - // fits_modify (header, "ZERO_PT", "%lf", 1, ZERO_POINT); - // fits_modify (header, "FWHM_X", "%lf", 1, FWHMx); - // fits_modify (header, "FWHM_Y", "%lf", 1, FWHMy); - // fits_modify (header, "ANGLE", "%lf", 1, angle); - // fits_modify (header, "FSATUR", "%lf", 1, saturate); // can use 0.0 - // fits_modify (header, "FLIMIT", "%lf", 1, complete); // can use 0.0 - // fits_modify (header, "NSTARS", "%d", 1, N); - // write sources with models first 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; - } + model = pmSourceMagnitudes (source, psf, RADIUS); if (model == NULL) continue; @@ -378,42 +232,7 @@ 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); + apResid = source->apMag - source->fitMag; lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); @@ -421,10 +240,10 @@ psLineAdd (line, "%6.1f ", PAR[2]); psLineAdd (line, "%6.1f ", PAR[3]); - psLineAdd (line, "%6.3f ", fitMag + 25.0); + psLineAdd (line, "%6.3f ", source->fitMag + ZERO_POINT); psLineAdd (line, "%03d ", (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 ", apMag + 25.0); // 25.0 should come from config + psLineAdd (line, "%6.3f ", apMag + ZERO_POINT); psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' @@ -452,4 +271,22 @@ bool result; + // find config information for output header + float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); + float ZERO_POINT = psMetadataLookupF32 (&result, config, "ZERO_POINT"); + char *PHOTCODE = psMetadataLookupF32 (&result, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_META_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT", PS_META_F32 | PS_META_REPLACE, "zero point", ZERO_POINT); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT", PS_META_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_META_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS", PS_META_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_META_STR | PS_META_REPLACE, "photometry code", PHOTCODE); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X", PS_META_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_META_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_META_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_META_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_META_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + // set NAXIS to 0 fprintf (stderr, "setting naxis\n"); @@ -458,79 +295,15 @@ mdi->type = PS_META_S32; - // psMetadataAdd (imdata->header, PS_LIST_HEAD, "NAXIS", PS_META_S32 | PS_META_REPLACE, "", 0); - table = psArrayAlloc (sources->n); sources->n = 0; - float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS"); - 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; - } + model = pmSourceMagnitudes (source, psf, RADIUS); 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]; @@ -541,10 +314,10 @@ psMetadataAdd (row, PS_LIST_HEAD, "X_PIX", PS_META_F32, "", PAR[2]); psMetadataAdd (row, PS_LIST_HEAD, "Y_PIX", PS_META_F32, "", PAR[3]); - psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", fitMag + 25.0); + psMetadataAdd (row, PS_LIST_HEAD, "MAG_RAW", PS_META_F32, "", source->fitMag + ZERO_POINT); psMetadataAdd (row, PS_LIST_HEAD, "MAG_ERR", PS_META_F32, "", (int)(1000*dmag)); psMetadataAdd (row, PS_LIST_HEAD, "MAG_GAL", PS_META_F32, "", type); psMetadataAdd (row, PS_LIST_HEAD, "MAG_AP", PS_META_F32, "", lsky); psMetadataAdd (row, PS_LIST_HEAD, "LOG_SKY", PS_META_F32, "", 32.0); // should be 'Mgal' - psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X", PS_META_F32, "", apMag + 25.0); // 25.0 should come from config + psMetadataAdd (row, PS_LIST_HEAD, "FWHM_X", PS_META_F32, "", apMag + ZERO_POINT); psMetadataAdd (row, PS_LIST_HEAD, "FWHM_Y", PS_META_F32, "", PAR[4]); // should be 'FHWM x' psMetadataAdd (row, PS_LIST_HEAD, "THETA", PS_META_F32, "", PAR[5]); // should be 'FHWM y'