Changeset 40490

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/psastro.h (modified) (1 diff)
• src/psastroLoadGhosts.c (modified) (4 diffs)
• src/psastroLoadGlints.c (modified) (3 diffs)
• src/psastroMaskUpdates.c (modified) (9 diffs)

trunk/psastro/src/psastro.h

@@ -60 +60 @@
 psastroGhost     *psastroGhostAlloc (void);
 bool              psastroLoadGhosts (pmConfig *config);
+bool              psastroLoadGhostsHSC (pmConfig *config);
 
 bool              psastroDataSave (pmConfig *config, psMetadata *stats);

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;
+}

trunk/psastro/src/psastroLoadGlints.c

@@ -44 +44 @@
     double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO");
     double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH");
+    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
+    
 
     // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items)
@@ -95 +97 @@
         psProject (star->TP, star->sky, fpa->toSky);
         psPlaneTransformApply (star->FP, fpa->fromTPA, star->TP);
-        fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f\n", star->Mag, star->FP->x, star->FP->y);
+        fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f (%f %f) %8.1f %8.1f\n", star->Mag, star->FP->x, star->FP->y, star->sky->r * PS_DEG_RAD, star->sky->d * PS_DEG_RAD, star->TP->x,star->TP->y);
 
         // find the GLINT.REGION this star lands in (if any)
@@ -298 +300 @@
                 }
             }
+            if (!strcasecmp(glintType, "HSC")) {
+              // It's inefficient to keep looking these up.
+              double GLINT_RADIUS_INNER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.INNER");
+              double GLINT_RADIUS_OUTER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.OUTER");
+              
+              // double R_FPA = sqrt(pow(star->FP->y,2) + pow(star->FP->x,2));
+              double R_FPA =  sqrt(pow(star->TP->y,2) + pow(star->TP->x,2));
+              //              R_FPA *= 1.30;
+              if (R_FPA < GLINT_RADIUS_INNER) { continue; }
+              if (R_FPA > GLINT_RADIUS_OUTER) { continue; }
+              psTrace("psastro.masks",4,"HSC_GLINT_STARS: %f %f %f\n",
+                      //                      star->FP->x,star->FP->y,star->Mag);
+                      star->TP->x,star->TP->y,star->Mag);
+              psVector *C_terms;
+              psVector *R_terms;
+              double scale_factor = 1.0;
+              
+              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");
+                scale_factor = psMetadataLookupF32(&status, refItem->data.md, "FACTOR");
+                double limit = psMetadataLookupF32(&status, refItem->data.md, "LIMIT");
+                limit += MagOffset;
+                if (star->Mag > limit) { continue; }
+              }
+              // check 
+              double glintPixelScale = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.PIXEL.SCALE");
+              //              psMetadata *GLINT_PARAMETER_SET = psMetadataLookupPtr(&status, glintItem->data.md, "GLINT.FILTER.TERM");
+              //              psMetadata *GLINT_PARAMETERS    = psMetadataLookupPtr(&status, GLINT_PARAMETER_SET, filter);
+              
+
+              //              double theta0 = atan2(star->FP->y,star->FP->x);
+              double theta0 = atan2(star->TP->y,star->TP->x);
+              double star_radius_deg = R_FPA * glintPixelScale;
+
+              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);
+              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);
+
+              double inner_edge_angle, outer_edge_angle;
+              if (star_radius_deg < 0.8831) {
+                outer_edge_angle = 0.0;
+              }
+              else if (star_radius_deg < 0.9368) {
+                outer_edge_angle = 46.2985 * sqrt(star_radius_deg - 0.8831);
+              }
+              else {
+                outer_edge_angle = -2.67 + 14.3 * star_radius_deg;
+              }
+              if (star_radius_deg < 0.964) {
+                inner_edge_angle = -243.866 * pow(star_radius_deg - 0.932,2) + 2.4636;
+              }
+              else {
+                inner_edge_angle = 66.2061 * sqrt(star_radius_deg - 0.9635);
+              }
+
+              inner_edge_angle = inner_edge_angle * PS_RAD_DEG;
+              outer_edge_angle = outer_edge_angle * PS_RAD_DEG;
+              
+              // These define the ends of a single arc in arc-centric coordinates
+              double x1 = C - R * cos(inner_edge_angle);
+              double y1 = R * sin(inner_edge_angle);
+              double x2 = C - R * cos(outer_edge_angle);
+              double y2 = R * sin(outer_edge_angle);
+
+              // Create the endpoints for the pair of arcs in device coordinates (in millimeters)
+              double theta1 = theta0 + M_PI / 2.0;
+              double x1s = (x1 * cos(theta1) - y1 * sin(theta1)) * scale_factor;
+              double y1s = (y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor;
+              double x1e = (x2 * cos(theta1) - y2 * sin(theta1)) * scale_factor;
+              double y1e = (y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor;
+              
+              double x2s = (x1 * cos(theta1) + y1 * sin(theta1)) * scale_factor;
+              double y2s = (-1.0 * y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor;
+              double x2e = (x2 * cos(theta1) + y2 * sin(theta1)) * scale_factor;
+              double y2e = (-1.0 * y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor;
+
+              psVector *x_start = psVectorAlloc(4,PS_TYPE_F32);
+              psVector *y_start = psVectorAlloc(4,PS_TYPE_F32);
+              psVector *x_end   = psVectorAlloc(4,PS_TYPE_F32);
+              psVector *y_end   = psVectorAlloc(4,PS_TYPE_F32);
+                      
+              // CZW: I know this looks like a typo, but trust me, it's not a typo.
+              psPlane *fp = psPlaneAlloc();
+              psPlane *tp = psPlaneAlloc();
+
+              fprintf(stderr,"HSC GLINT %f: x1x2(%f %f) (%f %f) x1sx1e (%f %f %f %f) x2sx2e (%f %f %f %f)\n",
+                      star->Mag,
+                      x1,y1,x2,y2,
+                      x1s,y1s,x1e,y1e,
+                      x2s,y2s,x2e,y2e);
+                      
+              
+              tp->x = 13.5 * (y1s / 0.015 + 12.78);
+              tp->y = 13.5 * (x1s / -0.015 + 57.74);
+              psPlaneTransformApply(fp,fpa->fromTPA, tp);
+              x_start->data.F32[0] = fp->x;
+              y_start->data.F32[0] = fp->y;
+              //              x_start->data.F32[0] = tp->x;
+              //              y_start->data.F32[0] = tp->y;
+
+              tp->x = 13.5 * (y1e / 0.015 + 12.78);
+              tp->y = 13.5 * (x1e / -0.015 + 57.74);
+              psPlaneTransformApply(fp,fpa->fromTPA, tp);
+              x_end->data.F32[0]   = fp->x;
+              y_end->data.F32[0]   = fp->y;
+              //              x_end->data.F32[0]   = tp->x;
+              //              y_end->data.F32[0]   = tp->y;
+
+              x_start->data.F32[1] = x_end->data.F32[0];  
+              y_start->data.F32[1] = y_end->data.F32[0];  
+              x_end->data.F32[1]   = x_start->data.F32[0];
+              y_end->data.F32[1]   = y_start->data.F32[0];
+
+              tp->x = 13.5 * (y2s / 0.015 + 12.78);
+              tp->y = 13.5 * (x2s / -0.015 + 57.74);
+              psPlaneTransformApply(fp,fpa->fromTPA, tp);
+              x_start->data.F32[2] = fp->x;
+              y_start->data.F32[2] = fp->y;
+              //              x_start->data.F32[2] = tp->x;
+              //              y_start->data.F32[2] = tp->y;
+
+              tp->x = 13.5 * (y2e / 0.015 + 12.78);
+              tp->y = 13.5 * (x2e / -0.015 + 57.74);
+              psPlaneTransformApply(fp,fpa->fromTPA, tp);
+              x_end->data.F32[2]   = fp->x;
+              y_end->data.F32[2]   = fp->y;
+              //              x_end->data.F32[2]   = tp->x;
+              //              y_end->data.F32[2]   = tp->y;
+              
+              x_start->data.F32[3] = x_end->data.F32[2];  
+              y_start->data.F32[3] = y_end->data.F32[2];  
+              x_end->data.F32[3]   = x_start->data.F32[2];
+              y_end->data.F32[3]   = y_start->data.F32[2];
+
+              psFree(fp);
+              psFree(tp);
+              
+              for (int glint_point = 0; glint_point < 4; glint_point++) {
+                for (int nChip = 0; nChip < fpa->chips->n; nChip++) {
+                  pmChip *chip = fpa->chips->data[nChip];
+                  if (!chip) continue;
+                  
+                  if (!psastroFindChipInYrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) {
+                    continue;
+                  }
+                  if (!psastroFindChipInXrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) {
+                    continue;
+                  }
+
+                  double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length;
+                  psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point]);
+                  psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, x_end->data.F32[glint_point], y_end->data.F32[glint_point]);
+
+                  chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0);
+                  glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2));
+
+                  // select the 0th readout of the 0th cell for this chip
+                  if (!chip->cells) continue;
+                  if (!chip->cells->n) continue;
+                  pmCell *glintCell = chip->cells->data[0];
+                  if (!glintCell) continue;
+                  if (!glintCell->readouts) continue;
+                  if (!glintCell->readouts->n) continue;
+                  pmReadout *glintReadout = glintCell->readouts->data[0];
+                  if (!glintReadout) continue;
+                  
+                  // save the glints on the readout->analysis metadata, creating if needed
+                  psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS.HSC");
+                  if (glints == NULL) { 
+                    glints = psArrayAllocEmpty (100);
+                    if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS.HSC", PS_DATA_ARRAY, "astrometry matches", glints)) {
+                      psWarning("failure to add glints to readout");
+                      psFree (glints);
+                      continue;
+                  }
+                    psFree (glints);
+                  }
+                  
+                  fprintf (stderr, "glint %s : %d %f,%f to %f,%f (%f %f %f)\n", glintType, nChip, xChip0, yChip0, xChip1, yChip1, glint_length, glintWidth, chip_angle);
+                  psTrace("psastro.masks",4,"HSC_GLINT: Star: %f %f Glint CR %f %f Parameters: %f %f %f Ends: %f %f -> %f %f Chip: %f %f -> %f %f @ %s %f\n",
+                          star->FP->x,star->FP->y,
+                          C,R,
+                          inner_edge_angle,outer_edge_angle,theta0,
+                          x_start->data.F32[glint_point],y_start->data.F32[glint_point],x_end->data.F32[glint_point],y_end->data.F32[glint_point],
+                          xChip0,yChip0,xChip1,yChip1,
+                          psMetadataLookupStr(&status,glintReadout->parent->parent->concepts,"CHIP.NAME"),chip_angle);
+                  psVector *glint = psVectorAlloc(5,PS_TYPE_F32);
+                  glint->data.F32[0] = xChip0;
+                  glint->data.F32[1] = yChip0;
+                  glint->data.F32[2] = glint_length * pixelScale / pixelScale;
+                  glint->data.F32[3] = glintWidth;
+                  glint->data.F32[4] = chip_angle;
+
+                  psArrayAdd (glints, 100, glint);
+
+                  psFree (glint);
+                } // End loop over chips
+              } // End loop over glint endpoints
+
+              psFree(x_start);
+              psFree(x_end);
+              psFree(y_start);
+              psFree(y_end);
+            } // End HSC glint block
+            
         }
     }

trunk/psastro/src/psastroMaskUpdates.c

@@ -129 +129 @@
     
     bool REFSTAR_MASK_BLEED                = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED");
-
+    bool REFSTAR_MASK_BLEED_ORIENTATION_X  = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_X");
+    bool REFSTAR_MASK_BLEED_ORIENTATION_Y  = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_Y");
+
+    
     double REFSTAR_MASK_MAX_MAG            = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG");
     double REFSTAR_MASK_SATSTAR_MAG_SLOPE  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_SLOPE");
@@ -143 +146 @@
     double REFSTAR_MASK_BLEED_MAG_SLOPE    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE");
 
+    double REFSTAR_MASK_BLEED_MAG_MAX_X    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_X");
+    double REFSTAR_MASK_BLEED_MAG_SLOPE_X  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_X");
+    double REFSTAR_MASK_BLEED_MAG_MAX_Y    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_Y");
+    double REFSTAR_MASK_BLEED_MAG_SLOPE_Y  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_Y");
+
     double REFSTAR_MASK_CROSSTALK_MAG_MAX  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX");
     double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE");
@@ -180 +188 @@
     }
 
+    psTrace("psastro.masks",2,"Configuration loaded for masking.");
+    
     // 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)
@@ -187 +197 @@
     REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset;
     REFSTAR_MASK_BLEED_MAG_MAX += MagOffset;
+    REFSTAR_MASK_BLEED_MAG_MAX_X += MagOffset;
+    REFSTAR_MASK_BLEED_MAG_MAX_Y += MagOffset;
+    
     GHOST_MAX_MAG += MagOffset;
 
+    psTrace("psastro.masks",2,"Magnitudes: max_mag: %f satstar: %f satspike: %f bleed: %f ghost %f\n",
+            REFSTAR_MASK_MAX_MAG,REFSTAR_MASK_SATSTAR_MAG_MAX,REFSTAR_MASK_SATSPIKE_MAG_MAX,REFSTAR_MASK_BLEED_MAG_MAX,GHOST_MAX_MAG);
+    
     // select the output mask image :: we mosaic to chip mosaic format
     pmFPAfile *outMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT.MASK");
@@ -287 +303 @@
                 if (!refstars) continue;
 
+                psTrace("psastro.masks",2,"In Readout loop.");
                 // we need to generate the following masks regions:
                 // 1) circle around the saturated stars (scaled by magnitude)
@@ -301 +318 @@
                         }
                     }
-
+                    psTrace("psastro.masks",4,"In refstar loop: %d/%ld %f %f\n",
+                            i,refstars->n,ref->Mag,REFSTAR_MASK_MAX_MAG);
                     if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;
 
@@ -334 +352 @@
 
                         // The length should also be a function of the image background level
-
+                        psTrace("psastro.masks",4,"Masking: Radius: %f Theta: %f Length: %f Width: %f\n",
+                                radius,Theta,spikeLength,spikeWidth);
+                        psTrace("psastro.masks",4,"Also: %f %f %f %f\n",ref->chip->x,ref->chip->y,ROTANGLE,REFSTAR_MASK_SATSTAR_POS_ZERO);
                         psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
                     }
@@ -349 +369 @@
                             float xCell = 0.0;
                             float yCell = 0.0;
-                            float xEnd = 0.0;
-                            float yEnd = 0.0;
+                            float xStart = 0.0;
+                            float xEnd   = 0.0;
+                            float yStart = 0.0;
+                            float yEnd   = 0.0;
+                            //                            float xEnd = 0.0;
+                            //                            float yEnd = 0.0;
+                            float width = 0.0;
+                            float length = 0.0;
                             // find coordinate of star on cell
                             pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
                             // find coordinate of end-point on chip
 
-                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
-                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
-
-                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
-                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
+                            if (REFSTAR_MASK_BLEED_ORIENTATION_X) {
+                              length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_X*(REFSTAR_MASK_BLEED_MAG_MAX_X - ref->Mag));
+                              width  = REFSTAR_MASK_BLEED_MAG_SLOPE  *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
+
+                              pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - length, yCell - 0.5 * width);
+                              pmChipCoordsForCell (&xEnd, &yEnd,     refCell, xCell + length, yCell + 0.5 * width);
+                              psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd, (int) yEnd + 1);
+                            }
+                            if (REFSTAR_MASK_BLEED_ORIENTATION_Y) {
+                              length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_Y*(REFSTAR_MASK_BLEED_MAG_MAX_Y - ref->Mag));
+                              width  = REFSTAR_MASK_BLEED_MAG_SLOPE  *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
+
+                              pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - 0.5 * width, yCell - length);
+                              pmChipCoordsForCell (&xEnd, &yEnd,     refCell, xCell + 0.5 * width, yCell + length);
+                              psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd + 1, (int) yEnd);
+                            }
+                            
                         }
                     }
@@ -392 +430 @@
                     }
                 }
-
+                // Because it's very GPC1-specific, handle alternate glint types here
+                psArray *hsc_glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS.HSC");
+                if (hsc_glints) {
+                  for (int i = 0; i < hsc_glints->n; i++) {
+                    // Contains (start_x, start_y, length, width, theta)
+                    psVector *glint = hsc_glints->data[i];
+
+                    psTrace("psastro.masks", 4, "Masking glint %d on chip %s with params (%f,%f,%f,%f,%f)\n",
+                            i, psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"),
+                            glint->data.F32[0], glint->data.F32[1],
+                            glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]);
+
+                    psastroMaskBox (readoutMask->mask, glintMaskValue, glint->data.F32[0], glint->data.F32[1],
+                                    glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]);
+                  }
+                }
+                
                 // Select the crosstalk artifact regions for this readout, and mask a circular region
                 // corresponding to the source star