- Timestamp:
- Oct 24, 2017, 4:20:46 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/czw_branch/20170908/psastro/src/psastroLoadGlints.c
r24645 r40179 44 44 double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO"); 45 45 double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH"); 46 double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); 47 46 48 47 49 // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items) … … 95 97 psProject (star->TP, star->sky, fpa->toSky); 96 98 psPlaneTransformApply (star->FP, fpa->fromTPA, star->TP); 97 fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f \n", star->Mag, star->FP->x, star->FP->y);99 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); 98 100 99 101 // find the GLINT.REGION this star lands in (if any) … … 298 300 } 299 301 } 302 if (!strcasecmp(glintType, "HSC")) { 303 // It's inefficient to keep looking these up. 304 double GLINT_RADIUS_INNER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.INNER"); 305 double GLINT_RADIUS_OUTER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.OUTER"); 306 307 // double R_FPA = sqrt(pow(star->FP->y,2) + pow(star->FP->x,2)); 308 double R_FPA = sqrt(pow(star->TP->y,2) + pow(star->TP->x,2)); 309 // R_FPA *= 1.30; 310 if (R_FPA < GLINT_RADIUS_INNER) { continue; } 311 if (R_FPA > GLINT_RADIUS_OUTER) { continue; } 312 psTrace("psastro.masks",4,"HSC_GLINT_STARS: %f %f %f\n", 313 // star->FP->x,star->FP->y,star->Mag); 314 star->TP->x,star->TP->y,star->Mag); 315 psVector *C_terms; 316 psVector *R_terms; 317 double scale_factor = 1.0; 318 319 char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); 320 psMetadataItem *item = psMetadataLookup(glintItem->data.md, "GLINT.FILTER.TERM"); 321 if (!item) { 322 psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM data missing"); 323 return false; 324 } 325 if (item->type != PS_DATA_METADATA_MULTI) { 326 psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM not multi"); 327 return false; 328 } 329 psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false); 330 psMetadataItem *refItem = NULL; 331 while ((refItem = psListGetAndIncrement(iter))) { 332 if (refItem->type != PS_DATA_METADATA) { 333 psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM entry not metadata"); 334 return false; 335 } 336 char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); 337 if (!status) { 338 // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); 339 continue; 340 } 341 if (strcmp(refFilter, filter)) continue; 342 343 C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS"); 344 R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS"); 345 scale_factor = psMetadataLookupF32(&status, refItem->data.md, "FACTOR"); 346 double limit = psMetadataLookupF32(&status, refItem->data.md, "LIMIT"); 347 if (star->Mag > limit) { continue; } 348 } 349 // check 350 double glintPixelScale = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.PIXEL.SCALE"); 351 // psMetadata *GLINT_PARAMETER_SET = psMetadataLookupPtr(&status, glintItem->data.md, "GLINT.FILTER.TERM"); 352 // psMetadata *GLINT_PARAMETERS = psMetadataLookupPtr(&status, GLINT_PARAMETER_SET, filter); 353 354 355 // double theta0 = atan2(star->FP->y,star->FP->x); 356 double theta0 = atan2(star->TP->y,star->TP->x); 357 double star_radius_deg = R_FPA * glintPixelScale; 358 359 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); 360 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); 361 362 double inner_edge_angle, outer_edge_angle; 363 if (star_radius_deg < 0.8831) { 364 outer_edge_angle = 0.0; 365 } 366 else if (star_radius_deg < 0.9368) { 367 outer_edge_angle = 46.2985 * sqrt(star_radius_deg - 0.8831); 368 } 369 else { 370 outer_edge_angle = -2.67 + 14.3 * star_radius_deg; 371 } 372 if (star_radius_deg < 0.964) { 373 inner_edge_angle = -243.866 * pow(star_radius_deg - 0.932,2) + 2.4636; 374 } 375 else { 376 inner_edge_angle = 66.2061 * sqrt(star_radius_deg - 0.9635); 377 } 378 379 inner_edge_angle = inner_edge_angle * PS_RAD_DEG; 380 outer_edge_angle = outer_edge_angle * PS_RAD_DEG; 381 382 // These define the ends of a single arc in arc-centric coordinates 383 double x1 = C - R * cos(inner_edge_angle); 384 double y1 = R * sin(inner_edge_angle); 385 double x2 = C - R * cos(outer_edge_angle); 386 double y2 = R * sin(outer_edge_angle); 387 388 // Create the endpoints for the pair of arcs in device coordinates (in millimeters) 389 double theta1 = theta0 + M_PI / 2.0; 390 double x1s = (x1 * cos(theta1) - y1 * sin(theta1)) * scale_factor; 391 double y1s = (y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; 392 double x1e = (x2 * cos(theta1) - y2 * sin(theta1)) * scale_factor; 393 double y1e = (y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; 394 395 double x2s = (x1 * cos(theta1) + y1 * sin(theta1)) * scale_factor; 396 double y2s = (-1.0 * y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; 397 double x2e = (x2 * cos(theta1) + y2 * sin(theta1)) * scale_factor; 398 double y2e = (-1.0 * y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; 399 400 psVector *x_start = psVectorAlloc(4,PS_TYPE_F32); 401 psVector *y_start = psVectorAlloc(4,PS_TYPE_F32); 402 psVector *x_end = psVectorAlloc(4,PS_TYPE_F32); 403 psVector *y_end = psVectorAlloc(4,PS_TYPE_F32); 404 405 // CZW: I know this looks like a typo, but trust me, it's not a typo. 406 psPlane *fp = psPlaneAlloc(); 407 psPlane *tp = psPlaneAlloc(); 408 409 fprintf(stderr,"HSC GLINT: x1x2(%f %f) (%f %f) x1sx1e (%f %f %f %f) x2sx2e (%f %f %f %f)\n", 410 x1,y1,x2,y2, 411 x1s,y1s,x1e,y1e, 412 x2s,y2s,x2e,y2e); 413 414 415 tp->x = 13.5 * (y1s / 0.015 + 12.78); 416 tp->y = 13.5 * (x1s / -0.015 + 57.74); 417 psPlaneTransformApply(fp,fpa->fromTPA, tp); 418 x_start->data.F32[0] = fp->x; 419 y_start->data.F32[0] = fp->y; 420 x_start->data.F32[0] = tp->x; 421 y_start->data.F32[0] = tp->y; 422 423 tp->x = 13.5 * (y1e / 0.015 + 12.78); 424 tp->y = 13.5 * (x1e / -0.015 + 57.74); 425 psPlaneTransformApply(fp,fpa->fromTPA, tp); 426 x_end->data.F32[0] = fp->x; 427 y_end->data.F32[0] = fp->y; 428 x_end->data.F32[0] = tp->x; 429 y_end->data.F32[0] = tp->y; 430 431 x_start->data.F32[1] = x_end->data.F32[0]; 432 y_start->data.F32[1] = y_end->data.F32[0]; 433 x_end->data.F32[1] = x_start->data.F32[0]; 434 y_end->data.F32[1] = y_start->data.F32[0]; 435 436 tp->x = 13.5 * (y2s / 0.015 + 12.78); 437 tp->y = 13.5 * (x2s / -0.015 + 57.74); 438 psPlaneTransformApply(fp,fpa->fromTPA, tp); 439 x_start->data.F32[2] = fp->x; 440 y_start->data.F32[2] = fp->y; 441 x_start->data.F32[2] = tp->x; 442 y_start->data.F32[2] = tp->y; 443 444 tp->x = 13.5 * (y2e / 0.015 + 12.78); 445 tp->y = 13.5 * (x2e / -0.015 + 57.74); 446 psPlaneTransformApply(fp,fpa->fromTPA, tp); 447 x_end->data.F32[2] = fp->x; 448 y_end->data.F32[2] = fp->y; 449 x_end->data.F32[2] = tp->x; 450 y_end->data.F32[2] = tp->y; 451 452 x_start->data.F32[3] = x_end->data.F32[2]; 453 y_start->data.F32[3] = y_end->data.F32[2]; 454 x_end->data.F32[3] = x_start->data.F32[2]; 455 y_end->data.F32[3] = y_start->data.F32[2]; 456 457 psFree(fp); 458 psFree(tp); 459 460 for (int glint_point = 0; glint_point < 4; glint_point++) { 461 for (int nChip = 0; nChip < fpa->chips->n; nChip++) { 462 pmChip *chip = fpa->chips->data[nChip]; 463 if (!chip) continue; 464 465 if (!psastroFindChipInYrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { 466 continue; 467 } 468 if (!psastroFindChipInXrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { 469 continue; 470 } 471 472 double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length; 473 psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point]); 474 psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, x_end->data.F32[glint_point], y_end->data.F32[glint_point]); 475 476 chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0); 477 glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2)); 478 479 // select the 0th readout of the 0th cell for this chip 480 if (!chip->cells) continue; 481 if (!chip->cells->n) continue; 482 pmCell *glintCell = chip->cells->data[0]; 483 if (!glintCell) continue; 484 if (!glintCell->readouts) continue; 485 if (!glintCell->readouts->n) continue; 486 pmReadout *glintReadout = glintCell->readouts->data[0]; 487 if (!glintReadout) continue; 488 489 // save the glints on the readout->analysis metadata, creating if needed 490 psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS.HSC"); 491 if (glints == NULL) { 492 glints = psArrayAllocEmpty (100); 493 if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS.HSC", PS_DATA_ARRAY, "astrometry matches", glints)) { 494 psWarning("failure to add glints to readout"); 495 psFree (glints); 496 continue; 497 } 498 psFree (glints); 499 } 500 501 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); 502 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", 503 star->FP->x,star->FP->y, 504 C,R, 505 inner_edge_angle,outer_edge_angle,theta0, 506 x_start->data.F32[glint_point],y_start->data.F32[glint_point],x_end->data.F32[glint_point],y_end->data.F32[glint_point], 507 xChip0,yChip0,xChip1,yChip1, 508 psMetadataLookupStr(&status,glintReadout->parent->parent->concepts,"CHIP.NAME"),chip_angle); 509 psVector *glint = psVectorAlloc(5,PS_TYPE_F32); 510 glint->data.F32[0] = xChip0; 511 glint->data.F32[1] = yChip0; 512 glint->data.F32[2] = glint_length * pixelScale / pixelScale; 513 glint->data.F32[3] = glintWidth; 514 glint->data.F32[4] = chip_angle; 515 516 psArrayAdd (glints, 100, glint); 517 518 psFree (glint); 519 } // End loop over chips 520 } // End loop over glint endpoints 521 522 psFree(x_start); 523 psFree(x_end); 524 psFree(y_start); 525 psFree(y_end); 526 } // End HSC glint block 527 300 528 } 301 529 }
Note:
See TracChangeset
for help on using the changeset viewer.
