Index: trunk/psastro/src/psastroBuildFPA.c =================================================================== --- trunk/psastro/src/psastroBuildFPA.c (revision 5560) +++ trunk/psastro/src/psastroBuildFPA.c (revision 5575) @@ -20,4 +20,13 @@ pmFPA *fpa = pmFPAAlloc (NULL, NULL); + // identify matrix for fpa to tpa + fpa->toTPA = psPlaneDistortAlloc (1, 1, 0, 0); + fpa->toTPA->x->coeff[1][0][0][0] = 1; + fpa->toTPA->x->mask[1][1][0][0] = 1; + + fpa->toTPA->y->coeff[0][1][0][0] = 1; + fpa->toTPA->y->mask[1][1][0][0] = 1; + fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane); + psArrayRealloc (fpa->chips, Nchips); for (int i = 0; i < Nchips; i++) { @@ -30,5 +39,13 @@ // XXX EAM : extend this function to take more than one header cell->header = header; - pmCellInterpretWCS (cell, header); + pmAstromReadWCS (&chip->toFPA, &fpa->toSky, header); + chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); + + // allocate identity matrix for cell to chip + cell->toChip = psPlaneTransformAlloc (1, 1); + cell->toChip->x->coeff[1][0] = 1; + cell->toChip->x->mask[1][1] = 1; + cell->toChip->y->coeff[0][1] = 1; + cell->toChip->y->mask[1][1] = 1; psArrayRealloc (cell->readouts, Nreadouts); @@ -52,11 +69,11 @@ } -// interpret header WCS -bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) { - - pmChip *chip; - pmFPA *fpa; +// interpret header WCS (only handles traditional WCS for the moment) +bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) { + + psPlaneTransform *toFPA; + psProjection *toSky; + psProjectionType type; bool status; - psProjectionType type; float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; float pc1_1, pc1_2, pc2_1, pc2_2; @@ -136,22 +153,4 @@ got_matrix: - // XXX EAM : these must already be set - chip = cell->parent; - fpa = chip->parent; - - // set toChip to identity as default - // XXX EAM : can we actually use higher order? - int nX = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NX"); - if (!status) nX = 1; - int nY = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NY"); - if (!status) nY = 1; - - cell->toChip = psPlaneTransformAlloc (1, 1); - cell->toChip->x->coeff[1][0] = 1; - cell->toChip->x->mask[1][1] = 1; - - cell->toChip->y->coeff[0][1] = 1; - cell->toChip->y->mask[1][1] = 1; - // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the // toFPA must be modified to match the crval(i), scale(i) and crpix(i) @@ -161,34 +160,67 @@ // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) // XXX EAM : I've fixed this in pslib eam_rel8_b2 - chip->toFPA = psPlaneTransformAlloc (1, 1); + toFPA = psPlaneTransformAlloc (1, 1); - chip->toFPA->x->coeff[0][0] = crpix1; - chip->toFPA->x->coeff[1][0] = pc1_1; - chip->toFPA->x->coeff[0][1] = pc1_2; - chip->toFPA->x->mask[1][1] = 1; - - chip->toFPA->y->coeff[0][0] = crpix2; - chip->toFPA->y->coeff[1][0] = pc2_1; - chip->toFPA->y->coeff[0][1] = pc2_2; - chip->toFPA->y->mask[1][1] = 1; - - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); - - // set toTPA to identity as default - // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581) - if (fpa->toTPA == NULL) { - fpa->toTPA = psPlaneDistortAlloc (1, 1, 0, 0); - fpa->toTPA->x->coeff[1][0][0][0] = 1; - fpa->toTPA->x->mask[1][1][0][0] = 1; - - fpa->toTPA->y->coeff[0][1][0][0] = 1; - fpa->toTPA->y->mask[1][1][0][0] = 1; - fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane); - } else { - psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n"); - } + toFPA->x->coeff[0][0] = crpix1; + toFPA->x->coeff[1][0] = pc1_1; + toFPA->x->coeff[0][1] = pc1_2; + toFPA->x->mask[1][1] = 1; + + toFPA->y->coeff[0][0] = crpix2; + toFPA->y->coeff[1][0] = pc2_1; + toFPA->y->coeff[0][1] = pc2_2; + toFPA->y->mask[1][1] = 1; // center of projection is (0,0) coordinate of TPA - fpa->toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); + toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); + + *toFPAOut = toFPA; + *toSkyOut = toSky; + return true; } + + +// convert toFPA / toSky components to traditional WCS +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) { + + switch (toSky->type) { + case PS_PROJ_SIN: + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN"); + break; + case PS_PROJ_TAN: + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN"); + break; + case PS_PROJ_AIT: + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT"); + break; + case PS_PROJ_PAR: + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR"); + break; + default: + psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", toSky->type); + return false; + } + + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD); + + psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]); + + // alternative representations use + // CD1_1 = PC001001*CDELT1, etc + // make these representations optional + + return true; +}