Changeset 5073 for trunk/psphot/src/psphotOutput.c

Timestamp:

Sep 19, 2005, 5:24:40 PM (21 years ago)

Author:

eugene

Message:

added CMP, CMF output formats

File:

• trunk/psphot/src/psphotOutput.c (modified) (12 diffs)

trunk/psphot/src/psphotOutput.c

@@ -73 +73 @@
     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 +117 @@
         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 +132 @@
         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 +146 @@
         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 +193 @@
 
     // 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 +272 @@
         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 +317 @@
     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 +353 @@
         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 +372 @@
         // 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 +383 @@
         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 +401 @@
         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 +412 @@
 
 // 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;