Index: /trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 28385) +++ /trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 28386) @@ -378,5 +378,5 @@ // XXX make it optional to write out CDi_j terms, or other versions // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity - if (!wcs->wcsCDkeys) { + if (!wcs->wcsCDkeys) { double cdelt1 = wcs->cdelt1; @@ -384,5 +384,5 @@ psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1); psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2); - + // test the PC00i00j varient: psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1 @@ -390,5 +390,5 @@ psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1 psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2 - + // Elixir-style polynomial terms // XXX currently, Elixir/DVO cannot accept mixed orders @@ -398,18 +398,18 @@ if (fitOrder > 1) { for (int i = 0; i <= fitOrder; i++) { - for (int j = 0; j <= fitOrder; j++) { - if (i + j < 2) - continue; - if (i + j > fitOrder) - continue; - sprintf (name, "PCA1X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - sprintf (name, "PCA2X%1dY%1d", i, j); - psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); - } + for (int j = 0; j <= fitOrder; j++) { + if (i + j < 2) + continue; + if (i + j > fitOrder) + continue; + sprintf (name, "PCA1X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + sprintf (name, "PCA2X%1dY%1d", i, j); + psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j)); + } } psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder); } - + // remove any existing 'CDi_j style' wcs keywords if (psMetadataLookup(header, "CD1_1")) { @@ -419,4 +419,28 @@ psMetadataRemoveKey(header, "CD2_2"); } + + // Remove 'CDi_jX' WCS keywords + psString cd11 = psStringCopy("CD1_1 "); + psString cd12 = psStringCopy("CD1_2 "); + psString cd21 = psStringCopy("CD2_1 "); + psString cd22 = psStringCopy("CD2_2 "); + for (char extra = 'A'; extra <= 'Z'; extra++) { + cd11[strlen(cd11)-1] = extra; + if (psMetadataLookup(header, cd11)) { + cd12[strlen(cd12)-1] = extra; + cd21[strlen(cd21)-1] = extra; + cd22[strlen(cd22)-1] = extra; + psMetadataRemoveKey(header, cd11); + psMetadataRemoveKey(header, cd12); + psMetadataRemoveKey(header, cd21); + psMetadataRemoveKey(header, cd22); + } + } + psFree(cd11); + psFree(cd12); + psFree(cd21); + psFree(cd22); + + } else { @@ -427,8 +451,8 @@ if (psMetadataLookup(header, "PC001001")) { - psMetadataRemoveKey(header, "PC001001"); - psMetadataRemoveKey(header, "PC001002"); - psMetadataRemoveKey(header, "PC002001"); - psMetadataRemoveKey(header, "PC002002"); + psMetadataRemoveKey(header, "PC001001"); + psMetadataRemoveKey(header, "PC001002"); + psMetadataRemoveKey(header, "PC002001"); + psMetadataRemoveKey(header, "PC002002"); } } @@ -685,5 +709,5 @@ } - // generate transform with the original orientation (does this rotate about 'center'?) + // generate transform with the original orientation (does this rotate about 'center'?) psPlaneTransform *tpa2 = psPlaneTransformRotate (NULL, tpa1, -1.0*angle); @@ -967,5 +991,5 @@ // outFPA projection must be defined as the goal - + // the output transformations are: // chip -> FPA : standard linear trans with needed rotation, etc @@ -990,23 +1014,23 @@ for (int i = 0; i < nSamples; i++) { - psSphere srcSky; - psPlane *srcChip = psPlaneAlloc(); - psPlane *dstTP = psPlaneAlloc(); + psSphere srcSky; + psPlane *srcChip = psPlaneAlloc(); + psPlane *dstTP = psPlaneAlloc(); srcChip->x = bounds->x0 + (i * deltaX / nSamples); srcChip->y = y; - psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip); - psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP); - psDeproject (&srcSky, &srcTP, inFPA->toSky); - - // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD); - - psProject (dstTP, &srcSky, outFPA->toSky); + psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip); + psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP); + psDeproject (&srcSky, &srcTP, inFPA->toSky); + + // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD); + + psProject (dstTP, &srcSky, outFPA->toSky); srcChip->x -= bounds->x0; srcChip->y -= bounds->y0; - psArrayAdd (src, 100, srcChip); - psArrayAdd (dst, 100, dstTP); + psArrayAdd (src, 100, srcChip); + psArrayAdd (dst, 100, dstTP); psFree(srcChip); // drop our refs to s and d @@ -1021,23 +1045,23 @@ if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) { psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed"); - psFree(src); - psFree(dst); + psFree(src); + psFree(dst); return NULL; } - + # if (0) for (int i = 0; i < src->n; i++) { - - psSphere srcSky, dstSky; - psPlane *srcChip = src->data[i]; - psPlane *dstTP = dst->data[i]; - - psPlaneTransformApply (&srcFP, newToFPA, srcChip); - psDeproject (&srcSky, &srcFP, outFPA->toSky); - psDeproject (&dstSky, dstTP, outFPA->toSky); - - double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0; - double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0; - fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD); + + psSphere srcSky, dstSky; + psPlane *srcChip = src->data[i]; + psPlane *dstTP = dst->data[i]; + + psPlaneTransformApply (&srcFP, newToFPA, srcChip); + psDeproject (&srcSky, &srcFP, outFPA->toSky); + psDeproject (&dstSky, dstTP, outFPA->toSky); + + double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0; + double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0; + fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD); }