Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 20720) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 21025) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-11-13 19:38:16 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-12-17 19:12:34 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -80,5 +80,5 @@ // convert toFPA / toSky components to traditional WCS -// we require the header to have NAXIS1,NAXIS2, the field of the FPA +// we require the header to have NAXIS1,NAXIS2, the field of the FPA // the center of the TPA/Sky projection is 0.5*(NAXIS1,NAXIS2) bool pmAstromReadBilevelMosaic (pmFPA *fpa, const psMetadata *header) @@ -96,10 +96,10 @@ if (!status1 || !status2) { - Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1"); - Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2"); + Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1"); + Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2"); } if (!status1 || !status2) { - psFree (wcs); + psFree (wcs); psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header"); return false; @@ -465,23 +465,23 @@ } - // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter - // XXX this only works if toTPA is at most a linear transformation + // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter + // XXX this only works if toTPA is at most a linear transformation psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY); - for (int i = 0; i <= toFPA->x->nX; i++) { - for (int j = 0; j <= toFPA->x->nY; j++) { - double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j]; - double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j]; - toFPAnew->x->coeff[i][j] = f1 + f2; - - double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j]; - double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j]; - toFPAnew->y->coeff[i][j] = g1 + g2; - } - } - toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0]; - toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0]; - - psFree (toFPA); - toFPA = toFPAnew; + for (int i = 0; i <= toFPA->x->nX; i++) { + for (int j = 0; j <= toFPA->x->nY; j++) { + double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j]; + double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j]; + toFPAnew->x->coeff[i][j] = f1 + f2; + + double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j]; + double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j]; + toFPAnew->y->coeff[i][j] = g1 + g2; + } + } + toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0]; + toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0]; + + psFree (toFPA); + toFPA = toFPAnew; // adjust reference pixel for new toSky reference coordinate @@ -492,11 +492,11 @@ psSphere *sky = psSphereAlloc(); - sky->r = toSky->R; - sky->d = toSky->D; + sky->r = toSky->R; + sky->d = toSky->D; psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection psPlaneTransformApply (fpOld, fpa->fromTPA, tp); - sky->r = fpa->toSky->R; - sky->d = fpa->toSky->D; + sky->r = fpa->toSky->R; + sky->d = fpa->toSky->D; psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection psPlaneTransformApply (fpNew, fpa->fromTPA, tp); @@ -595,4 +595,9 @@ pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, double tol) { + PS_ASSERT_PTR_NON_NULL(fpa, NULL); + PS_ASSERT_PTR_NON_NULL(chip, NULL); + PS_ASSERT_PTR_NON_NULL(fpa->toFPA, NULL); + PS_ASSERT_PTR_NON_NULL(fpa->toTPA, NULL); + // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY @@ -607,5 +612,5 @@ // create a temporary transform which combines toTPA and toFPA. allow toTPA to have 0th // and 1st order terms - + // XXX require fpa->toTPA->x->nX == fpa->toTPA->x->nY // XXX require fpa->toTPA->y->nX == fpa->toTPA->y->nY @@ -618,13 +623,13 @@ for (int i = 0; i <= toTPA->x->nX; i++) { - for (int j = 0; j <= toTPA->x->nY; j++) { - double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j]; - double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j]; - toTPA->x->coeff[i][j] = f1 + f2; - - double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j]; - double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j]; - toTPA->y->coeff[i][j] = g1 + g2; - } + for (int j = 0; j <= toTPA->x->nY; j++) { + double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j]; + double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j]; + toTPA->x->coeff[i][j] = f1 + f2; + + double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j]; + double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j]; + toTPA->y->coeff[i][j] = g1 + g2; + } } toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0]; @@ -641,8 +646,8 @@ psPlane *center = psPlaneTransformGetCenter (toTPA, tol); if (!center) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); - psFree (toTPA); - psFree (wcs); - return NULL; + psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); + psFree (toTPA); + psFree (wcs); + return NULL; } @@ -703,7 +708,7 @@ psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol); if (!center) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); - psFree (wcs); - return NULL; + psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); + psFree (wcs); + return NULL; } @@ -744,7 +749,7 @@ psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol); if (!center) { - psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); - psFree (wcs); - return NULL; + psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center."); + psFree (wcs); + return NULL; } @@ -801,5 +806,5 @@ /***** - + 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 @@ -814,5 +819,5 @@ 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 - + *****/ @@ -821,6 +826,6 @@ L,M: coord on the focal plane (pixels) P,Q: coord in the tangent plane (microns or mm?) - R,D: coord on the sky - + 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 @@ -830,12 +835,12 @@ /* at this point, we have extracted from the header the WCS terms in the form of a polynomial, - * wcs->trans, which will convert X,Y in pixels to L,M in degrees. we also have the following + * wcs->trans, which will convert X,Y in pixels to L,M in degrees. we also have the following * elements defined: * type (projection type) * crval1,2 (in RA,DEC degrees) - * crpix1,2 + * crpix1,2 * cdelt1,2 (in degrees / pixel) * pixelScale (microns / pixel) - * + * * now we convert wcs->trans to toFPA, which is different from wcs->trans in 3 important ways: * 1) the output is in microns (not degrees): divide by cdelt1,2