Index: /trunk/psModules/src/astrom/Makefile.am =================================================================== --- /trunk/psModules/src/astrom/Makefile.am (revision 10598) +++ /trunk/psModules/src/astrom/Makefile.am (revision 10599) @@ -5,9 +5,11 @@ libpsmodulesastrom_la_SOURCES = \ pmAstrometryObjects.c \ - pmAstrometryDistortion.c + pmAstrometryDistortion.c \ + pmAstrometryWCS.c pkginclude_HEADERS = \ pmAstrometryObjects.h \ - pmAstrometryDistortion.h + pmAstrometryDistortion.h \ + pmAstrometryWCS.h CLEANFILES = *~ Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10599) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10599) @@ -0,0 +1,641 @@ +/** @file pmAstrometryWCS.c +* +* @brief functions to convert FITS WCS keywords to / from pmFPA structures +* +* @ingroup Astrometry +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-09 21:06:44 $ +* +* Copyright 2006 Institute for Astronomy, University of Hawaii +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "pmFPA.h" +#include "pmAstrometryWCS.h" + +// interpret header WCS (only handles traditional WCS for the moment) +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) +{ + + psProjectionType type; + bool status, pcKeys, cdKeys; + 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 + // XXX add the Elixir DIS / WRP two-layer projection here + 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; + } + + // what type of WCS keywords are available? + psMetadataLookupF32 (&pcKeys, header, "PC001001"); + psMetadataLookupF32 (&cdKeys, header, "CD1_1"); + + if (cdKeys && pcKeys) { + psLogMsg ("psastro", 2, "warning: both CDi_j and PC00i00j defined in headers, using CDi_j terms\n"); + } + if (!cdKeys && !pcKeys) { + psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms"); + // XXX we could default here to RA, DEC, ROTANGLE + 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 + if (pcKeys) { + cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1"); + 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 + if (cdKeys) { + pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1"); + 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 + double 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: + + /***** + + For mosaic astrometry, we need to have a starting set of projection terms in which the + chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or + microns). This set of projections, coupled with an identity toTPA (ie, no distortion) will + result in substantial errors between the observed and predicted star positions on the focal + plane: this is the measurement of the optical distortion in the camera. At the same time, + we need to carry around the transformations which allow us to make an accurate calculation + of the position of the stars based on the input (per-chip) astrometry. These + transformations will allow us to match the raw and ref stars robustly. To convert the + per-chip astrometry (which may have been calculated with a different plate scale for each + chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust + the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for + each chip (eg, cdelt1). Thus, we need to carry around both the + + *****/ + + { + double rX = 1.0; + double rY = 1.0; + + // XXX free an existing toFPA? + psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1); + + // basic transformation from chip to FPA + toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2); + 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] = -(pc2_1*crpix1 + pc2_2*crpix2); + toFPA->y->coeff[1][0] = pc2_1; + toFPA->y->coeff[0][1] = pc2_2; + toFPA->y->mask[1][1] = 1; + + // projection from TPA to SKY + psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); + + if (fpa->toSky == NULL) + { + // XXX for now, use the identity for TPA <--> FPA + fpa->toTangentPlane = psPlaneDistortIdentity (1); + fpa->fromTangentPlane = psPlaneDistortIdentity (1); + fpa->toSky = toSky; + } else + { + if (fpa->toTangentPlane == NULL) + psAbort ("wcs", "projection defined, tangent-plane not defined"); + if (fpa->fromTangentPlane == NULL) + psAbort ("wcs", "projection defined, tangent-plane not defined"); + + // adjust for common toSky for mosaic: + // ignore the TPA since toTPA is identity + // find the FPA coordinate of 0,0 for this chip. + psPlane *fp = psPlaneAlloc(); + psPlane *chip = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + chip->x = chip->y = 0; + + psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate + p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate + p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + + toFPA->x->coeff[0][0] = fp->x; + toFPA->y->coeff[0][0] = fp->y; + + rX = toSky->Xs/fpa->toSky->Xs; + rY = toSky->Ys/fpa->toSky->Ys; + + // correct to common plate scale + toFPA->x->coeff[1][0] *= rX; + toFPA->x->coeff[0][1] *= rX; + toFPA->y->coeff[1][0] *= rY; + toFPA->y->coeff[0][1] *= rY; + + psFree (fp); + psFree (sky); + psFree (chip); + psFree (toSky); + } + + chip->toFPA = toFPA; + chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA); + + while (fpa->toSky->R < 0) + fpa->toSky->R += 2.0*M_PI; + while (fpa->toSky->R > 2.0*M_PI) + fpa->toSky->R -= 2.0*M_PI; + + // remove the correction to the common plate scale + if (isMosaic) + { + chip->toFPA->x->coeff[1][0] /= rX; + chip->toFPA->x->coeff[0][1] /= rX; + chip->toFPA->y->coeff[1][0] /= rY; + chip->toFPA->y->coeff[0][1] /= rY; + } + + psTrace ("psastro", 5, "toFPA: %f %f (%f,%f),(%f,%f)\n", + chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], + chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1], + chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]); + + psTrace ("psastro", 5, "frFPA: %f %f (%f,%f),(%f,%f)\n", + chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0], + chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1], + chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]); + + psLogMsg ("psastro", 3, "field center: %f,%f, plate scale: %f,%f (arcsec/pixel)\n", + DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D, + 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys); + } + return true; +} + +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +// this requires toTP to be the identity transformation +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale) +{ + + // techinically, we can have a plate scale here (toTPA:dx,dy != 1) + if (!psPlaneDistortIsIdentity (toTPA)) + psAbort ("psastro", "invalid TPA transformation"); + + // XXX require toFPA->x->nX == toFPA->x->nY + // XXX require toFPA->y->nX == toFPA->y->nY + // XXX require toFPA->x->nX == toFPA->y->nX + // XXX require toFPA->nX == 1,2,3 + + switch (toSky->type) { + case PS_PROJ_SIN: + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--SIN"); + break; + case PS_PROJ_TAN: + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--TAN"); + break; + case PS_PROJ_AIT: + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--AIT"); + break; + case PS_PROJ_PAR: + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR"); + break; + default: + psLogMsg ("psastro", 2, "warning: unknown projection type %d\n", toSky->type); + return false; + } + + # if (0) + // XXX not really right: needs to deal with non-identity to coeffs + // XXX actually, totally wrong. fix the conversions + // XXX need to handle the plateScale + + /* discussion of the coord transformations: + X,Y: coord on a chip in pixels + L,M: coord on the focal plane (pixels) + P,Q: coord in the tangent plane (microns or mm?) + R,D: coord on the sky + + this function creates WCS terms which convert directly from chip to sky. + this function requires a linear, unrotated toTPA distortion term + toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel) + tpSky->Xs,Ys defines the plate scale (radians / micron) + */ + + // solve for CDELT1,2 (degrees / pixel) + cdelt1 = DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0]; + cdelt2 = DEG_RAD*toSky->Ys*toTPA->y->coeff[0][1][0][0]; + + // L,M = toFPA(X,Y) + // solve for CRPIX1,2 (Xo,Yo) : L,M(Xo,Yo) = 0,0 + + // linear solution for Xo,Yo: + xcoeff = toFPA->x->coeff; + ycoeff = toFPA->y->coeff; + R = (xcoeff[1][0]*ycoeff[0][1] - xcoeff[0][1]*ycoeff[1][0]); + Xo = det*(ycoeff[0][0]*xcoeff[0][1] - xcoeff[0][0]*ycoeff[0][1]); + Yo = det*(xcoeff[0][0]*ycoeff[1][0] - ycoeff[0][0]*xcoeff[1][0]); + + if (toFPA->x->nX > 1) { + + psPolynomial2D *XdX = psPolynomial2D_dX(toFPA->x); + psPolynomial2D *XdY = psPolynomial2D_dY(toFPA->x); + + psPolynomial2D *YdX = psPolynomial2D_dX(toFPA->y); + psPolynomial2D *YdY = psPolynomial2D_dY(toFPA->y); + + psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32); + psVector *Beta = psVectorAlloc (2, PS_DATA_F32); + + // XXX this loop is rather arbitrary in length... + // XXX measure the error and use as criterion + for (int i = 0; i < 10; i++) { + // NOTE: order is: [y][x] + Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo); + Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo); + Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo); + Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo); + + Beta->data.F32[0] = psPolynomial2DEval (toFPA->x, Xo, Yo); + Beta->data.F32[1] = psPolynomial2DEval (toFPA->y, Xo, Yo); + + psMatrixGJSolveF32 (Alpha, Beta); + + Xo += Beta->data.F32[0]; + Yo += Beta->data.F32[1]; + } + } + + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", Xo); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", Yo); + + psPolynomial2D *xWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->x->nX, toFPA->x->nY); + psPolynomial2D *yWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->y->nX, toFPA->y->nY); + + psPolynomial2D *xPx = psPolynomial2DCopy (toFPA->x); + psPolynomial2D *yPx = psPolynomial2DCopy (toFPA->y); + + // skip the zero order terms + // XXX double check that these relationships are correct + for (int i = 0; i < toFPA->x->nX; i++) { + psPolynomial2D *xPy = psPolynomial2DCopy (xPx); + psPolynomial2D *yPy = psPolynomial2DCopy (yPx); + for (int j = 0; j < toFPA->x->nY; j++) { + xWCS->coords[i][j] = psPolynomial2DEval (xPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j); + yWCS->coords[i][j] = psPolynomial2DEval (yPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j); + psPolynomial2D_dY(xPy, xPy); + psPolynomial2D_dY(yPy, yPy); + } + psPolynomial2D_dX(xPx, xPx); + psPolynomial2D_dX(yPx, yPx); + } + + while (coords[0].crval1 < 0) + coords[0].crval1 += 360.0; + while (coords[0].crval1 > 360.0) + coords[0].crval1 -= 360.0; + + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", Xo); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", Yo); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", xWCS->coeff[1][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", xWCS->coeff[0][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", yWCS->coeff[1][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", yWCS->coeff[0][1]); + + // XXX respect the masks + for (int i = 0; i < xWCS->nX; i++) { + for (int j = 0; j < xWCS->nX; j++) { + if (i + j < 2) + continue; + sprintf (name, "PCA1dX%1dY%1d", i, j); + psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", xWCS->coeff[i][j]); + } + } + for (int i = 0; i < yWCS->nX; i++) { + for (int j = 0; j < yWCS->nX; j++) { + if (i + j < 2) + continue; + sprintf (name, "PCA2dX%1dY%1d", i, j); + psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", yWCS->coeff[i][j]); + } + } + + # else + + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale); + + # endif + + // alternative representations use + // CD1_1 = PC001001*CDELT1, etc + // make these representations optional + + return true; +} + +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) +{ + + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP"); + + // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", 1.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", 1.0); + + if (toFPA->x->nX != toFPA->x->nY) + psAbort ("psastro", "mis-matched tangent plane orders (1)"); + if (toFPA->x->nX != toFPA->y->nX) + psAbort ("psastro", "mis-matched tangent plane orders (2)"); + if (toFPA->x->nX != toFPA->y->nY) + psAbort ("psastro", "mis-matched tangent plane orders (3)"); + + switch (toFPA->x->nX) { + case 3: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]); + + case 2: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]); + + case 1: + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale); + break; + + case 0: + psAbort ("psastro", "invalid tangent plane order"); + } + return true; +} + +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) +{ + + psMetadata *header = psMetadataAlloc (); + + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS"); + + // XXX need to handle the plateScale?? + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale); + + if (toTP->x->nX != toTP->x->nY) + psAbort ("psastro", "mis-matched tangent plane orders (1)"); + if (toTP->x->nX != toTP->y->nX) + psAbort ("psastro", "mis-matched tangent plane orders (2)"); + if (toTP->x->nX != toTP->y->nY) + psAbort ("psastro", "mis-matched tangent plane orders (3)"); + + switch (toTP->x->nX) { + case 3: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]); + + case 2: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]); + + case 1: + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]); + break; + + case 0: + psAbort ("psastro", "invalid tangent plane order"); + } + return header; +} + +// construct a psPlaneDistort which is the identify transformation +psPlaneDistort *psPlaneDistortIdentity (int order) +{ + + psPlaneDistort *distort; + + if (order < 1) + psAbort ("psastro", "invalid order"); + if (order > 3) + psAbort ("psastro", "invalid order"); + + // all coeffs and masks initially set to 0 + distort = psPlaneDistortAlloc (order, order, 0, 0); + + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + distort->x->mask [i][j][0][0] = 1; + distort->y->mask [i][j][0][0] = 1; + } + } + } + distort->x->coeff[1][0][0][0] = 1; + distort->y->coeff[0][1][0][0] = 1; + + return distort; +} + +// check that the given psPlaneDistort is the identity +bool psPlaneDistortIsIdentity (psPlaneDistort *distort) +{ + + int order; + bool status; + + // we currently only support up to 3rd order polynomials + if (distort->x->nX < 1) + return false; + if (distort->x->nY < 1) + return false; + if (distort->y->nX < 1) + return false; + if (distort->y->nY < 1) + return false; + + if (distort->x->nX > 3) + return false; + if (distort->x->nY > 3) + return false; + if (distort->y->nX > 3) + return false; + if (distort->y->nY > 3) + return false; + + if (distort->x->nZ > 0) + return false; + if (distort->x->nT > 0) + return false; + if (distort->y->nZ > 0) + return false; + if (distort->y->nT > 0) + return false; + + if (distort->x->nX != distort->x->nY) + return false; + if (distort->y->nX != distort->y->nY) + return false; + + status = true; + order = distort->x->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= !distort->x->mask[i][j][0][0]; + } else { + if (i + j != 1) { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (distort->x->coeff[i][j][0][0] == 0.0); + } else { + // linear, diagonal terms must be 1.0 + status &= (distort->x->coeff[i][j][0][0] == 1.0); + } + } + } + } + + order = distort->y->nX; + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + // high-order cross terms must be masked (eg, x^3 y^2) + status &= !distort->y->mask[i][j][0][0]; + } else { + if (i + j != 1) { + // non-linear and off-diagonal terms must be 0 (eg, x^2, x y) + status &= (distort->y->coeff[i][j][0][0] == 0.0); + } else { + // linear, diagonal terms must be 1.0 + status &= (distort->y->coeff[i][j][0][0] == 1.0); + } + } + } + } + return status; +} Index: /trunk/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10599) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10599) @@ -0,0 +1,27 @@ +/** @file pmAstrometryDistortion.h +* +* @brief functions to convert FITS WCS keywords to / from pmFPA structures +* +* @ingroup Astrometry +* +* @author EAM, IfA +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-09 21:06:44 $ +* +* Copyright 2006 Institute for Astronomy, University of Hawaii +*/ + +#ifndef PM_ASTROMETRY_WCS_H +#define PM_ASTROMETRY_WCS_H + +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic); +bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale); +bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale); +psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale); + +// move to pslib +psPlaneDistort *psPlaneDistortIdentity (int order); +bool psPlaneDistortIsIdentity (psPlaneDistort *distort); + +#endif // PM_ASTROMETRY_WCS_H Index: /trunk/psModules/test/astrom/tap_psAstrometryWCS.c =================================================================== --- /trunk/psModules/test/astrom/tap_psAstrometryWCS.c (revision 10599) +++ /trunk/psModules/test/astrom/tap_psAstrometryWCS.c (revision 10599) @@ -0,0 +1,20 @@ +#include +#include +#include +#include + +#include "tap.h" +#include "pstap.h" + +int main (void) +{ + + pmModelGroupInit (); + + plan_tests(240); + + diag("pmSourcePhotometry tests"); + + +} +