Index: trunk/psModules/src/astrom/pmAstrometryTable.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryTable.c (revision 19432) +++ (revision ) @@ -1,754 +1,0 @@ -/** @file pmAstrometryTable.c - * - * @brief Functions to read and write astrometric information using the generic astrometry - * model. - * - * The generic model does not specify the location of the boresite on the sky, and it includes - * a model for the rotator and motion of the boresite. - * - * @ingroup AstroImage - * - * @author EAM, IfA - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-08-08 19:41:56 $ - * - * Copyright 2007 Institute for Astronomy, University of Hawaii - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -/******************************************************************************/ -/* INCLUDE FILES */ -/******************************************************************************/ -#include -#include -#include -#include -#include -#include // for unlink -#include - -#include "pmConfig.h" -#include "pmDetrendDB.h" -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPALevel.h" -#include "pmFPAview.h" -#include "pmFPAfile.h" -#include "pmFPAExtent.h" -#include "pmAstrometryWCS.h" -#include "pmAstrometryUtils.h" -#include "pmAstrometryRegions.h" -#include "pmAstrometryTable.h" - -# define REQUIRE(TEST,MESSAGE){ if (!(TEST)) { psAbort (MESSAGE); }} - -/********************* CheckDataStatus functions *****************************/ - -bool pmAstromCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file) { - - pmFPA *fpa = file->fpa; - - if (view->chip == -1) { - bool exists = pmAstromCheckDataStatusForFPA (fpa); - return exists; - } - if (view->chip >= fpa->chips->n) { - psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n); - return false; - } - pmChip *chip = fpa->chips->data[view->chip]; - - if (view->cell == -1) { - bool exists = pmAstromCheckDataStatusForChip (chip); - return exists; - } - if (view->cell >= chip->cells->n) { - psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n); - return false; - } - psError(PS_ERR_IO, false, "Astrometry only valid at the chip level"); - return false; -} - -bool pmAstromCheckDataStatusForFPA (const pmFPA *fpa) { - - if (!fpa->toTPA) return false; - if (!fpa->fromTPA) return false; - if (!fpa->toSky) return false; - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - if (pmAstromCheckDataStatusForChip (chip)) return true; - } - return false; -} - -bool pmAstromCheckDataStatusForChip (const pmChip *chip) { - - if (!chip->toFPA) return false; - if (!chip->fromFPA) return false; // XXX not strictly needed? - return true; -} - -/********************* Write Data functions *****************************/ - -bool pmAstromWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { - - // write the full table in one pass: require the level to be FPA - if (view->chip != -1) { - psError(PS_ERR_IO, false, "Astrometry must be written at the FPA level"); - return false; - } - - if (!pmAstromWriteFPA (file)) { - psError(PS_ERR_IO, false, "Failed to write Astrometry for fpa"); - return false; - } - return true; -} - -// write out all chip-level Astrometry data for this FPA -bool pmAstromWriteFPA (pmFPAfile *file) { - - if (!pmAstromWritePHU (file)) { - psError(PS_ERR_IO, false, "Failed to write PHU for Astrometry table"); - return false; - } - - if (!pmAstromWriteChips (file)) { - psError(PS_ERR_IO, false, "Failed to write Astrometry for chips"); - return false; - } - - if (!pmAstromWriteFP (file)) { - psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table"); - return false; - } - - if (!pmAstromWriteTP (file)) { - psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table"); - return false; - } - - if (!pmAstromWriteSky (file)) { - psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table"); - return false; - } - - return true; -} - -bool pmAstromWritePHU (pmFPAfile *file) { - - // output header data - psMetadata *outhead = psMetadataAlloc(); - - // use the FPA phu to generate the PHU header - pmHDU *phu = file->fpa->hdu; - - // if there is no FPA PHU, this is a single header+image (extension-less) file. This could be - // the case for an input SPLIT set of files being written out as a MEF. if there is a PHU, - // write it out as a 'blank' - if (phu) { - psMetadataCopy (outhead, phu->header); - } else { - pmConfigConformHeader (outhead, file->format); - } - - psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); - psFitsWriteBlank (file->fits, outhead, ""); - file->wrote_phu = true; - - psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); - psFree (outhead); - - return true; -} - -// fourth layer holds the chips -bool pmAstromWriteChips (pmFPAfile *file) { - - psMetadata *header = psMetadataAlloc(); - psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "CHIPS"); - psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "FOCAL_PLANE"); - psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); - psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "POLYNOMIAL"); - - psArray *table = psArrayAllocEmpty (1); - - pmFPAview *view = pmFPAviewAlloc (0); - - pmChip *chip = NULL; - while ((chip = pmFPAviewNextChip (view, file->fpa, 1)) != NULL) { - - if (!chip->toFPA) continue; - assert (chip->toFPA->x); - assert (chip->toFPA->y); - - psRegion *region = pmChipPixels (chip); - - // set the chip name - char *chiprule = psStringCopy ("{CHIP.NAME}"); - char *chipname = pmFPAfileNameFromRule (chiprule, file, view); - - for (int i = 0; i <= chip->toFPA->x->nX; i++) { - for (int j = 0; j <= chip->toFPA->x->nY; j++) { - psMetadata *row = psMetadataAlloc (); - - psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", chipname); - psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "FOCAL_PLANE"); - psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); - psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); - - psMetadataAddS32(row, PS_LIST_TAIL, "XORDER", PS_META_REPLACE, "", i); - psMetadataAddS32(row, PS_LIST_TAIL, "YORDER", PS_META_REPLACE, "", j); - psMetadataAddS32(row, PS_LIST_TAIL, "NXORDER", PS_META_REPLACE, "", chip->toFPA->x->nX); - psMetadataAddS32(row, PS_LIST_TAIL, "NYORDER", PS_META_REPLACE, "", chip->toFPA->x->nY); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", chip->toFPA->x->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", chip->toFPA->y->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", chip->toFPA->x->coeffErr[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", chip->toFPA->y->coeffErr[i][j]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_X", PS_META_REPLACE, "", chip->toFPA->x->coeffMask[i][j]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_Y", PS_META_REPLACE, "", chip->toFPA->y->coeffMask[i][j]); - psArrayAdd (table, 100, row); - psFree (row); - } - } - psFree (chiprule); - psFree (chipname); - psFree (region); - } - - if (!psFitsWriteTable (file->fits, header, table, "CHIPS")) { - psError(PS_ERR_IO, false, "writing sky data\n"); - psFree(table); - return false; - } - - psFree (view); - psFree (table); - psFree (header); - return true; -} - -// third layer is the focal plane -bool pmAstromWriteFP (pmFPAfile *file) { - - bool status; - - psMetadata *header = psMetadataAlloc(); - psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "FOCAL_PLANE"); - psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); - psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); - psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "POLYNOMIAL"); - - psArray *table = psArrayAllocEmpty (1); - - // XXX is this or the tpa region correct? - psRegion *region = pmAstromFPAExtent (file->fpa); - - // rotate the toTPA to have 0.0 posangle - float posangle = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.POSANGLE"); - - // the to/from TPA transform currently has rotation of posangle; remove it & create the toTPA version - psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, -posangle); - psPlaneTransform *toTPA = psPlaneTransformInvert(NULL, fromTPA, *region, 50); - - for (int i = 0; i <= toTPA->x->nX; i++) { - for (int j = 0; j <= toTPA->x->nY; j++) { - psMetadata *row = psMetadataAlloc (); - psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "FOCAL_PLANE"); - psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "TANGENT_PLANE"); - psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); - psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); - - psMetadataAddS32(row, PS_LIST_TAIL, "XORDER", PS_META_REPLACE, "", i); - psMetadataAddS32(row, PS_LIST_TAIL, "YORDER", PS_META_REPLACE, "", j); - psMetadataAddS32(row, PS_LIST_TAIL, "NXORDER", PS_META_REPLACE, "", toTPA->x->nX); - psMetadataAddS32(row, PS_LIST_TAIL, "NYORDER", PS_META_REPLACE, "", toTPA->x->nY); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", toTPA->x->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", toTPA->y->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", toTPA->x->coeffErr[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", toTPA->y->coeffErr[i][j]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_X", PS_META_REPLACE, "", toTPA->x->coeffMask[i][j]); - psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_Y", PS_META_REPLACE, "", toTPA->y->coeffMask[i][j]); - - psArrayAdd (table, 100, row); - psFree (row); - } - } - - if (!psFitsWriteTable (file->fits, header, table, "FP")) { - psError(PS_ERR_IO, false, "writing sky data\n"); - psFree(table); - psFree (header); - psFree (region); - return false; - } - - psFree (table); - psFree (header); - psFree (region); - return true; -} - -// second layer is the tangent plane -bool pmAstromWriteTP (pmFPAfile *file) { - - bool status; - - // get the boresite model parameters - float Xo = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.X0"); - float Yo = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.Y0"); - float RX = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.RX"); - float RY = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.RY"); - float To = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.T0"); - float Po = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.P0"); - float PosZero = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.POS_ZERO"); /// XXX be consistent with degrees v radians - char *refChip = psMetadataLookupStr (&status, file->fpa->concepts, "FPA.REF.CHIP"); - - psMetadata *header = psMetadataAlloc(); - psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "TANGENT_PLANE"); - psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "SKY"); - psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); - psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "PROJECTION"); - - psArray *table = psArrayAllocEmpty (1); - psMetadata *row = psMetadataAlloc (); - psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "TANGENT_PLANE"); - psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "SKY"); - - psRegion *region = pmAstromFPAExtent (file->fpa); - psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", region->x0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", region->x1); - psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", region->y0); - psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); - - // psMetadataAddF32(row, PS_LIST_TAIL, "XREF", PS_META_REPLACE, "", file->fpa->toSky->R * PM_DEG_RAD); - // psMetadataAddF32(row, PS_LIST_TAIL, "YREF", PS_META_REPLACE, "", file->fpa->toSky->D * PM_DEG_RAD); - - psMetadataAddF32(row, PS_LIST_TAIL, "XSCALE", PS_META_REPLACE, "", file->fpa->toSky->Xs * PM_DEG_RAD); - psMetadataAddF32(row, PS_LIST_TAIL, "YSCALE", PS_META_REPLACE, "", file->fpa->toSky->Ys * PM_DEG_RAD); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_X0", PS_META_REPLACE, "boresite parameter", Xo); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_Y0", PS_META_REPLACE, "boresite parameter", Yo); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_RX", PS_META_REPLACE, "boresite parameter", RX); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_RY", PS_META_REPLACE, "boresite parameter", RY); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_T0", PS_META_REPLACE, "boresite parameter", To); - psMetadataAddF32(row, PS_LIST_TAIL, "BORE_P0", PS_META_REPLACE, "boresite parameter", Po); - psMetadataAddF32(row, PS_LIST_TAIL, "POS_ZERO", PS_META_REPLACE, "boresite parameter", PosZero); - psMetadataAddStr(row, PS_LIST_TAIL, "REF_CHIP", PS_META_REPLACE, "boresite parameter", refChip); - - psArrayAdd (table, 100, row); - psFree (row); - - if (!psFitsWriteTable (file->fits, header, table, "TP")) { - psError(PS_ERR_IO, false, "writing sky data\n"); - psFree (region); - psFree (table); - psFree (header); - return false; - } - - psFree (region); - psFree (table); - psFree (header); - return (true); -} - -// first layer is the sky -bool pmAstromWriteSky (pmFPAfile *file) { - - psMetadata *header = psMetadataAlloc(); - psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "SKY"); - psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "NONE"); - psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "NONE"); - psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "NONE"); - - psArray *table = psArrayAllocEmpty (1); - psMetadata *row = psMetadataAlloc (); - psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", "SKY"); - psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "NONE"); - - psArrayAdd (table, 100, row); - psFree (row); - - if (!psFitsWriteTable (file->fits, header, table, "SKY")) { - psError(PS_ERR_IO, false, "writing sky data\n"); - psFree(table); - return false; - } - - psFree (table); - psFree (header); - return (true); -} - -/********************* Read Data functions *****************************/ - -bool pmAstromReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) - { - - // write the full table in one pass: require the level to be FPA - if (view->chip != -1) { - psError(PS_ERR_IO, false, "Astrometry must be read at the FPA level"); - return false; - } - - if (!pmAstromReadFPA (file)) { - psError(PS_ERR_IO, false, "Failed to read Astrometry for fpa"); - return false; - } - return true; - } - -// read out all chip-level Astrometry data for this FPA -bool pmAstromReadFPA (pmFPAfile *file) { - - if (!pmAstromReadPHU (file)) { - psError(PS_ERR_IO, false, "Failed to read PHU for Astrometry table"); - return false; - } - - if (!pmAstromReadChips (file)) { - psError(PS_ERR_IO, false, "Failed to read Astrometry for chips"); - return false; - } - - if (!pmAstromReadFP (file)) { - psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table"); - return false; - } - - // NOTE : TP must come after FP as it applies the POS, ROT boresite corrections to the - // transformation determined in FP - if (!pmAstromReadTP (file)) { - psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table"); - return false; - } - - if (!pmAstromReadSky (file)) { - psError(PS_ERR_IO, false, "Failed to read Sky for Astrometry table"); - return false; - } - - return true; -} - -bool pmAstromReadPHU (pmFPAfile *file) { - - fprintf (stderr, "not sure we need to read PHU\n"); - return true; -} - -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (i); - } - return -1; -} - -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (chip); - } - return NULL; -} - -// first layer converts Chip to Focal Plane -bool pmAstromReadChips (pmFPAfile *file) { - - bool status; - - // set FITS cursor - if (!psFitsMoveExtName (file->fits, "CHIPS")) { - psError(PS_ERR_IO, false, "missing CHIPS extension in astrometry table\n"); - return false; - } - - // free exising tranformations in prep for new alloc below - for (int i = 0; i < file->fpa->chips->n; i++) { - pmChip *chip = file->fpa->chips->data[i]; - psFree (chip->toFPA); - chip->toFPA = NULL; - } - - // XXX do I need anything from the header? - psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header - if (!header) psAbort("cannot read table header"); - - // load the full table in one shot - psArray *table = psFitsReadTable (file->fits); - if (!table) psAbort("cannot read table"); - fprintf (stderr, "read %ld rows from FP\n", table->n); - - // parse the table entries - for (int i = 0; i < table->n; i++) { - psMetadata *row = table->data[i]; - - // name of the chip for this row. - char *chipname = psMetadataLookupStr (&status, row, "SEGMENT"); - - // get chip from name - pmChip *chip = pmConceptsChipFromName (file->fpa, chipname); - REQUIRE (chip, "invalid chip name"); - - // define the toFPA transform if not already defined - int nX = psMetadataLookupS32(&status, row, "NXORDER"); REQUIRE (status, "missing NXORDER"); - int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER"); - if (chip->toFPA == NULL) { - chip->toFPA = psPlaneTransformAlloc(nX, nY); - } else { - REQUIRE (chip->toFPA->x->nX == nX, "mismatch in chip order"); - REQUIRE (chip->toFPA->x->nY == nY, "mismatch in chip order"); - REQUIRE (chip->toFPA->y->nX == nX, "mismatch in chip order"); - REQUIRE (chip->toFPA->y->nY == nY, "mismatch in chip order"); - } - - int ix = psMetadataLookupS32(&status, row, "XORDER"); REQUIRE (status, "missing XORDER"); - int iy = psMetadataLookupS32(&status, row, "YORDER"); REQUIRE (status, "missing YORDER"); - - // XXX do I need to use this or can i rely on the camera concepts? - // minX = psMetadataLookupF32(&status, row, "MINX"); - // maxX = psMetadataLookupF32(&status, row, "MAXX"); - // minY = psMetadataLookupF32(&status, row, "MINY"); - // maxY = psMetadataLookupF32(&status, row, "MAXY"); - - // XXX need to include mask values - chip->toFPA->x->coeff[ix][iy] = psMetadataLookupF32(&status, row, "POLY_X"); - chip->toFPA->y->coeff[ix][iy] = psMetadataLookupF32(&status, row, "POLY_Y"); - chip->toFPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X"); - chip->toFPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y"); - chip->toFPA->x->coeffMask[ix][iy] = psMetadataLookupU8(&status, row, "MASK_X"); - chip->toFPA->y->coeffMask[ix][iy] = psMetadataLookupU8(&status, row, "MASK_Y"); - } - - // convert the toFPA transfomations to fromFPA transformations - for (int i = 0; i < file->fpa->chips->n; i++) { - pmChip *chip = file->fpa->chips->data[i]; - if (!chip->toFPA) continue; - psRegion *region = pmChipPixels (chip); - psFree (chip->fromFPA); - chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); - psFree (region); - } - - psFree (table); - psFree (header); - return true; -} - -// second layer converts Focal Plane to Tangent Plane (unrotated) -bool pmAstromReadFP (pmFPAfile *file) { - - bool status; - - if (!psFitsMoveExtName (file->fits, "FP")) { - psError(PS_ERR_IO, false, "missing FP extension in astrometry table\n"); - return false; - } - - psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header - if (!header) psAbort("cannot read table header"); - - // free the old - psFree (file->fpa->toTPA); - file->fpa->toTPA = NULL; - - // read the complete table data at one shot - psArray *table = psFitsReadTable (file->fits); - fprintf (stderr, "read %ld rows from FP\n", table->n); - - // parse the table - for (int i = 0; i < table->n; i++) { - psMetadata *row = table->data[i]; - - // there is only one transformation in this table; the order is defined in the header - int nX = psMetadataLookupS32(&status, row, "NXORDER"); REQUIRE (status, "missing NXORDER"); - int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER"); - if (file->fpa->toTPA == NULL) { - // allocate the new transformation - file->fpa->toTPA = psPlaneTransformAlloc(nX, nY); - } else { - REQUIRE (file->fpa->toTPA->x->nX == nX, "mismatch in chip order"); - REQUIRE (file->fpa->toTPA->x->nY == nY, "mismatch in chip order"); - REQUIRE (file->fpa->toTPA->y->nX == nX, "mismatch in chip order"); - REQUIRE (file->fpa->toTPA->y->nY == nY, "mismatch in chip order"); - } - - int ix = psMetadataLookupS32(&status, row, "XORDER"); REQUIRE (status, "missing XORDER"); - int iy = psMetadataLookupS32(&status, row, "YORDER"); REQUIRE (status, "missing YORDER"); - file->fpa->toTPA->x->coeff[ix][iy] = psMetadataLookupF32(&status, row, "POLY_X"); REQUIRE (status, "missing POLY_X"); - file->fpa->toTPA->y->coeff[ix][iy] = psMetadataLookupF32(&status, row, "POLY_Y"); REQUIRE (status, "missing POLY_Y"); - file->fpa->toTPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X"); REQUIRE (status, "missing ERROR_X"); - file->fpa->toTPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y"); REQUIRE (status, "missing ERROR_Y"); - file->fpa->toTPA->x->coeffMask[ix][iy] = psMetadataLookupU8 (&status, row, "MASK_X"); REQUIRE (status, "missing MASK_X"); - file->fpa->toTPA->y->coeffMask[ix][iy] = psMetadataLookupU8 (&status, row, "MASK_Y"); REQUIRE (status, "missing MASK_Y"); - } - - psRegion *region = pmAstromFPAExtent (file->fpa); - psFree (file->fpa->fromTPA); - file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50); - - psFree (table); - psFree (header); - psFree (region); - return true; -} - -# define TRANSFER(TO,FROM,NAME) { \ - psMetadataItem *item = psMetadataLookup(FROM,NAME); \ - if (!item) psAbort ("cannot find %s", NAME); \ - psMetadataItem *newItem = psMetadataItemCopy(item); \ - if (!psMetadataAddItem(TO, newItem, PS_LIST_TAIL, PS_META_REPLACE)) { \ - psAbort ("cannot copy %s", NAME); \ - } \ - psFree (newItem); } - -// third layer applies boresite corrections and converts tangent plane to sky -bool pmAstromReadTP (pmFPAfile *file) { - - if (!psFitsMoveExtName (file->fits, "TP")) { - psError(PS_ERR_IO, false, "missing TP extension in astrometry table\n"); - return false; - } - - psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header - if (!header) psAbort("cannot read table header"); - - psArray *table = psFitsReadTable (file->fits); - fprintf (stderr, "read %ld rows from TP\n", table->n); - if (table->n != 1) psAbort("invalid number of rows in TP table (%ld)", table->n); - - psMetadata *row = table->data[0]; - - // move needed items to the concepts - TRANSFER (file->fpa->concepts, row, "XSCALE"); - TRANSFER (file->fpa->concepts, row, "XSCALE"); - TRANSFER (file->fpa->concepts, row, "YSCALE"); - TRANSFER (file->fpa->concepts, row, "BORE_X0"); - TRANSFER (file->fpa->concepts, row, "BORE_Y0"); - TRANSFER (file->fpa->concepts, row, "BORE_RX"); - TRANSFER (file->fpa->concepts, row, "BORE_RY"); - TRANSFER (file->fpa->concepts, row, "BORE_T0"); - TRANSFER (file->fpa->concepts, row, "BORE_P0"); - TRANSFER (file->fpa->concepts, row, "POS_ZERO"); - TRANSFER (file->fpa->concepts, row, "REF_CHIP"); - - psFree (table); - psFree (header); - return (true); -} - -// first layer is the sky -bool pmAstromReadSky (pmFPAfile *file) { - - if (!psFitsMoveExtName (file->fits, "SKY")) { - psError(PS_ERR_IO, false, "missing SKY extension in astrometry table\n"); - return false; - } - - psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header - if (!header) psAbort("cannot read table header"); - - psArray *table = psFitsReadTable (file->fits); - fprintf (stderr, "read %ld rows from SKY\n", table->n); - if (table->n != 1) psAbort("invalid number of rows in SKY table (%ld)", table->n); - - // XXX not much information of interest in this table... - - psFree (table); - psFree (header); - return (true); -} - -// third layer applies boresite corrections and converts tangent plane to sky -bool pmAstromSetTP (pmFPAfile *file, psMetadata *concepts) { - - bool status; - - // these externally supplied values are used to set the final transformation terms - double RA = psMetadataLookupF64 (&status, concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA"); - double DEC = psMetadataLookupF64 (&status, concepts, "FPA.DEC"); REQUIRE (status, "missing FPA.DEC"); - double POS = PM_RAD_DEG * psMetadataLookupF64 (&status, concepts, "FPA.POSANGLE"); REQUIRE (status, "missing FPA.POSANGLE"); - // double ROT = psMetadataLookupF64 (&status, concepts, "FPA.ROTANGLE"); REQUIRE (status, "missing FPA.ROTANGLE"); - - // get projection scale; center is supplied - float Xs = psMetadataLookupF32(&status, file->fpa->concepts, "XSCALE") * PM_RAD_DEG; REQUIRE (status, "missing XSCALE"); - float Ys = psMetadataLookupF32(&status, file->fpa->concepts, "YSCALE") * PM_RAD_DEG; REQUIRE (status, "missing YSCALE"); - - // allocate a new toSky projection using the reported position - psFree (file->fpa->toSky); - file->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_DIS); - - // get boresite correction terms. RX,RY,To,Po define an ellipse along which the boresite travels - double Xo = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_X0"); REQUIRE (status, "missing "); - double Yo = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_Y0"); REQUIRE (status, "missing "); - double RX = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_RX"); REQUIRE (status, "missing "); - double RY = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_RY"); REQUIRE (status, "missing "); - double To = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_T0"); REQUIRE (status, "missing "); - double Po = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_P0"); REQUIRE (status, "missing "); - - // the true rotation of the instrument is POSANGLE - POS_ZERO - double PosZero = PM_RAD_DEG * psMetadataLookupF32(&status, file->fpa->concepts, "POS_ZERO"); REQUIRE (status, "missing POS_ZERO"); - char *refChip = psMetadataLookupStr(&status, file->fpa->concepts, "REF_CHIP"); REQUIRE (status, "missing REF_CHIP"); - - // apply true posangle = -(POS - POS_ZERO) - psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, (POS - PosZero)); - psFree (file->fpa->fromTPA); - file->fpa->fromTPA = fromTPA; - - psRegion *region = pmAstromFPAExtent (file->fpa); - psFree (file->fpa->toTPA); - file->fpa->toTPA = psPlaneTransformInvert(NULL, file->fpa->fromTPA, *region, 50); - - // current position of the nominal boresite in refChip coordinates - double X = Xo + RX*cos(POS - To)*cos(Po) + RY*sin(POS - To)*sin(Po); - double Y = Yo + RY*sin(POS - To)*cos(Po) - RX*cos(POS - To)*sin(Po); - - // get reference chip from name - pmChip *chip = pmConceptsChipFromName (file->fpa, refChip); - if (!chip) psAbort ("invalid chip name for reference"); - - psPlane *boreCH = psPlaneAlloc(); - psPlane *boreFP = psPlaneAlloc(); - psPlane *boreTP = psPlaneAlloc(); - psSphere *boreSky = psSphereAlloc(); - - // find the true RA,DEC coord of the mirror of the reported boresite - boreCH->x = -X; - boreCH->y = -Y; - psPlaneTransformApply (boreFP, chip->toFPA, boreCH); - psPlaneTransformApply (boreTP, file->fpa->toTPA, boreFP); - psDeproject (boreSky, boreTP, file->fpa->toSky); - - // modify the projection to account for offset between true and reported boresite - file->fpa->toSky->R = boreSky->r; - file->fpa->toSky->D = boreSky->d; - - psFree (boreCH); - psFree (boreFP); - psFree (boreTP); - psFree (boreSky); - - psFree (region); - return (true); -} - Index: trunk/psModules/src/astrom/pmAstrometryTable.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryTable.h (revision 19432) +++ (revision ) @@ -1,43 +1,0 @@ -/* @file pmAstrometryTable.h - * @brief Astrometry table I/O functions - * - * @author EAM, IfA - * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-22 17:53:27 $ - * Copyright 2006 Institute for Astronomy, University of Hawaii - */ - -#ifndef PM_ASTROMETRY_TABLE_H -#define PM_ASTROMETRY_TABLE_H - -/// @addtogroup Astrometry -/// @{ - -bool pmAstromCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file); -bool pmAstromCheckDataStatusForFPA (const pmFPA *fpa); -bool pmAstromCheckDataStatusForChip (const pmChip *chip); - -bool pmAstromWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmAstromWriteFPA (pmFPAfile *file); -bool pmAstromWritePHU (pmFPAfile *file); -bool pmAstromWriteSky (pmFPAfile *file); -bool pmAstromWriteTP (pmFPAfile *file); -bool pmAstromWriteFP (pmFPAfile *file); -bool pmAstromWriteChips (pmFPAfile *file); - -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); - -bool pmAstromReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); -bool pmAstromReadFPA (pmFPAfile *file); -bool pmAstromReadPHU (pmFPAfile *file); -bool pmAstromReadChips (pmFPAfile *file); -bool pmAstromReadFP (pmFPAfile *file); -bool pmAstromReadTP (pmFPAfile *file); -bool pmAstromReadSky (pmFPAfile *file); - -bool pmAstromSetTP (pmFPAfile *file, psMetadata *concepts); - -/// @} -#endif