Changeset 40490 for trunk/psastro/src/psastroLoadGhosts.c

Timestamp:

Jun 27, 2018, 4:11:05 PM (8 years ago)

Author:

eugene

Message:

add options and elements for HSC ghosts and glints

Location:

trunk/psastro

Files:

• . (modified) (1 prop)
• src/psastroLoadGhosts.c (modified) (4 diffs)

trunk/psastro/src/psastroLoadGhosts.c

@@ -40 +40 @@
  * calculate ghost FPA and Chip positions for the stars loaded on the FPA
  */
+
+                    
+
 bool psastroLoadGhosts (pmConfig *config) {
 
@@ -58 +61 @@
     psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions");
 
+    // XXX this needs to be camera-dependent
+    if (0) {
+      psWarning("Doing the hard coded switch to HSC ghosts still.");
+      return (psastroLoadGhostsHSC(config));
+    } 
+
+    psWarning("hard coded switch to GPC ghosts.");
+
     // select the current recipe
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
@@ -68 +79 @@
     if (!REFSTAR_MASK_GHOST) return true;
 
+    
     char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL");
     if (!strcasecmp(ghostFile, "NONE")) return true;
@@ -256 +268 @@
     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;
+  psLogMsg("psastro", PS_LOG_INFO, "ghost file %s",ghostFile);
+  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(ghostModel, "GHOST.MODEL.HSC");
+  if (!item) {
+    psLogMsg ("psastro", PS_LOG_INFO, "GHOST.MODEL.HSC data missing");
+    return false;
+  }
+  if (item->type != PS_DATA_METADATA_MULTI) {
+    psLogMsg ("psastro", PS_LOG_INFO, "GHOST.MODEL.HSC 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, "GHOST.MODEL.HSC 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");
+    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];
+          if (ref->Mag > MAX_MAG) continue;
+          
+          psastroGhost *ghost = psastroGhostAlloc ();
+          
+          double R_TPA  = sqrt(pow(ref->TP->y,2) + pow(ref->TP->x,2));
+          double theta0 = atan2(ref->TP->y,ref->TP->x);
+          double star_radius_deg = R_TPA * glintPixelScale;
+          
+          //      psTrace("psastro.ghost",5,
+          psLogMsg("psastro",PS_LOG_INFO,
+                  "Begin ghost %d/%ld: MAX_MAG: %g; ref_mag: %g @ (%.10g,%.10g) TPA: (%f %f) R: %f %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->TP->y,
+                  R_TPA,star_radius_deg,theta0);
+          
+          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 = R / 0.015;
+          
+          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) {
+            psLogMsg("psastro",PS_LOG_INFO,
+                     //     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"),
+                    C,
+                    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;
+}