Changeset 36148 for branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c

Timestamp:

Sep 24, 2013, 12:09:27 PM (13 years ago)

Author:

eugene

Message:

first version of psphotFullForce which builds

Location:

branches/eam_branches/ipp-20130904/psphot/src

Files:

• . (modified) (1 prop)
• psphotPetroFlux.c (modified) (13 diffs)

branches/eam_branches/ipp-20130904/psphot/src

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce

branches/eam_branches/ipp-20130904/psphot/src/psphotPetroFlux.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
+bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule)
 {
     bool status = true;
@@ -35 +35 @@
         psAssert (sources, "missing sources?");
 
-        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-        // psAssert (psf, "missing psf?");
-
-        if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
+        if (!psphotPetroFluxReadout (config, recipe, view, filerule, readout, sources)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -46 +43 @@
 }
 
-bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
+bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources) {
 
     bool status = false;
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping kron fluxes");
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping petrosian fluxes");
         return true;
     }
 
-    psTimerStart ("psphot.kron");
+    psTimerStart ("psphot.petro");
 
     // determine the number of allowed threads
@@ -61 +58 @@
     if (!status) {
         nThreads = 0;
-    }
-
-    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
-    if (!status) {
-        KRON_APPLY_WEIGHT = true;
-    }
-
-    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
-    if (!status) {
-        KRON_APPLY_WINDOW = false;
     }
 
@@ -105 +92 @@
 
             // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("PSPHOT_KRON_ITERATE");
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_PETRO_FLUX");
 
             psArrayAdd(job->args, 1, readout);
@@ -111 +98 @@
             PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
-            PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
 
 // set this to 0 to run without threading
@@ -146 +132 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot.kron", PS_LOG_WARN, "measure kron fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
+    psLogMsg ("psphot.petro", PS_LOG_WARN, "measure petro fluxes : %f sec for %ld objects\n", psTimerMark ("psphot.petro"), sources->n);
     return true;
 }
@@ -153 +139 @@
 
     pmReadout *readout              = job->args->data[0];
-    psArray *sources                = job->args->data[2];
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+    psArray *sources                = job->args->data[1];
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
 
     for (int i = 0; i < sources->n; i++) {
-
-        // XXX what do i need for the kron flux (given a supplied kron radius?)
 
         pmSource *source = sources->data[i];
@@ -169 +153 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
 
-        // XXX where are we storing the supplied kron radius?
+        // XXX where are we storing the supplied petro radius?
         if (!isfinite(source->moments->Mrf)) continue;
 
@@ -180 +164 @@
         }
 
-        // is the right??
-        float windowRadius = 6.0*source->moments->Mrf ;
+        // is the right??  XXX need the correct location for Rpet
+        float windowRadius = 2.0*source->moments->Mrf ;
 
         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-        bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
         psAssert (source->pixels, "WTF?");
 
@@ -211 +195 @@
     float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
 
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
     // Calculate the Petro magnitude (make this block optional?)
     // XXX set the aperture here
-    float radPetro  = 2.5*Mrf;
+    float radPetro  = 2.0*source->moments->Mrf;
     float radPetro2 = radPetro*radPetro;
 
@@ -222 +211 @@
     // set vPix to the source pixels (it may have been set to the
     // smoothed image above)
-    vPix = source->pixels->data.F32;
-    vWgt = source->variance->data.F32;
-    vMsk = source->mask->data.F32;
+    psF32 **vPix = source->pixels->data.F32;
+    psF32 **vWgt = source->variance->data.F32;
+    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -253 +242 @@
 
     // return the flux and error to parameters?
-    source->moments->PetroFlux    = Sum;
-    source->moments->PetroFluxErr = sqrt(Var);
-
-    return true;
-}
+    // source->moments->PetroFlux    = Sum;
+    // source->moments->PetroFluxErr = sqrt(Var);
+
+    fprintf (stderr, "petro flux: %f +/- %f\n", Sum, sqrt(Var));
+
+    return true;
+}

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce