Changeset 27838 for branches/tap_branches/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (modified) (20 diffs)

branches/tap_branches/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

@@ -45 +45 @@
 // followed by a zero-size matrix, followed by the table data
 
-bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources,
-                                psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
+bool pmSourcesWrite_CMF_PS1_V2 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
@@ -54 +53 @@
     psArray *table;
     psMetadata *row;
-    int i;
     psF32 *PAR, *dPAR;
     psEllipseAxes axes;
@@ -69 +67 @@
     float magOffset = NAN;
     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
+    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
+    if (!isfinite(zeropt)) {
+        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
+    }
     if (status1 && status2 && (exptime > 0.0)) {
         magOffset = zeropt + 2.5*log10(exptime);
     }
+    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -80 +81 @@
         pmSource *source = (pmSource *) sources->data[0];
         if (source->seq == -1) {
-          // let's write these out in S/N order
-          sources = psArraySort (sources, pmSourceSortBySN);
+            // let's write these out in S/N order
+            sources = psArraySort (sources, pmSourceSortBySN);
         } else {
-          sources = psArraySort (sources, pmSourceSortBySeq);
+            sources = psArraySort (sources, pmSourceSortBySeq);
         }
     }
@@ -91 +92 @@
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
-    for (i = 0; i < sources->n; i++) {
+    for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
 
@@ -99 +100 @@
         // generated on Alloc, and would thus be wrong for read in sources.
         if (source->seq == -1) {
-          source->seq = i;
+            source->seq = i;
         }
 
@@ -111 +112 @@
             yPos = PAR[PM_PAR_YPOS];
             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-              xErr = dPAR[PM_PAR_XPOS];
-              yErr = dPAR[PM_PAR_YPOS];
+                xErr = dPAR[PM_PAR_XPOS];
+                yErr = dPAR[PM_PAR_YPOS];
             } else {
-              // in linear-fit mode, there is no error on the centroid
-              xErr = source->peak->dx;
-              yErr = source->peak->dy;
+                // in linear-fit mode, there is no error on the centroid
+                xErr = source->peak->dx;
+                yErr = source->peak->dy;
             }
             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
@@ -209 +210 @@
     }
 
+    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
     psMetadata *header = psMetadataCopy(NULL, tableHeader);
     pmSourceMasksHeader(header);
 
     if (table->n == 0) {
-        psFitsWriteBlank(fits, header, extname);
+        if (!psFitsWriteBlank(fits, header, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write blank sources file.");
+            psFree(table);
+            psFree(header);
+            return false;
+        }
         psFree(table);
         psFree(header);
@@ -221 +228 @@
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable(fits, header, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
         psFree(table);
         psFree(header);
@@ -261 +268 @@
     for (int i = 0; i < numSources; i++) {
         psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i);
+            psFree(sources);
+            return NULL;
+        }
 
         pmSource *source = pmSourceAlloc ();
@@ -304 +316 @@
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
+        source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
 
         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -330 +343 @@
 }
 
-// XXX this layout is still the same as PS1_DEV_1
-bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
+bool pmSourcesWrite_CMF_PS1_V2_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
 {
 
@@ -340 +352 @@
     psF32 xPos, yPos;
     psF32 xErr, yErr;
+    int nRow = -1;
 
     // create a header to hold the output data
@@ -347 +360 @@
     psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
 
+    pmChip *chip = readout->parent->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    // zero point corrections
+    bool status1 = false;
+    bool status2 = false;
+    float magOffset = 0.0;
+    float exptime   = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE");
+    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
+    if (!isfinite(zeropt)) {
+        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
+    }
+    if (status1 && status2 && (exptime > 0.0)) {
+        magOffset = zeropt + 2.5*log10(exptime);
+    }
+
     // let's write these out in S/N order
     sources = psArraySort (sources, pmSourceSortBySN);
@@ -353 +382 @@
 
     // which extended source analyses should we perform?
-    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
-    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
     // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
 
-    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
-    psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    assert (radialBinsLower->n == radialBinsUpper->n);
+    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
+
+    // int nRadialBins = 0;
+    // if (doAnnuli) {
+    //  // get the max count of radial bins
+    //  for (int i = 0; i < sources->n; i++) {
+    //      pmSource *source = sources->data[i];
+    //      if (!source->extpars) continue;
+    //      if (!source->extpars->radProfile ) continue;
+    //         if (!source->extpars->radProfile->binSB) continue;
+    //      nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
+    //  }
+    // }
 
     // we write out all sources, regardless of quality.  the source flags tell us the state
@@ -396 +437 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
 
-# if (0)
+        float AxialRatio = NAN;
+        float AxialTheta = NAN;
+        pmSourceExtendedPars *extpars = source->extpars;
+        if (extpars) {
+            AxialRatio = extpars->axes.minor / extpars->axes.major;
+            AxialTheta = extpars->axes.theta;
+        }
+        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
+        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
+
         // Petrosian measurements
         // XXX insert header data: petrosian ref radius, flux ratio
+        // XXX check flags to see if Pet was measured
         if (doPetrosian) {
-            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
-            if (petrosian) {
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
+            pmSourceExtendedPars *extpars = source->extpars;
+            if (extpars) {
+                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point 
+                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
+                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
             } else {
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
@@ -411 +469 @@
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); 
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); 
             }
         }
 
+# if (0)
         // Kron measurements
         if (doKron) {
@@ -446 +509 @@
             }
         }
-
-        // Flux Annuli
+# endif
+
+        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
         if (doAnnuli) {
-            pmSourceAnnuli *annuli = source->extpars->annuli;
-            if (annuli) {
-                psVector *fluxVal = annuli->flux;
-                psVector *fluxErr = annuli->fluxErr;
-                psVector *fluxVar = annuli->fluxVar;
-
-                for (int j = 0; j < fluxVal->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
-                }
-            } else {
-                for (int j = 0; j < radialBinsLower->n; j++) {
-                    char name[32];
-                    sprintf (name, "FLUX_VAL_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
-                    sprintf (name, "FLUX_ERR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
-                    sprintf (name, "FLUX_VAR_R_%02d", j);
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
-                }
-            }
-        }
-
-# endif
-        psArrayAdd (table, 100, row);
-        psFree (row);
-    }
-
+            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVectorInit (radSB, NAN);
+            psVectorInit (radFlux, NAN);
+            psVectorInit (radFill, NAN);
+            if (!source->extpars) goto empty_annuli;
+            if (!source->extpars->radProfile) goto empty_annuli;
+            if (!source->extpars->radProfile->binSB) goto empty_annuli;
+            psAssert (source->extpars->radProfile->binSum, "programming error");
+            psAssert (source->extpars->radProfile->binFill, "programming error");
+            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
+
+            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
+            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
+                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
+                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
+                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
+            }
+
+        empty_annuli:
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
+            psFree (radSB);
+            psFree (radFlux);
+            psFree (radFill);
+        }
+        if (nRow < 0) {
+            nRow = row->list->n;
+        } else {
+            psAssert (nRow == row->list->n, "inconsistent row lengths");
+        }
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+    
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
-        psFree (outhead);
-        psFree (table);
-        return true;
-    }
-
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+    
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree (outhead);
+    psFree(table);
+    return false;
     }
     psFree (outhead);
     psFree (table);
-
+    
     return true;
 }
 
 // XXX this layout is still the same as PS1_DEV_1
-bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname)
+bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, pmReadout *readout, psArray *sources, char *extname)
 {
 
@@ -609 +684 @@
 
     if (table->n == 0) {
-        psFitsWriteBlank (fits, outhead, extname);
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
         psFree (outhead);
         psFree (table);
@@ -617 +697 @@
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
         psFree (outhead);
         psFree(table);