Changeset 21106 for branches/pap_branch_20090108/ppStack/src/ppStackSources.c

Timestamp:

Jan 12, 2009, 1:03:29 PM (17 years ago)

Author:

Paul Price

Message:

Adding more debugging output. Think I've solved the problem of the poor photometry: the background being added in to the fake images produces the bad photometry. Without it, the photometry s.d. goes from 0.08 mag to 0.02 mag. Want to replace it with the proper psphot background model (not done yet). These files are currently in an intermediate state.

File:

• branches/pap_branch_20090108/ppStack/src/ppStackSources.c (modified) (6 diffs)

branches/pap_branch_20090108/ppStack/src/ppStackSources.c

@@ -9 +9 @@
 //#define TESTING                         // Enable debugging output
 
+// Size of fake image; set by hand because it's trouble to get it from other places
+#define FAKE_COLS 4861
+#define FAKE_ROWS 4913
+
+
+
+#define MIN_MAG -16                     // Minimum magnitude to consider
 
 float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
@@ -59 +66 @@
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
+
+#if 1
+        // Excise bright sources that might be saturated
+        psArray *sources = sourceLists->data[i]; // Sources for image
+        for (int j = 0; j < sources->n; j++) {
+            pmSource *source = sources->data[j]; // Source of interest
+            if (source && source->psfMag < MIN_MAG) {
+                psFree(source);
+                sources->data[j] = NULL;
+            }
+        }
+#endif
+
+
+#if 1
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "start_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
 
         float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
@@ -105 +139 @@
         FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
         psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
+        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
         for (int i = 0; i < matches->n; i++) {
             pmSourceMatch *match = matches->data[i]; // Match of interest
             psVectorInit(mag, NAN);
+            psVectorInit(magErr, NAN);
             for (int j = 0; j < match->num; j++) {
                 if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
@@ -114 +150 @@
                 int index = match->image->data.U32[j]; // Image index
                 mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
+                magErr->data.F32[index] = match->magErr->data.F32[j];
             }
             for (int j = 0; j < num; j++) {
-                fprintf(outMatches, "%f ", mag->data.F32[j]);
+                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
             }
             fprintf(outMatches, "\n");
         }
         psFree(mag);
+        psFree(magErr);
         fclose(outMatches);
     }
@@ -146 +184 @@
                 continue;
             }
+
+#if 1
+            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS],
+                y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; // Coordinates of source
+            if (x > 2187 && x < 2191 && y > 2280 && y < 2284) {
+                fprintf(stderr, "Image %d source: (%f,%f) %f --> %f\n",
+                        i, x, y, source->psfMag, source->psfMag + magCorr);
+            }
+#endif
+
             source->psfMag += magCorr;
         }
+
+#if 1
+// Size of fake image; set by hand because it's trouble to get it from other places
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sources, NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "corrected_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
     }
     psFree(trans);
@@ -162 +225 @@
         psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
                                                transThresh, starRej, starSys);
+        for (int i = 0; i < num; i++) {
+            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+        }
         psFree(trans);
     }