Index: trunk/psModules/src/astrom/pmAstrometryTable.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15884) +++ trunk/psModules/src/astrom/pmAstrometryTable.c (revision 15897) @@ -10,6 +10,6 @@ * * @author EAM, IfA - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-19 18:57:05 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-22 17:53:15 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -120,4 +120,19 @@ } + 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"); @@ -125,18 +140,4 @@ } - if (!pmAstromWriteTP (file)) { - psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table"); - return false; - } - - if (!pmAstromWriteFP (file)) { - psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table"); - return false; - } - - if (!pmAstromWriteChips (file)) { - psError(PS_ERR_IO, false, "Failed to write Astrometry for chips"); - return false; - } return true; } @@ -174,122 +175,4 @@ psFree (outhead); - 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); -} - -// second layer is the tangent plane -bool pmAstromWriteTP (pmFPAfile *file) { - - 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, "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, "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); - - 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); -} - -// third layer is the focal plane -bool pmAstromWriteFP (pmFPAfile *file) { - - 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); - - for (int i = 0; i <= file->fpa->toTPA->x->nX; i++) { - for (int j = 0; j <= file->fpa->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, "NX", PS_META_REPLACE, "", file->fpa->toTPA->x->nX); - psMetadataAddS32(row, PS_LIST_TAIL, "NY", PS_META_REPLACE, "", file->fpa->toTPA->x->nY); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", file->fpa->toTPA->x->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", file->fpa->toTPA->y->coeff[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", file->fpa->toTPA->x->coeffErr[i][j]); - psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", file->fpa->toTPA->y->coeffErr[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); - return false; - } - - psFree (table); - psFree (header); - psFree (region); return true; } @@ -332,10 +215,14 @@ psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", region->y1); - psMetadataAddS32(row, PS_LIST_TAIL, "NX", PS_META_REPLACE, "", i); - psMetadataAddS32(row, PS_LIST_TAIL, "NY", PS_META_REPLACE, "", j); + 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); @@ -357,4 +244,158 @@ 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); } @@ -504,4 +545,6 @@ 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"); } @@ -509,4 +552,5 @@ 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); @@ -533,16 +577,7 @@ if (!header) psAbort("cannot read table header"); - // there is only one transformation in this table; the order is defined in the header - int nX = psMetadataLookupS32(&status, header, "NXORDER"); REQUIRE (status, "missing NXORDER"); - int nY = psMetadataLookupS32(&status, header, "NYORDER"); REQUIRE (status, "missing NYORDER"); - - // allocate the new transformation + // free the old psFree (file->fpa->toTPA); - file->fpa->toTPA = psPlaneTransformAlloc(nX, nY); - - // free & DON'T allocate the new transformation fromTPA; this is created after the boresite - // is applied in pmAstromReadTP - psFree (file->fpa->fromTPA); - file->fpa->fromTPA = NULL; + file->fpa->toTPA = NULL; // read the complete table data at one shot @@ -553,28 +588,50 @@ 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->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); - return true; -} + 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) { - bool status; - - // these externally supplied values are used to set the final transformation terms - double RA = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA"); - double DEC = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.DEC"); REQUIRE (status, "missing FPA.DEC"); - double POS = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.POSANGLE"); REQUIRE (status, "missing FPA.POSANGLE"); - double ROT = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.ROTANGLE"); REQUIRE (status, "missing FPA.ROTANGLE"); - if (!psFitsMoveExtName (file->fits, "TP")) { psError(PS_ERR_IO, false, "missing TP extension in astrometry table\n"); @@ -591,38 +648,17 @@ psMetadata *row = table->data[0]; - // get projection scale; center is supplied - float Xs = psMetadataLookupF32(&status, row, "XSCALE") * PM_RAD_DEG; REQUIRE (status, "missing XSCALE"); - float Ys = psMetadataLookupF32(&status, row, "YSCALE") * PM_RAD_DEG; REQUIRE (status, "missing YSCALE"); - - // allocate a new toSky projection - psFree (file->fpa->toSky); - file->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_WRP); - - // get boresite correction terms. L,M,T define an ellipse along which the boresite travels - double R_L = psMetadataLookupF32(&status, row, "ROTATOR_L"); REQUIRE (status, "missing ROTATOR_L"); - double R_M = psMetadataLookupF32(&status, row, "ROTATOR_M"); REQUIRE (status, "missing ROTATOR_M"); - double R_T = psMetadataLookupF32(&status, row, "ROTATOR_T"); REQUIRE (status, "missing ROTATOR_T"); - - // the position of the boresite on the ellipse is the ROTANGLE - ROT_ZERO - double R_0 = psMetadataLookupF32(&status, row, "ROT_ZERO"); REQUIRE (status, "missing ROT_ZERO"); - - // find the current position of the boresite for this image - double Lo = R_L*cos(ROT - R_0)*cos(R_T) + R_M*sin(ROT - R_0)*sin(R_T); - double Mo = R_M*sin(ROT - R_0)*cos(R_T) - R_L*cos(ROT - R_0)*sin(R_T); - - // apply new reference pixel to transformation - psPlaneTransform *tmpToTPA = psPlaneTransformSetCenter (NULL, file->fpa->toTPA, Lo, Mo); - psFree (file->fpa->toTPA); - - // apply POSANGLE - double Po = psMetadataLookupF32(&status, row, "POS_ZERO"); REQUIRE (status, "missing POS_ZERO"); - file->fpa->toTPA = psPlaneTransformRotate (NULL, tmpToTPA, POS - Po); - psFree (tmpToTPA); - - psRegion *region = pmAstromFPAExtent (file->fpa); - psFree (file->fpa->fromTPA); - file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50); - - psFree (region); + // 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); @@ -652,2 +688,75 @@ } +// 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_WRP); + + // 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); +} +