Index: /trunk/psastro/doc/notes.txt =================================================================== --- /trunk/psastro/doc/notes.txt (revision 6175) +++ /trunk/psastro/doc/notes.txt (revision 6176) @@ -10,4 +10,5 @@ - in the first stage, the per-chip results are applied to the chip.toFPA parameter set. + - in the second stage, the chip results are collectively used to modify the fpa.toSky terms. Index: /trunk/psastro/src/psastroBuildFPA.c =================================================================== --- /trunk/psastro/src/psastroBuildFPA.c (revision 6175) +++ /trunk/psastro/src/psastroBuildFPA.c (revision 6176) @@ -68,159 +68,2 @@ return (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; - float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; - float pc1_1, pc1_2, pc2_1, pc2_2; - - // *** interpret header data, convert to crval(i), etc - char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2"); - if (!status) { - psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n"); - return false; - } - - // determine projection type - type = PS_PROJ_NTYPE; - if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN; - if (!strcmp (&ctype[4], "-TAN")) type = PS_PROJ_TAN; - if (!strcmp (&ctype[4], "-AIT")) type = PS_PROJ_AIT; - if (!strcmp (&ctype[4], "-PAR")) type = PS_PROJ_PAR; - if (type == PS_PROJ_NTYPE) { - psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); - return false; - } - - crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); - crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); - crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); - crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); - - // test the CDELTi varient - cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); - if (status) { - cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); - - // test the CROTAi varient: - double rotate = psMetadataLookupF32 (&status, header, "CROTA2"); - if (status) { - double Lambda = cdelt2 / cdelt1; - pc1_1 = cos(rotate*RAD_DEG); - pc1_2 = -sin(rotate*RAD_DEG) * Lambda; - pc2_1 = sin(rotate*RAD_DEG) / Lambda; - pc2_2 = cos(rotate*RAD_DEG); - goto got_matrix; - } - - // test the PC00i00j varient: - pc1_1 = psMetadataLookupF32 (&status, header, "PC001001"); - if (status) { - pc1_2 = psMetadataLookupF32 (&status, header, "PC001002"); - pc2_1 = psMetadataLookupF32 (&status, header, "PC002001"); - pc2_2 = psMetadataLookupF32 (&status, header, "PC002002"); - - // XXX EAM : add Elixir polynomial terms here eventually - goto got_matrix; - } - psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; - } - - // test the CDi_j varient - pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1"); - if (status) { - pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2"); - pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1"); - pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2"); - - // renormalize to cdelt1, cdelt2, etc - float scale = hypot (pc1_1, pc1_2); - cdelt1 = cdelt2 = scale; - pc1_1 /= scale; - pc1_2 /= scale; - pc2_1 /= scale; - pc2_2 /= scale; - goto got_matrix; - } - psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); - return false; - -got_matrix: - - // 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) - // scale = scale(i)/scale(0) (i == chip #) - // project crval1(0),crval2(0 using projection - - // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) - // XXX EAM : I've fixed this in pslib eam_rel8_b2 - toFPA = psPlaneTransformAlloc (1, 1); - - 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 - 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/psastroChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroChipAstrom.c (revision 6176) +++ /trunk/psastro/src/psastroChipAstrom.c (revision 6176) @@ -0,0 +1,73 @@ +# include "psastro.h" + +// measure per-chip astrometry terms +// this function requires a chip from a well-formed pmFPA +// with PSASTRO.OBJECTS on the readouts +bool psastroChipAstrom (pmChip *chip, psMetadata *config) { + + bool status; + psArray *match; + pmAstromStats stats; + + // cells->n > 1 is not well-defined + if (chip->cells->n > 1) { + psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1"); + return false; + } + pmCell *cell = chip->cells->data[0]; + + // readouts->n > 1 is not well-defined + if (cell->readouts->n > 1) { + psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1"); + return false; + } + pmReadout *readout = cell->readouts->data[0]; + psMetadata *header = pmReadoutGetHeader (readout); + + // load the corresponding reference data (DVO command) + // XXX EAM : this needs to be improved, perhaps moved to a higher level + psArray *refstars = psastroLoadReference (header, config); + + // we require PASTRO.OBJECTS (pmAstromObj) on the analysis for this readout + psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.OBJECTS"); + + // project the rawstars to the current best guess astrometry + psastroProjectRawstars (rawstars, readout); + + // use the header & config info to project refstars onto the focal plane + psastroProjectRefstars (refstars, readout); + + // testWriteRaw ("ref.inp", refstars); + // testWriteRaw ("raw.inp", rawstars); + + // fprintf (stderr, "rawstars:\n"); + // psastroDumpStars (rawstars); + // fprintf (stderr, "refstars:\n"); + // psastroDumpStars (refstars); + + // find initial offset / rotation + stats = pmAstromGridMatch (rawstars, refstars, config); + + // adjust the chip.toFPA terms only + pmAstromGridApply (chip->toFPA, stats); + chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); + + // use fit result to re-project rawstars + psastroProjectRawstars (rawstars, readout); + psastroProjectRefstars (refstars, readout); + + // testWriteRaw ("ref.dat", refstars); + // testWriteRaw ("raw.dat", rawstars); + + // use small radius to match stars + match = pmAstromRadiusMatch (rawstars, refstars, config); + + // improved fit for astrometric terms + pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config); + chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); + + // per-chip astrometry has been updated + pmAstromWriteWCS (chip->toFPA, fpa->toSky, header); + + return true; +} Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 6176) +++ /trunk/psastro/src/psastroConvert.c (revision 6176) @@ -0,0 +1,7 @@ +# include "psastro.h" + +psArray *pmSourceToAstromObj (psArray *sources) { + + + +} Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 6175) +++ /trunk/psastro/src/psastroUtils.c (revision 6176) @@ -32,10 +32,5 @@ } -// measure per-chip astrometry terms -bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) { - - bool status; - psArray *match; - pmAstromStats stats; +bool psastroNormFPA (pmFPA *fpa, psMetadata *config) { // save the raw astrometry for later reference @@ -45,67 +40,5 @@ for (int i = 0; i < fpa->chips->n; i++) { pmChip *chip = fpa->chips->data[i]; - - // cells->n > 1 is not yet well-defined - if (chip->cells->n > 1) { - psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1"); - return false; - } - - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; - - // readouts->n > 1 is not yet well-defined - if (cell->readouts->n > 1) { - psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1"); - return false; - } - - // load the corresponding reference data (DVO command) - psArray *refstars = psastroLoadReference (cell->header, config); - - for (int k = 0; k < cell->readouts->n; k++) { - - pmReadout *readout = cell->readouts->data[k]; - - // pull out the SUBSET rawstars (a view) - psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "STARS.SUBSET"); - - // project the rawstars to the current best guess astrometry - psastroProjectRawstars (rawstars, readout); - - // use the header & config info to project refstars onto the focal plane - psastroProjectRefstars (refstars, readout); - - // testWriteRaw ("ref.inp", refstars); - // testWriteRaw ("raw.inp", rawstars); - - // fprintf (stderr, "rawstars:\n"); - // psastroDumpStars (rawstars); - // fprintf (stderr, "refstars:\n"); - // psastroDumpStars (refstars); - - // find initial offset / rotation - stats = pmAstromGridMatch (rawstars, refstars, config); - - // adjust the chip.toFPA terms only - pmAstromGridApply (chip->toFPA, stats); - chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); - - // use fit result to re-project rawstars - psastroProjectRawstars (rawstars, readout); - psastroProjectRefstars (refstars, readout); - - // testWriteRaw ("ref.dat", refstars); - // testWriteRaw ("raw.dat", rawstars); - - // use small radius to match stars - match = pmAstromRadiusMatch (rawstars, refstars, config); - - // 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); - } + psastroChipAstrom (chip, config); } Index: /trunk/psastro/src/psastroWCS.c =================================================================== --- /trunk/psastro/src/psastroWCS.c (revision 6176) +++ /trunk/psastro/src/psastroWCS.c (revision 6176) @@ -0,0 +1,158 @@ +# include "psastro.h" + +// 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; + float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2; + float pc1_1, pc1_2, pc2_1, pc2_2; + + // *** interpret header data, convert to crval(i), etc + char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2"); + if (!status) { + psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n"); + return false; + } + + // determine projection type + type = PS_PROJ_NTYPE; + if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN; + if (!strcmp (&ctype[4], "-TAN")) type = PS_PROJ_TAN; + if (!strcmp (&ctype[4], "-AIT")) type = PS_PROJ_AIT; + if (!strcmp (&ctype[4], "-PAR")) type = PS_PROJ_PAR; + if (type == PS_PROJ_NTYPE) { + psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype); + return false; + } + + crval1 = psMetadataLookupF32 (&status, header, "CRVAL1"); + crval2 = psMetadataLookupF32 (&status, header, "CRVAL2"); + crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1"); + crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2"); + + // test the CDELTi varient + cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); + if (status) { + cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2"); + + // test the CROTAi varient: + double rotate = psMetadataLookupF32 (&status, header, "CROTA2"); + if (status) { + double Lambda = cdelt2 / cdelt1; + pc1_1 = cos(rotate*RAD_DEG); + pc1_2 = -sin(rotate*RAD_DEG) * Lambda; + pc2_1 = sin(rotate*RAD_DEG) / Lambda; + pc2_2 = cos(rotate*RAD_DEG); + goto got_matrix; + } + + // test the PC00i00j varient: + pc1_1 = psMetadataLookupF32 (&status, header, "PC001001"); + if (status) { + pc1_2 = psMetadataLookupF32 (&status, header, "PC001002"); + pc2_1 = psMetadataLookupF32 (&status, header, "PC002001"); + pc2_2 = psMetadataLookupF32 (&status, header, "PC002002"); + + // XXX EAM : add Elixir polynomial terms here eventually + goto got_matrix; + } + psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); + return false; + } + + // test the CDi_j varient + pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1"); + if (status) { + pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2"); + pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1"); + pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2"); + + // renormalize to cdelt1, cdelt2, etc + float scale = hypot (pc1_1, pc1_2); + cdelt1 = cdelt2 = scale; + pc1_1 /= scale; + pc1_2 /= scale; + pc2_1 /= scale; + pc2_2 /= scale; + goto got_matrix; + } + psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n"); + return false; + +got_matrix: + + // 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) + // scale = scale(i)/scale(0) (i == chip #) + // project crval1(0),crval2(0 using projection + + // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581) + // XXX EAM : I've fixed this in pslib eam_rel8_b2 + toFPA = psPlaneTransformAlloc (1, 1); + + 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 + 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; +}