Changeset 15400

Timestamp:

Oct 28, 2007, 3:41:13 PM (19 years ago)

Author:

eugene

Message:

adding extended source measurements

Location:

branches/eam_branch_20071023/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotAnnuli.c (modified) (2 diffs)
• psphotIsophotal.c (modified) (4 diffs)
• psphotMakeResiduals.c (modified) (3 diffs)
• psphotPetrosian.c (modified) (5 diffs)
• psphotRadialProfile.c (modified) (3 diffs)

branches/eam_branch_20071023/psphot/src/Makefile.am

@@ -47 +47 @@
         psphotModelWithPSF.c     \
         psphotExtendedSources.c  \
+        psphotRadialProfile.c    \
         psphotPetrosian.c        \
         psphotIsophotal.c        \

branches/eam_branch_20071023/psphot/src/psphot.h

@@ -116 +116 @@
 psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
 
+bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool            psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal);

branches/eam_branch_20071023/psphot/src/psphotAnnuli.c

@@ -8 +8 @@
   assert (source->extpars->profile->flux);
 
+  bool status;
+
   psVector *radius = source->extpars->profile->radius;
   psVector *weight = source->extpars->profile->weight;
@@ -14 +16 @@
   // XXX how do I define the radii?  we can put a vector in the recipe...
   // radialBins defines the bounds or start and stop (we can skip some that way...
-  psVector *radialBinsLower = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_LOWER);
-  psVector *radialBinsUpper = psMetadataLookupVector (&status, recipe, RADIAL_ANNULAR_BINS_UPPER);
+  psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+  psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
   assert (radialBinsLower->n == radialBinsUpper->n);
 

branches/eam_branch_20071023/psphot/src/psphotIsophotal.c

@@ -7 +7 @@
   assert (source->extpars->profile->radius);
   assert (source->extpars->profile->flux);
+
+  bool status;
 
   psVector *radius = source->extpars->profile->radius;
@@ -19 +21 @@
   // XXX do I need to worry about crazy outliers?  
   // XXX should i be smoothing or fitting the curve?
-  firstBelow = -1;
-  lastAbove = -1;
+  int firstBelow = -1;
+  int lastAbove = -1;
   for (int i = 0; i < flux->n; i++) {
     if (flux->data.F32[i] > ISOPHOT_FLUX) lastAbove = i;
@@ -31 +33 @@
   for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
     if (i <= firstBelow) {
-      isophotFluxFirst += flux->data.F32[i];
+      isophotalFluxFirst += flux->data.F32[i];
     }
     if (i <= lastAbove) {
-      isophotFluxLast += flux->data.F32[i];
+      isophotalFluxLast += flux->data.F32[i];
     }
   }
-  float isophotalFlux    = 0.5*(isophotFluxLast + isophotFluxFirst);
-  float isophotalFluxErr = 0.5*fabs(isophotFluxLast - isophotFluxFirst);
+  float isophotalFlux    = 0.5*(isophotalFluxLast + isophotalFluxFirst);
+  float isophotalFluxErr = 0.5*fabs(isophotalFluxLast - isophotalFluxFirst);
 
   float isophotalRad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
@@ -57 +59 @@
 
 }
-
-# if (0)
-  // XXX cache the tmpScalar 
-  psScalar fluxScalar;
-  fluxScalar.type.type = PS_TYPE_F32;
-  fluxScalar.data.F32  = ISOPHOT_FLUX;
-
-  // radius and flux must be sorted by radius
-
-  // XXX in general, flux decreases monotonically with radius.  however, since exceptions are
-  // possible we need to do something to smooth or otherwise handle the flux variations
-  
-  // get the index of the flux-sorted vector
-  psVector *fluxIndex = psVectorSortIndex (flux);
-
-  // XXX need to write psVectorIndexBinaryDisect, which operates on a 
-
-  // use disection to get in the right vicinity
-  binS = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_FIRST);
-  binE = psVectorIndexBinaryDisect (&result, flux, fluxIndex, fluxScalar, PS_BISECT_LAST);
-
-  // find the min and max radius bins in this range
-
-  // XXX do something to choose a more accurate radial bin
-
-# endif

branches/eam_branch_20071023/psphot/src/psphotMakeResiduals.c

@@ -28 +28 @@
 
     float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
+    PS_ASSERT (status, false);
+
+    float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
     PS_ASSERT (status, false);
 
@@ -195 +198 @@
                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
 
-                // if value < NRESID_SIGMA*sigma, mask pixel in resid map
-                if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*fluxStats->sampleStdev) {
-                  resid->mask->data.U8[oy][ox] = RESID_SIG_MASK;
+                if (resid->Ro->data.F32[oy][ox] < pixelSN*fluxStats->sampleStdev) {
+                  resid->mask->data.U8[oy][ox] = 1;
                 }
 
@@ -234 +236 @@
                 resid->Ry->data.F32[oy][ox] = B->data.F64[2];
 
-                dRo = sqrt(A->data.F32[0][0]);
-                if (resid->Ro->data.F32[oy][ox] < NRESID_SIGMA*dRo) {
-                  resid->mask->data.U8[oy][ox] = RESID_SIG_MASK;
+                float dRo = sqrt(A->data.F32[0][0]);
+                if (resid->Ro->data.F32[oy][ox] < pixelSN*dRo) {
+                  resid->mask->data.U8[oy][ox] = 1;
                 }
                 //resid->weight->data.F32[oy][ox] = XXX;

branches/eam_branch_20071023/psphot/src/psphotPetrosian.c

@@ -2 +2 @@
 
 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+
+  bool status;
 
   assert (source->extpars);
@@ -17 +19 @@
 
   // first find flux at R0
-  firstBelow = -1;
-  lastAbove = -1;
+  int firstBelow = -1;
+  int lastAbove = -1;
   for (int i = 0; i < radius->n; i++) {
     if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
@@ -25 +27 @@
 
   // average flux in this range
-  fluxR0 = 0.0;
-  fluxRn = 0;
+  float fluxR0 = 0.0;
+  int fluxRn = 0;
   for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {
     fluxR0 += flux->data.F32[i];
     fluxRn ++;
   }
-  fluxR0 = fluxR0 / (float)(fluxRn);
+  fluxR0 /= (float)(fluxRn);
   
   // target flux for petrosian radius
-  fluxRP = fluxR0 * PETROSIAN_RF;
+  float fluxRP = fluxR0 * PETROSIAN_RF;
 
   // find the first bin below the flux level and the last above the level
@@ -58 +60 @@
     }
   }
-  float flux    = 0.5*(fluxLast + fluxFirst);
-  float fluxErr = 0.5*fabs(fluxLast - fluxFirst);
+  float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
+  float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
   // XXX need to use the weight appropriately here...
 
@@ -70 +72 @@
   
   // these are uncalibrated: instrumental mags and pixel units
-  source->extpars->petrosian->mag    = -2.5*log10(flux);
-  source->extpars->petrosian->magErr = fluxErr / flux;
+  source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
+  source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;
 
   source->extpars->petrosian->rad    = rad;

branches/eam_branch_20071023/psphot/src/psphotRadialProfile.c

@@ -5 +5 @@
     // allocate pmSourceExtendedParameters, if not already defined
     if (!source->extpars) {
-        source->extpars = pmSourceExtendedParametersAlloc ();
+        source->extpars = pmSourceExtendedParsAlloc ();
     }
 
     if (!source->extpars->profile) {
-        source->extpars->profile = pmSourceRadialProfile (); 
+        source->extpars->profile = pmSourceRadialProfileAlloc (); 
     }    
     
@@ -19 +19 @@
     psVector *radius = source->extpars->profile->radius;
     psVector *flux   = source->extpars->profile->flux;
-    psVector *weight = source->extpars->profile->radius;
+    psVector *weight = source->extpars->profile->weight;
 
     // XXX use the extended source model here for Xo, Yo?
@@ -40 +40 @@
     flux->n = n;
 
-    SortVectorsByRadius (radius, flux, weight);
+    // XXX need to sort here
+    // SortVectorsByRadius (radius, flux, weight);
 
     return true;