Changeset 17935 for branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c

Timestamp:

Jun 5, 2008, 3:29:09 AM (18 years ago)

Author:

eugene

Message:

added forced photometry code; compiles, not tested

File:

• branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c (modified) (1 diff)

branches/eam_branch_20080511/ppSim/src/ppSimLoadForceSources.c

@@ -1 +1 @@
 # include "ppSim.h"
+
+// Reset a pointer: free and set to NULL
+#define RESET(PTR) \
+    psFree(PTR); \
+    PTR = NULL;
 
 psArray *ppSimLoadForceSources(pmConfig *config, const pmFPAview *view) {
 
-    psArray *result = psArrayAlloc(0);
-    return result;
+    bool status;
 
-# if (0)
+    pmFPAfile *outFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.OUTPUT");
+    psAssert(outFile, "missing PPSIM.OUTPUT");
+    psAssert(outFile->fpa, "missing fpr for PPSIM.OUTPUT");
 
-    // de-activate all files except PSASTRO.REFSTARS
-    pmFPAfileActivate (config->files, false, NULL);
-    pmFPAfileActivate (config->files, true, "PSASTRO.OUT.REFSTARS");
+    pmChip *outChip = pmFPAviewThisChip (view, outFile->fpa);
+    pmPSF *psf = psMetadataLookupPtr (&status, outChip->analysis, "PSPHOT.PSF");
+    assert (psf);
 
-    readout = pmFPAviewThisReadout (view, fpa);
-    if (! readout->data_exists) { continue; }
+    psArray *spots = psMetadataLookupPtr(&status, outFile->fpa->analysis, "FORCED.SPOTS");
+    PS_ASSERT_PTR_NON_NULL (spots, NULL);
 
-                // read WCS data from the corresponding header
-                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
+    // in this program, we are looping over the output image, rather than the input as in ppImage
+    pmFPAfile *srcFile = psMetadataLookupPtr(NULL, config->files, "PPSIM.REAL.SOURCES");
+    psAssert(srcFile, "missing PPSIM.REAL.SOURCES");
 
-                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
-                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");
+    pmFPA     *fpa     = srcFile->fpa;
+    pmChip    *chip    = pmFPAviewThisChip (view, fpa);
+    // pmCell    *cell    = pmFPAviewThisCell (view, fpa);
+    pmReadout *readout = pmFPAviewThisReadout (view, fpa);
 
-                float minX = -fieldPadding*Nx;
-                float maxX = (1+fieldPadding)*Nx;
-                float minY = -fieldPadding*Ny;
-                float maxY = (1+fieldPadding)*Ny;
+    // XXX we should have only one cell and readout per chip, right?
 
-                // the refstars is a subset within range of this chip
-                psArray *refstars = psArrayAllocEmpty (100);
+    // read WCS from existing output from psphot & psastro
+    // XXX can we be more flexible?  ie, use the header of PPSIM.CHIP?
+    bool bilevelAstrometry = false;
+    status = psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
+    psAssert (status, "failed to set astrometry");
 
-                // select the reference objects within range of this readout
-                // project the reference objects to this chip
-                for (int i = 0; i < refs->n; i++) {
-                    pmAstromObj *ref = pmAstromObjCopy(refs->data[index->data.S32[i]]);
+    // the pixel scale is only used to set the focal-plane coordinate system
+    // setting this value to 1 makes this system have the units of pixels
+    float pixelScale = 1.0;
+    status = psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry);
+    psAssert (status, "failed to set astrometry");
 
-                    psProject (ref->TP, ref->sky, fpa->toSky);
-                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
-                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+    float minX = 0.0;
+    float maxX = readout->image->numCols;
+    float minY = 0.0;
+    float maxY = readout->image->numRows;
 
-                    // limit the X,Y range of the refs to the selected chip
-                    if (ref->chip->x < minX) goto skip;
-                    if (ref->chip->x > maxX) goto skip;
-                    if (ref->chip->y < minY) goto skip;
-                    if (ref->chip->y > maxY) goto skip;
+    // the sources is a subset within range of this chip
+    psArray *sources = psArrayAllocEmpty (100);
 
-                    psArrayAdd (refstars, 100, ref);
-                skip:
-                    psFree (ref);
+    // Select the spots within range of this readout.  Project the spots to this chip
+    for (int i = 0; i < spots->n; i++) {
+        pmAstromObj *obj = spots->data[i];
 
-                    if (nMax && (refstars->n >= nMax)) break;
-                }
-                psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
+        psProject (obj->TP, obj->sky, fpa->toSky);
+        psPlaneTransformApply (obj->FP, fpa->fromTPA, obj->TP);
+        psPlaneTransformApply (obj->chip, chip->fromFPA, obj->FP);
 
-                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
-                psFree (refstars);
+        // limit the X,Y range of the objs to the selected chip
+        if (obj->chip->x < minX) continue;
+        if (obj->chip->x > maxX) continue;
+        if (obj->chip->y < minY) continue;
+        if (obj->chip->y > maxY) continue;
 
-                if (matchLumFunc) {
-                    // in this case, no PSASTRO.REFSTARS is added to readout->analysis
-                    if (!psastroRefstarSubset (readout)) {
-                        psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n");
-                        psFree (index);
-                        psFree (view);
-                        return false;
-                    }
-                }
-            }
-        }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
+        // convert the pmAstromObj to pmSource
+
+        // instantiate a model for the PSF at this location, Io = 1.0
+        pmModel *model = pmModelFromPSFforXY (psf, obj->chip->x, obj->chip->y, 1.0);
+
+        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
+        // XXX use a fixed radius??
+        // float radius = model->modelRadius (model->params, roughNoise);
+        // radius = PS_MAX (radius, 1.0);
+        float radius = 5.0;
+
+        // construct a source, with model flux pixels set, based on the model
+        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+
+        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
+        // XXX i should be applying the gain and the correct effective area
+        // psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+        // psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+
+        // these are not really needed since we will be fitting for them
+        source->psfMag = NAN;
+        source->errMag = NAN;
+
+        // insert the source flux in the image
+        // XXX not sure the offset is really 0,0
+        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, 0.0, 0.0);
+        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, 0.0, 0.0);
+
+        // Blow away the image parts of the source, which makes the memory explode
+        RESET(source->pixels);
+        RESET(source->weight);
+        RESET(source->maskObj);
+        RESET(source->maskView);
+        RESET(source->modelFlux);
+        RESET(source->psfFlux);
+        RESET(source->blends);
+
+        psArrayAdd (sources, 100, source);
+        psFree(source);                 // Drop local reference
     }
-    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
+    psTrace ("psastro", 4, "Added %ld sources\n", sources->n);
 
-    // activate all files except PSASTRO.OUTPUT
-
-    pmFPAfileActivate (config->files, true, NULL);
-    pmFPAfileActivate (config->files, false, "PSASTRO.OUT.REFSTARS");
-# endif
-
-    return NULL;
+    return sources;
 }