Index: trunk/psastro/src/psastroAstromGuess.c =================================================================== --- trunk/psastro/src/psastroAstromGuess.c (revision 21409) +++ trunk/psastro/src/psastroAstromGuess.c (revision 21422) @@ -6,6 +6,6 @@ * * @author IfA - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-07 02:03:34 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-09 21:25:34 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -42,6 +42,6 @@ psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE); if (!recipe) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!"); + return false; } @@ -49,5 +49,5 @@ bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); if (!status) { - useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); } @@ -55,6 +55,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -62,7 +62,7 @@ double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); - return false; - } + psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); + return false; + } psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32); @@ -80,5 +80,5 @@ bool bilevelAstrometry = false; if (!useModel) { - psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); + psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); } @@ -88,35 +88,35 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - if (!useModel) { - if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; - } + if (!useModel) { + if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; + } if (newFPA) { newFPA = false; - while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; - while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; + while (fpa->toSky->R < 0) fpa->toSky->R += 2.0*M_PI; + while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI; RAminSky = fpa->toSky->R - M_PI; RAmaxSky = fpa->toSky->R + M_PI; } - // report and save the current best guess for the chip 0,0 pixel coordinates - { - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - psVectorAppend (cornerL, ptFP.x); - psVectorAppend (cornerM, ptFP.y); - psVectorAppend (cornerP, ptTP.x); - psVectorAppend (cornerQ, ptTP.y); - psVectorAppend (cornerR, ptSky.r); - psVectorAppend (cornerD, ptSky.d); - } + // report and save the current best guess for the chip 0,0 pixel coordinates + { + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + psVectorAppend (cornerL, ptFP.x); + psVectorAppend (cornerM, ptFP.y); + psVectorAppend (cornerP, ptTP.x); + psVectorAppend (cornerQ, ptTP.y); + psVectorAppend (cornerR, ptSky.r); + psVectorAppend (cornerD, ptSky.d); + } // apply the new WCS guess data to all of the data in the readouts @@ -132,5 +132,5 @@ if (rawstars == NULL) { continue; } - *nStars += rawstars->n; + *nStars += rawstars->n; for (int i = 0; i < rawstars->n; i++) { pmAstromObj *raw = rawstars->data[i]; @@ -151,14 +151,14 @@ } - // dump or plot the resulting projected positions - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRawstars (rawstars, fpa, chip); - } - - psastroVisualPlotRawStars(rawstars, fpa, chip, recipe); - - if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotRawstars (rawstars, fpa, chip, recipe); - } + // dump or plot the resulting projected positions + if (psTraceGetLevel("psastro.dump") > 0) { + psastroDumpRawstars (rawstars, fpa, chip); + } + + pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe); + + if (psTraceGetLevel("psastro.plot") > 0) { + psastroPlotRawstars (rawstars, fpa, chip, recipe); + } } } @@ -170,12 +170,12 @@ psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR", PS_META_REPLACE, "", *nStars); if (*nStars == 0) { - psLogMsg ("psastro", 2, "no sources available for astrometry\n"); - psFree (view); - return true; - } - - psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", - DEG_RAD*RAmin, DEG_RAD*DECmin, - DEG_RAD*RAmax, DEG_RAD*DECmax); + psLogMsg ("psastro", 2, "no sources available for astrometry\n"); + psFree (view); + return true; + } + + psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", + DEG_RAD*RAmin, DEG_RAD*DECmin, + DEG_RAD*RAmax, DEG_RAD*DECmax); psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN", PS_META_REPLACE, "", RAmin); @@ -216,13 +216,13 @@ pmHDU *hdu = pmFPAviewThisHDU (view, fpa); if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } else { - if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - return false; - } + if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } } return true; @@ -238,12 +238,12 @@ // load mosaic-level astrometry? if (phu) { - char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); - if (ctype) { - *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } } if (*bilevelAstrometry) { - pmAstromReadBilevelMosaic (fpa, phu->header); - } + pmAstromReadBilevelMosaic (fpa, phu->header); + } psFree (view); return true; @@ -258,6 +258,6 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; } @@ -290,55 +290,55 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - // XXX we are currently inconsistent with marking the good vs the bad data - // psastroChipAstrom sets data_exists to false if the fit is bad. this is - // probably wrong since it implies there is no data! - - // skip chips for which the astrometry failed (NASTRO == 0) - if (!chip->cells->n) goto skip_chip; - pmCell *cell = chip->cells->data[0]; - if (!cell) goto skip_chip; - - if (!cell->readouts->n) goto skip_chip; - pmReadout *readout = cell->readouts->data[0]; - if (!readout) goto skip_chip; - - psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - if (!updates) goto skip_chip; - - int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); - if (!nAstro) goto skip_chip; - - float astError = psMetadataLookupF32 (&status, updates, "CERROR"); - if (fabs(astError) < 1e-6) goto skip_chip; - - psPlane ptCH, ptFP, ptTP; - psSphere ptSky; - - ptCH.x = 0; - ptCH.y = 0; - psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); - psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); - psDeproject (&ptSky, &ptTP, fpa->toSky); - psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); - - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, ptFP.x); - psVectorAppend (cornerMn, ptFP.y); - psVectorAppend (cornerPn, ptTP.x); - psVectorAppend (cornerQn, ptTP.y); - psVectorAppend (cornerRn, ptSky.r); - psVectorAppend (cornerDn, ptSky.d); - psVectorAppend (cornerMK, 0); - continue; + // XXX we are currently inconsistent with marking the good vs the bad data + // psastroChipAstrom sets data_exists to false if the fit is bad. this is + // probably wrong since it implies there is no data! + + // skip chips for which the astrometry failed (NASTRO == 0) + if (!chip->cells->n) goto skip_chip; + pmCell *cell = chip->cells->data[0]; + if (!cell) goto skip_chip; + + if (!cell->readouts->n) goto skip_chip; + pmReadout *readout = cell->readouts->data[0]; + if (!readout) goto skip_chip; + + psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); + if (!updates) goto skip_chip; + + int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO"); + if (!nAstro) goto skip_chip; + + float astError = psMetadataLookupF32 (&status, updates, "CERROR"); + if (fabs(astError) < 1e-6) goto skip_chip; + + psPlane ptCH, ptFP, ptTP; + psSphere ptSky; + + ptCH.x = 0; + ptCH.y = 0; + psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH); + psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP); + psDeproject (&ptSky, &ptTP, fpa->toSky); + psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d); + + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, ptFP.x); + psVectorAppend (cornerMn, ptFP.y); + psVectorAppend (cornerPn, ptTP.x); + psVectorAppend (cornerQn, ptTP.y); + psVectorAppend (cornerRn, ptSky.r); + psVectorAppend (cornerDn, ptSky.d); + psVectorAppend (cornerMK, 0); + continue; skip_chip: - // new corner locations based on the calibrated astrometry - psVectorAppend (cornerLn, 0.0); - psVectorAppend (cornerMn, 0.0); - psVectorAppend (cornerPn, 0.0); - psVectorAppend (cornerQn, 0.0); - psVectorAppend (cornerRn, 0.0); - psVectorAppend (cornerDn, 0.0); - psVectorAppend (cornerMK, 1); + // new corner locations based on the calibrated astrometry + psVectorAppend (cornerLn, 0.0); + psVectorAppend (cornerMn, 0.0); + psVectorAppend (cornerPn, 0.0); + psVectorAppend (cornerQn, 0.0); + psVectorAppend (cornerRn, 0.0); + psVectorAppend (cornerDn, 0.0); + psVectorAppend (cornerMK, 1); } @@ -349,14 +349,14 @@ for (int i = 0; i < cornerRo->n; i++) { - - psPlane ptTP; - psSphere ptSky; - - ptSky.r = cornerRo->data.F32[i]; - ptSky.d = cornerDo->data.F32[i]; - - psProject (&ptTP, &ptSky, fpa->toSky); - psVectorAppend (cornerPs, ptTP.x); - psVectorAppend (cornerQs, ptTP.y); + + psPlane ptTP; + psSphere ptSky; + + ptSky.r = cornerRo->data.F32[i]; + ptSky.d = cornerDo->data.F32[i]; + + psProject (&ptTP, &ptSky, fpa->toSky); + psVectorAppend (cornerPs, ptTP.x); + psVectorAppend (cornerQs, ptTP.y); } @@ -364,9 +364,9 @@ map->x->coeffMask[1][1] = PS_POLY_MASK_SET; map->y->coeffMask[1][1] = PS_POLY_MASK_SET; - + // fit the valid chips, mask the invalid chips psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs); psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs); - + // apply the linear fit... psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs); @@ -377,5 +377,5 @@ psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf); - psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); + pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd); psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); @@ -387,5 +387,5 @@ float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]); float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]); - + psLogMsg ("psastro", 3, "boresite offset : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]); psLogMsg ("psastro", 3, "boresite angle : %f, scale: %f", angle*PS_DEG_RAD, scale); @@ -395,7 +395,7 @@ psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER"); if (!header) { - header = psMetadataAlloc(); - psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); - psFree (header); // drop this reference + header = psMetadataAlloc(); + psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); + psFree (header); // drop this reference } @@ -408,11 +408,11 @@ if (DEBUG) { - FILE *f = fopen ("corners.dat", "w"); - for (int i = 0; i < cornerRo->n; i++) { - fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", - cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], - cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); - } - fclose (f); + FILE *f = fopen ("corners.dat", "w"); + for (int i = 0; i < cornerRo->n; i++) { + fprintf (f, "%10.6f %10.6f %9.2f %9.2f %9.2f %9.2f | %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f\n", + cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], + cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]); + } + fclose (f); } @@ -437,5 +437,5 @@ psFree (map); psFree (view); - + return true;