Index: trunk/psastro/src/psastroLoadGlints.c
===================================================================
--- trunk/psastro/src/psastroLoadGlints.c	(revision 24645)
+++ trunk/psastro/src/psastroLoadGlints.c	(revision 40490)
@@ -44,4 +44,6 @@
     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,5 +97,5 @@
 	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,4 +300,232 @@
 		}
 	    }
+	    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
+	    
 	}
     }