Changeset 40490 for trunk/psastro/src/psastroLoadGlints.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/psastroLoadGlints.c (modified) (3 diffs)

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