Changeset 40179 for branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c

Timestamp:

Oct 24, 2017, 4:20:46 PM (9 years ago)

Author:

watersc1

Message:

HSC changes that haven't been committed yet. Includes a number of configuration changes, improved glint masking, unfinished ghost masking, and the PSF model that I forgot to include previously.

File:

• branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c (modified) (2 diffs)

branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c

@@ -40 +40 @@
  * calculate ghost FPA and Chip positions for the stars loaded on the FPA
  */
+
+                    
+
 bool psastroLoadGhosts (pmConfig *config) {
 
@@ -256 +259 @@
     return false;
 }
+
+
+ 
+bool psastroLoadGhostsHSC (pmConfig *config) {
+  bool status;
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+    float zeropt, exptime, MAX_MAG;
+    psMetadata *ghostModel = NULL;
+    psVector *C_terms;
+    psVector *R_terms;
+    float glintPixelScale;
+    
+    psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions");
+
+    // select the current recipe
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
+    if (!recipe) {
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe");
+        return false;
+    }
+
+    bool REFSTAR_MASK_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST");
+    if (!REFSTAR_MASK_GHOST) return true;
+
+    // Load model from file
+    char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL");
+    if (!strcasecmp(ghostFile, "NONE")) return true;
+
+    if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) {
+        psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model");
+        return false;
+    }
+
+    // Allocate FPA and refstars
+    pmFPAview *view = pmFPAviewAlloc (0);
+    // select the input astrometry data (also carries the refstars)
+    pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
+    if (!astrom) {
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
+        goto escape;
+    }
+    pmFPA *fpa = astrom->fpa;
+
+    char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID");
+    psMetadataItem *item = psMetadataLookup(glintItem->data.md, "GLINT.FILTER.TERM");
+    if (!item) {
+      psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM data missing");
+      return false;
+    }
+    if (item->type != PS_DATA_METADATA_MULTI) {
+      psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM not multi");
+      return false;
+    }
+    psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false);
+    psMetadataItem *refItem = NULL;
+    while ((refItem = psListGetAndIncrement(iter))) {
+      if (refItem->type != PS_DATA_METADATA) {
+        psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM entry not metadata");
+        return false;
+      }
+      char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
+      if (!status) {
+        // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
+        continue;
+      }
+      if (strcmp(refFilter, filter)) continue;
+
+      C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS");
+      R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS");
+      MAX_MAG = psMetadataLookupF32(&status, refItem->data.md, "LIMIT");
+      glintPixelScale = psMetadataLookupF32(&status, refItem->data.md, "SCALE");
+    }
+
+    // really error-out here?  or just skip?
+    //    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa, recipe)) {
+    //        psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
+    //  goto escape;
+    //    }
+    
+    // recipe values are given in instrumental magnitudes
+    // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime)
+    //    float MagOffset = zeropt + 2.5*log10(exptime);
+    //    MAX_MAG += MagOffset;
+
+    // this loop selects the matched stars for all chips
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+        if (!chip->fromFPA) { continue; }
+
+        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+            // process each of the readouts
+            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+                if (! readout->data_exists) { continue; }
+
+                // select the raw objects for this readout (loaded in psastroChooseRefstars.c)
+                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
+                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only 
+                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
+                if (refstars == NULL) { continue; }
+
+                // identify the bright stars of interest
+                for (int i = 0; i < refstars->n; i++) {
+                    pmAstromObj *ref = refstars->data[i];
+                    psTrace("psastro.ghost",5,"Begin ghost %d/%ld: MAX_MAG: %g; ref_mag: %g @ (%.10g,%.10g) TPA: (%f %f)",
+                            i,refstars->n,MAX_MAG,ref->Mag,ref->sky->r * 180 / M_PI,ref->sky->d * 180.0 / M_PI,ref->TP->x,ref->RP->y);
+                    if (ref->Mag > MAX_MAG) continue;
+                    
+                    psastroGhost *ghost = psastroGhostAlloc ();
+
+                    double R_TPA  = sqrt(pow(star->TP->y,2) + pow(star->TP->x,2));
+                    double theta0 = atan2(star->TP->y,star->TP->x);
+                    double star_radius_deg = R_TPA * glintPixelScale;
+                    if (star_radius_deg > 0.5) { continue; }
+
+                    // Calculate expected position.
+                    double C = C_terms->data.F32[0] + C_terms->data.F32[1] * star_radius_deg + C_terms->data.F32[2] * pow(star_radius_deg,2) +
+                      C_terms->data.F32[3] * pow(star_radius_deg,3) + C_terms->data.F32[4] * pow(star_radius_deg,4) +
+                      C_terms->data.F32[5] * pow(star_radius_deg,5) + C_terms->data.F32[6] * pow(star_radius_deg,6) +
+                      C_terms->data.F32[7] * pow(star_radius_deg,7);
+                    double R = R_terms->data.F32[0] + R_terms->data.F32[1] * star_radius_deg + R_terms->data.F32[2] * pow(star_radius_deg,2) +
+                      R_terms->data.F32[3] * pow(star_radius_deg,3) + R_terms->data.F32[4] * pow(star_radius_deg,4) +
+                      R_terms->data.F32[5] * pow(star_radius_deg,5) + R_terms->data.F32[6] * pow(star_radius_deg,6) +
+                      R_terms->data.F32[7] * pow(star_radius_deg,7);
+                    R = 13.5 * (R / 0.015); // Convert to TPA units.
+
+                    psPlane *fp = psPlaneAlloc();
+                    psPlane *tp = psPlaneAlloc();
+
+                    tp->x = 13.5 * (C * cos(theta0) / 0.015 + 12.78);
+                    tp->y = 13.5 * (C * sin(theta0) / 0.015 + 57.74);
+                    psPlaneTransformApply(fp,fpa->fromTPA, tp);
+                    
+                    ghost->srcFP->x = ref->FP->x; 
+                    ghost->srcFP->y = ref->FP->y;
+                    ghost->FP->x    = fp->x;
+                    ghost->FP->y    = fp->y;
+
+                    ghost->inner.major = 0.0;
+                    ghost->inner.minor = 0.0;
+                    ghost->outer.major = R;
+                    ghost->outer.minor = R;
+                    
+                    pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y);
+                    
+
+                    
+                    if (ghostChip) {
+                      psTrace("psastro.ghost",5,"   in ghost: %d/%ld: ref: (%f,%f) or (%s) ghost: (%f,%f) or (%f,%f,%s) %f inner: (%f,%f,%f) outer: (%f,%f,%f)",
+                              i,refstars->n,
+                              ref->FP->x,ref->FP->y,
+                              psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"),
+                              ghost->FP->x,ghost->FP->y,
+                              ghost->chip->x,ghost->chip->y,psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME"),
+                              rSrc,
+                              ghost->inner.major,ghost->inner.minor,ghost->inner.theta,
+                              ghost->outer.major,ghost->outer.minor,ghost->outer.theta);
+/*                    psWarning("   GHOST_DATA: %d/%ld: ref: (%f,%f) or (%f,%f,%s) ghost: (%f,%f) or (%f,%f,%s) %f inner: (%f,%f,%f) outer: (%f,%f,%f)", */
+/*                            i,refstars->n, */
+/*                            ref->FP->x,ref->FP->y, */
+/*                            ref_chip_x,ref_chip_y,psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"), */
+/*                            ghost->FP->x,ghost->FP->y, */
+/*                            ghost->chip->x,ghost->chip->y,psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME"), */
+/*                            rSrc, */
+/*                            ghost->inner.major,ghost->inner.minor,ghost->inner.theta, */
+/*                            ghost->outer.major,ghost->outer.minor,ghost->outer.theta); */
+                    }
+                    else {
+                      psTrace("psastro.ghost",5,"   in ghost: %d/%ld: ref: (%f,%f) ghost: (%f,%f) or (%f,%f,%s) inner: (%f,%f,%f) outer: (%f,%f,%f)",
+                              i,refstars->n,
+                              ref->FP->x,ref->FP->y,
+                              ghost->FP->x,ghost->FP->y,
+                              ghost->chip->x,ghost->chip->y,"NONE",
+                              ghost->inner.major,ghost->inner.minor,ghost->inner.theta,
+                              ghost->outer.major,ghost->outer.minor,ghost->outer.theta);
+                    }                 
+                      
+
+                            
+                    if (!ghostChip) goto skip;
+                    if (!ghostChip->cells) goto skip;
+                    if (!ghostChip->cells->n) goto skip;
+
+                    
+                    pmCell *ghostCell = ghostChip->cells->data[0];
+                    if (!ghostCell) goto skip;
+                    if (!ghostCell->readouts) goto skip;
+                    if (!ghostCell->readouts->n) goto skip;
+                    pmReadout *ghostReadout = ghostCell->readouts->data[0];
+                    if (!ghostReadout) goto skip;
+
+                    psArray *ghosts = psMetadataLookupPtr (&status, ghostReadout->analysis, "PSASTRO.GHOSTS");
+                    if (ghosts == NULL) { 
+                        ghosts = psArrayAllocEmpty (100);
+                        if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) {
+                          psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout");
+                          goto escape;
+                        }
+                        psFree (ghosts);
+                    }
+
+                    psArrayAdd (ghosts, 100, ghost);
+
+                skip:
+                    
+                    psFree (ghost);
+                }
+            }
+        }
+    }
+
+    psastroExtractFreeChipBounds();
+
+    psFree (C_terms);
+    psFree (R_terms);
+    psFree (ghostModel);
+    psFree (view);
+    return true;
+
+escape:
+    psFree (C_terms);
+    psFree (R_terms);
+    psFree (ghostModel);
+    psFree (view);
+    return false;
+}