Changeset 5773

Timestamp:

Dec 13, 2005, 10:03:27 AM (20 years ago)

Author:

eugene

Message:

working updates

Location:

trunk/psphot/src

Files:

• psphot.c (modified) (2 diffs)
• psphotFullFit.c (modified) (2 diffs)
• psphotMagnitudes.c (modified) (2 diffs)
• psphotOutput.c (modified) (9 diffs)

trunk/psphot/src/psphot.c

@@ -75 +75 @@
       case 0:
         psphotEnsemblePSF (imdata, config, sources, psf, sky);
-        sources = psArraySort (sources, psphotSortBySN);
         psphotFullFit (imdata, config, sources, psf, sky);
 //      psphotReplaceUnfit (sources);
@@ -99 +98 @@
 
     // write out data in appropriate format
+    // measure aperture correction
+    // psphotApertureCorrection (sources, config, sky);
     psphotOutput (imdata, config, sources, psf, sky);
     exit (0);

trunk/psphot/src/psphotFullFit.c

@@ -79 +79 @@
             continue;
         }
+
+        // add the double-PSF fit mode
+        // how do we compare the results?  chisq? 
         
         // recalculate the source moments using the larger galaxy moments radius
@@ -113 +116 @@
         }
 
-        // subtract object, leave local sky
+        // subtract PSF for object, leave local sky
+        // XXX : we should keep the non-linear PSF fit if object was marked extended but
+        //       PSF fit is OK otherwise.
         psFree (source->modelPSF);
         source->modelPSF = save;

trunk/psphot/src/psphotMagnitudes.c

@@ -1 +1 @@
 # include "psphot.h"
 
-// XXX EAM : this aperture correction business is invalid (& wrong) for galaxies
+// XXX EAM : the apMag should only be calculated for the brighter sources
+// XXX EAM : SN limit set by user?
+// XXX EAM : masked region should be (optionally) elliptical
 pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) {
 
     int status;
     float x, y;
-    float sky, rflux, apMag, fitMag;
+    float rflux, apMag, fitMag;
+    pmModel *model;
 
     // use the correct model (PSF vs FLT)
-    pmModel *model = pmSourceSelectModel (source);
+    if (psf != NULL) {
+      model = source->modelPSF;
+    } else {
+      model = source->modelFLT;
+    }
     if (model == NULL) return NULL;
 
@@ -20 +27 @@
     psImageKeepCircle (source->mask, x, y, apRadius, "OR", PSPHOT_MASK_MARKED);
 
-    // save local sky for later
-    sky = model->params->data.F32[0];
-
     // replace source flux
     pmSourceAddModel (source->pixels, source->mask, model, false, false);
 
     // measure object photometry
-    status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+    status = pmSourcePhotometry (&source->fitMag, &source->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);
-    source->apMag  = apMag  - rflux * psf->skyBias * (M_PI * PS_SQR(apRadius));
-    source->fitMag = fitMag + psf->ApResid;
+    // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
+    if (psf != NULL) {
+      rflux   = pow (10.0, 0.4*source->fitMag);
+      source->apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius));
+      source->fitMag += psf->ApResid;
+    }
 
     // subtract object, leave local sky

trunk/psphot/src/psphotOutput.c

@@ -367 +367 @@
 }
 
-// dump the sources to an output file
+// write the PSF sources to an output file
 bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf) {
 
@@ -389 +389 @@
         pmSource *source = (pmSource *) sources->data[i];
 
-        // write out sources fitted as PSFs
         if (source->type != PM_SOURCE_STAR) continue;
-        if (source->mode & PM_SOURCE_FAIL) continue;
-        if (source->mode & PM_SOURCE_POOR) continue;
-
-        model = pmSourceMagnitudes (source, psf, RADIUS);
+        // if (source->mode & PM_SOURCE_FAIL) continue;
+        // if (source->mode & PM_SOURCE_POOR) continue;
+
+        model = pmSourceMagnitudes (source, psf, RADIUS, true);
         if (model == NULL) continue;
 
@@ -400 +399 @@
         dPAR = model->dparams->data.F32;
         
+        // dP is positional error
         dP = 0;
         dP += PS_SQR(dPAR[2]);
@@ -410 +410 @@
                  PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
 
-        for (j = 0; j < model->params->n - 4; j++) {
-            fprintf (f, "%9.6f ", PAR[j+4]);
+        for (j = 4; j < model->params->n; j++) {
+            fprintf (f, "%9.6f ", PAR[j]);
         }
         fprintf (f, " : ");
-        for (j = 0; j < model->params->n - 4; j++) {
-            fprintf (f, "%9.6f ", dPAR[j+4]);
+        for (j = 4; j < model->params->n; j++) {
+            fprintf (f, "%9.6f ", dPAR[j]);
         }
-        fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
-                 source[0].type, 
+        fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+                 source[0].type, source[0].mode, 
                  log10(model[0].chisq/model[0].nDOF), 
                  source[0].moments->SN, 
@@ -448 +448 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        model = (pmModel  *) source->modelFLT;
+
+        if (source->type != PM_SOURCE_GALAXY) continue;
+        // if (source->mode & PM_SOURCE_FAIL) continue;
+        // if (source->mode & PM_SOURCE_POOR) continue;
+        
+        model = pmSourceMagnitudes (source->modelFLT, NULL, model->radius);
         if (model == NULL) continue;
-        if (source->type != PM_SOURCE_GALAXY) continue;
-        if (source->mode & PM_SOURCE_FAIL) continue;
-        if (source->mode & PM_SOURCE_POOR) continue;
+
+        PAR  = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
         
-        params = model->params;
-        dparams = model->dparams;
-
-        // XXX these are hardwired for SGAUSS : this should be pushed into the
-        // model functions as an abstract function 
+        // dP is shape error
+        // XXX these are hardwired for SGAUSS
         dP = 0;
-        dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]);
-        dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]);
-        dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]);
+        dP += PS_SQR(dPAR[4] / PAR[4]);
+        dP += PS_SQR(dPAR[5] / PAR[5]);
+        dP += PS_SQR(dPAR[7] / PAR[7]);
         dP = sqrt (dP);
 
+        dmag = dPAR[1] / PAR[1];
+
         fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
-                 params[0].data.F32[2], 
-                 params[0].data.F32[3], 
-                 params[0].data.F32[0], 
-                 -2.5*log10(params[0].data.F32[1]), 
-                 (dparams[0].data.F32[1]/params[0].data.F32[1]),
-                 dP);
+                 PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
 
         for (j = 4; j < model->params->n; j++) {
-            fprintf (f, "%9.6f ", params[0].data.F32[j]);
+            fprintf (f, "%9.6f ", PAR[j]);
         }
         fprintf (f, " : ");
         for (j = 4; j < model->params->n; j++) {
-            fprintf (f, "%9.6f ", dparams[0].data.F32[j]);
+            fprintf (f, "%9.6f ", dPAR[j]);
         }
-        fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
-                 source[0].type, log10(model[0].chisq/model[0].nDOF), 
+        fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+                 source[0].type, source[0].mode,
+                 log10(model[0].chisq/model[0].nDOF), 
                  source[0].moments->SN, 
                  model[0].radius, 
@@ -497 +497 @@
     int i;
     FILE *f;
+    pmMoments *moment = NULL;
+    pmSource *source = NULL;
 
     f = fopen (filename, "w");
@@ -508 +510 @@
     // write sources with models first 
     for (i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
+        source = sources->data[i];
 
         // skip these sources (in PSF or FLT)
-        if (source->type == PM_SOURCE_STAR) {
-            if (source->mode & PM_SOURCE_FAIL) goto isNULL;
-            if (source->mode & PM_SOURCE_POOR) goto isNULL;
-            continue;
+        if (source->type == PM_SOURCE_STAR) continue;
+        if (source->type == PM_SOURCE_GALAXY) continue;
+            
+        if (source->moments == NULL) {
+          moment = empty;
+        } else {
+          moment = source->moments;
         }
-        if (source->type == PM_SOURCE_GALAXY) {
-            if (source->mode & PM_SOURCE_FAIL) goto isNULL;
-            if (source->mode & PM_SOURCE_POOR) goto isNULL;
-            continue;
-        }
-            
-    isNULL:
-        
-        if (source->moments == NULL) {
-          fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
-                   source->peak->x, source->peak->y, source->peak->counts, 
-                   empty->x, empty->y, 
-                   empty->Sx, empty->Sy, 
-                   empty->Sum, empty->Peak, 
-                   empty->Sky, empty->SN, 
-                   empty->nPixels, source->type); 
-        } else {
-          fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
-                   source->peak->x, source->peak->y, source->peak->counts, 
-                   source->moments->x, source->moments->y, 
-                   source->moments->Sx, source->moments->Sy, 
-                   source->moments->Sum, source->moments->Peak, 
-                   source->moments->Sky, source->moments->SN, 
-                   source->moments->nPixels, source->type); 
-        }
+
+        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+                 source->peak->x, source->peak->y, source->peak->counts, 
+                 source->moments->x, source->moments->y, 
+                 source->moments->Sx, source->moments->Sy, 
+                 source->moments->Sum, source->moments->Peak, 
+                 source->moments->Sky, source->moments->SN, 
+                 source->moments->nPixels, source->type); 
     }
     fclose (f);
@@ -552 +540 @@
     int i;
     FILE *f;
-    pmSource *source;
+    pmSource *source = NULL;
 
     f = fopen (filename, "w");
@@ -561 +549 @@
 
     for (i = 0; i < sources->n; i++) {
-        source = (pmSource *) sources->data[i];
+        source = sources->data[i];
         if (source->moments == NULL)
             continue;