Index: /trunk/psastro/doc/mktest.txt =================================================================== --- /trunk/psastro/doc/mktest.txt (revision 5574) +++ /trunk/psastro/doc/mktest.txt (revision 5575) @@ -22,2 +22,10 @@ within floating point errors. +examples: + +build a fake dataset +# psastro-mktest doc/psastro-mktest.conf test.cmp test.cat test.raw test.ref + +run psastrom on the fake data: +# psastro doc/psastro.conf test.cmp test.dat + Index: /trunk/psastro/src/psastro-mktest.c =================================================================== --- /trunk/psastro/src/psastro-mktest.c (revision 5574) +++ /trunk/psastro/src/psastro-mktest.c (revision 5575) @@ -80,9 +80,4 @@ chip->toFPA->y->coeff[1][1] = 0; - psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]); - psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]); - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); } @@ -136,12 +131,16 @@ // write output header WCS info - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "foo", RA*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "foo", DEC*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN"); - psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN"); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.x); - psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.y); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN"); + psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN"); + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "foo", RA*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "foo", DEC*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.x); + psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "foo", tmp1.y); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "foo", PS*DEG_RAD); + psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]); + psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]); } Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 5574) +++ /trunk/psastro/src/psastro.h (revision 5575) @@ -11,5 +11,7 @@ # define toSky projection -bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) ; +bool pmAstromReadWCS (psPlaneTransform **toFPA, psProjection **toSky, psMetadata *header); +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header); + pmFPA *pmFPACopyAstrom (pmFPA *inFPA); psMetadata *psastroArguments (int *argc, char **argv); @@ -21,13 +23,13 @@ bool psastroWriteCMP (pmFPA *fpa, char *filename); -psMetadata *testArguments (int *argc, char **argv); -psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout); -psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout); -psPlane* psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell); +psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout); +psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout); +psPlane *psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell); -bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); -bool testWriteRef (char *filename, psArray *sources); -bool testWriteRaw (char *filename, psArray *sources); +bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); +bool testWriteRef (char *filename, psArray *sources); +bool testWriteRaw (char *filename, psArray *sources); +psMetadata *testArguments (int *argc, char **argv); bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky); @@ -35,10 +37,4 @@ bool psastroProjectRefstars (psArray *stars, pmReadout *readout); -psPlane* psCoordChipToCell_EAM( - psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created - const psPlane* in, ///< the Chip coordinate - const pmCell* cell ///< the cell of interest -); - psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform); -pmFPA *pmFPACopyAstrom (pmFPA *inFPA); +pmFPA *pmFPACopyAstrom (pmFPA *inFPA); Index: /trunk/psastro/src/psastroBuildFPA.c =================================================================== --- /trunk/psastro/src/psastroBuildFPA.c (revision 5574) +++ /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; +} Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 5574) +++ /trunk/psastro/src/psastroUtils.c (revision 5575) @@ -102,7 +102,9 @@ match = pmAstromRadiusMatch (rawstars, refstars, config); - // fit astrometric terms + // improved fit for astrometric terms pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config); + chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); } + pmAstromWriteWCS (chip->toFPA, fpa->toSky, cell->header); } }