Index: trunk/psModules/src/astrom/pmAstrometryWCS.c =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10728) +++ trunk/psModules/src/astrom/pmAstrometryWCS.c (revision 10775) @@ -7,6 +7,6 @@ * @author EAM, IfA * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-14 06:26:15 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-16 05:34:54 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -190,5 +190,5 @@ wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2 wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1 - wcs->trans->y->coeff[1][0] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2 + wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2 return wcs; } @@ -207,6 +207,9 @@ if (i + j < 2) continue; - if (i + j > fitOrder) + if (i + j > fitOrder) { + wcs->trans->x->mask[i][j] = 1; + wcs->trans->y->mask[i][j] = 1; continue; + } sprintf (name, "PCA1X%1dY%1d", i, j); wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name); @@ -487,5 +490,5 @@ { // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1) - if (!psPlaneDistortIsIdentity (fpa->toTPA)) + if (!psPlaneDistortIsDiagonal (fpa->toTPA)) psAbort ("psastro", "invalid TPA transformation"); @@ -505,7 +508,7 @@ // crpix1,2 = X,Y(crval1,2) // start with linear solution for Xo,Yo: - double R = (chip->toFPA->x->coeff[1][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->x->coeff[1][0]); - double Xo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1])/R; - double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R; + double R = (chip->toFPA->x->coeff[1][0]*chip->toFPA->y->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->y->coeff[0][1]); + double Xo = (chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[0][1])/R; + double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[1][0] - chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R; // iterate to actual solution: requires small non-linear terms @@ -582,4 +585,6 @@ double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD; double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD; + wcs->cdelt1 = cdelt1; + wcs->cdelt2 = cdelt2; // convert wcs->trans to a matrix which yields L,M in pixels @@ -741,6 +746,6 @@ } -// check that the given psPlaneDistort is the identity -bool psPlaneDistortIsIdentity (psPlaneDistort *distort) +// check that the given psPlaneDistort is the identity * (Xs,Ys) +bool psPlaneDistortIsDiagonal (psPlaneDistort *distort) { @@ -787,12 +792,13 @@ 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 { status &= !distort->x->mask[i][j][0][0]; - } else { - if (i + j != 1) { + if ((i == 1) && (i + j == 1)) { + // linear, diagonal terms must be 1.0 + status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON); + } else { // 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); + status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON); } } @@ -805,12 +811,13 @@ 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 { status &= !distort->y->mask[i][j][0][0]; - } else { - if (i + j != 1) { + if ((j == 1) && (i + j == 1)) { + // linear, diagonal terms must be 1.0 + status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON); + } else { // 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); + status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON); } } Index: trunk/psModules/src/astrom/pmAstrometryWCS.h =================================================================== --- trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10728) +++ trunk/psModules/src/astrom/pmAstrometryWCS.h (revision 10775) @@ -7,6 +7,6 @@ * @author EAM, IfA * -* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-14 06:26:15 $ +* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-16 05:34:54 $ * * Copyright 2006 Institute for Astronomy, University of Hawaii @@ -42,5 +42,5 @@ // move to pslib psPlaneDistort *psPlaneDistortIdentity (int order); -bool psPlaneDistortIsIdentity (psPlaneDistort *distort); +bool psPlaneDistortIsDiagonal (psPlaneDistort *distort); # define PM_DEG_RAD 57.295779513082322 Index: trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c =================================================================== --- trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c (revision 10728) +++ trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c (revision 10775) @@ -29,6 +29,6 @@ plan_tests(984); - diag("pmAstrometryWCS tests compared with DVO coords routines"); - diag("this file tests pmAstromWCS <-> Header representations"); + // diag("pmAstrometryWCS tests compared with DVO coords routines"); + // diag("this file tests pmAstromWCS <-> Header representations"); test1in(); @@ -44,5 +44,5 @@ void test3in() { - diag("test pmAstromWCStoHeader"); + // diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -122,5 +122,5 @@ void test2in() { - diag("test pmAstromWCStoHeader"); + // diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -198,5 +198,5 @@ void test1in() { - diag("test pmAstromWCStoHeader"); + // diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -271,5 +271,5 @@ void test3 () { - diag("test pmAstromWCSfromHeader "); + // diag("test pmAstromWCSfromHeader "); psMemId id = psMemGetId(); @@ -332,5 +332,5 @@ void test2 () { - diag("test pmAstromWCSfromHeader"); + // diag("test pmAstromWCSfromHeader"); psMemId id = psMemGetId(); @@ -391,5 +391,5 @@ void test1() { - diag("test pmAstromWCSfromHeader"); + // diag("test pmAstromWCSfromHeader"); psMemId id = psMemGetId(); @@ -450,5 +450,5 @@ void testA() { - diag("test coord allocs"); + // diag("test coord allocs"); psMemId id = psMemGetId(); @@ -460,5 +460,5 @@ void testB() { - diag("test coord allocs"); + // diag("test coord allocs"); psMemId id = psMemGetId(); @@ -525,5 +525,5 @@ plan_tests(0); - diag("pmAstrometryWCS tests compared with DVO coords routines : SKIPPED (libdvo not available)"); + // diag("pmAstrometryWCS tests compared with DVO coords routines : SKIPPED (libdvo not available)"); return exit_status(); Index: trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c =================================================================== --- trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c (revision 10728) +++ trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c (revision 10775) @@ -24,5 +24,5 @@ int main (void) { - plan_tests(984); + plan_tests(1472); diag("pmAstromReadWCS tests compared with DVO coords routines"); @@ -474,4 +474,7 @@ { diag("test the inversion of the non-linear polynomial for toFPA -> fromFPA in pmAstromReadWCS"); + diag("note that the tolerance for these tests are rather loose"); + diag("a 2nd order polynomial is not a great approximate to 1 over a 2nd order polynomial"); + diag("unless the non-linear terms are quite small"); psMemId id = psMemGetId(); @@ -533,13 +536,14 @@ psPlaneTransformApply (bChip, chip->fromFPA, bFPA); - ok_float(aChip->x, bChip->x, "coordinate match: %f vs %f (delta = %f)", aChip->x, bChip->x, aChip->x - bChip->x); - ok_float(aChip->y, bChip->y, "coordinate match: %f vs %f (delta = %f)", aChip->y, bChip->y, aChip->y - bChip->y); - - // ok_float(aFPA->x, bFPA->x, "coordinate match: %f vs %f (delta = %f)", aFPA->x, bFPA->x, aFPA->x - bFPA->x); - // ok_float(aFPA->y, bFPA->y, "coordinate match: %f vs %f (delta = %f)", aFPA->y, bFPA->y, aFPA->y - bFPA->y); + // calculate appropriate tol values as f(x,y) + ok_float_tol(aChip->x, bChip->x, 1.0, "coordinate match: %f vs %f (delta = %f)", aChip->x, bChip->x, aChip->x - bChip->x); + ok_float_tol(aChip->y, bChip->y, 1.0, "coordinate match: %f vs %f (delta = %f)", aChip->y, bChip->y, aChip->y - bChip->y); + + ok_float(aFPA->x, bFPA->x, "coordinate match: %f vs %f (delta = %f)", aFPA->x, bFPA->x, aFPA->x - bFPA->x); + ok_float(aFPA->y, bFPA->y, "coordinate match: %f vs %f (delta = %f)", aFPA->y, bFPA->y, aFPA->y - bFPA->y); // in this example, TPA coordinates are 10 arcsec/mm; the tol. below represent 1nano-arcsec - // ok_float_tol(aTPA->x, bTPA->x, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->x, bTPA->x, aTPA->x - bTPA->x); - // ok_float_tol(aTPA->y, bTPA->y, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->y, bTPA->y, aTPA->y - bTPA->y); + ok_float_tol(aTPA->x, bTPA->x, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->x, bTPA->x, aTPA->x - bTPA->x); + ok_float_tol(aTPA->y, bTPA->y, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->y, bTPA->y, aTPA->y - bTPA->y); } } Index: trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c =================================================================== --- trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c (revision 10728) +++ trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c (revision 10775) @@ -28,9 +28,9 @@ diag("pmAstromWriteWCS tests compared with DVO coords routines"); - test1(); + // test1(); // test2(); // test3(); // test1x(); - // test2x(); + test2x(); // test3x(); // test3inv(); @@ -101,4 +101,476 @@ while (aSky->r < 0) aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; + + ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD); + ok_float(aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, aSky->d*PM_DEG_RAD - bSky->d*PM_DEG_RAD); + } + } + psFree (aSky); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + + psFree (bSky); + psFree (bChip); + + psFree (wcs); + psFree (header2); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test2() +{ + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + // build a DVO-style coordinate system + Coords coords; + strcpy (coords.ctype, "RA---TAN"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 0.0; + coords.crpix2 = 0.0; + coords.cdelt1 = 1.0/3600; + coords.cdelt2 = 1.0/3600; + coords.pc1_1 = 0.9; + coords.pc1_2 = 0.1; + coords.pc2_1 = -0.1; + coords.pc2_2 = 0.9; + coords.Npolyterms = 0; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + + psMetadata *header1 = WriteCoordsToHeader (&coords); + psMetadataConfigWrite (header1, "head1.md"); + + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false); + + ok (status, "converted WCS keywords to WCS astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + status = pmAstromWriteWCS (header2, fpa, chip); + psMetadataConfigWrite (header2, "head2.md"); + + pmAstromWCS *wcs = pmAstromWCSfromHeader(header2); + + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psSphere *aSky = psSphereAlloc(); + + psPlane *bChip = psPlaneAlloc(); + psSphere *bSky = psSphereAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + aChip->x = x; + aChip->y = y; + bChip->x = x; + bChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (aSky, aTPA, fpa->toSky); + + pmAstromWCStoSky (bSky, wcs, bChip); + + while (aSky->r > 2*M_PI) + aSky->r -= 2*M_PI; + while (aSky->r < 0) + aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; + + ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD); + ok_float(aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, aSky->d*PM_DEG_RAD - bSky->d*PM_DEG_RAD); + } + } + psFree (aSky); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + + psFree (bSky); + psFree (bChip); + + psFree (wcs); + psFree (header2); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test3() +{ + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + // build a DVO-style coordinate system + Coords coords; + strcpy (coords.ctype, "RA---TAN"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 0.0; + coords.crpix2 = 0.0; + coords.cdelt1 = 1.0/3600; + coords.cdelt2 = 1.0/3600; + coords.pc1_1 = 1.0; + coords.pc1_2 = 0.0; + coords.pc2_1 = 0.0; + coords.pc2_2 = 1.0; + coords.Npolyterms = 2; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + coords.polyterms[0][0] = 0.01; // L vs X^2 + coords.polyterms[2][1] = 0.01; // M vs Y^2 + + psMetadata *header1 = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false); + psMetadataConfigWrite (header1, "head1.md"); + + ok (status, "converted WCS keywords to WCS astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + status = pmAstromWriteWCS (header2, fpa, chip); + psMetadataConfigWrite (header2, "head2.md"); + + pmAstromWCS *wcs = pmAstromWCSfromHeader(header2); + + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psSphere *aSky = psSphereAlloc(); + + psPlane *bChip = psPlaneAlloc(); + psSphere *bSky = psSphereAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + aChip->x = x; + aChip->y = y; + bChip->x = x; + bChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (aSky, aTPA, fpa->toSky); + + pmAstromWCStoSky (bSky, wcs, bChip); + + while (aSky->r > 2*M_PI) + aSky->r -= 2*M_PI; + while (aSky->r < 0) + aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; + + ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD); + ok_float(aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, aSky->d*PM_DEG_RAD - bSky->d*PM_DEG_RAD); + } + } + psFree (aSky); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + + psFree (bSky); + psFree (bChip); + + psFree (wcs); + psFree (header2); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test1x() +{ + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + // build a DVO-style coordinate system + Coords coords; + strcpy (coords.ctype, "RA---TAN"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 20.0; + coords.crpix2 = 50.0; + coords.cdelt1 = 1.0/3600; + coords.cdelt2 = 1.0/3600; + coords.pc1_1 = 1.0; + coords.pc1_2 = 0.0; + coords.pc2_1 = 0.0; + coords.pc2_2 = 1.0; + coords.Npolyterms = 0; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + + psMetadata *header1 = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false); + + ok (status, "converted WCS keywords to WCS astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + status = pmAstromWriteWCS (header2, fpa, chip); + pmAstromWCS *wcs = pmAstromWCSfromHeader(header2); + + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psSphere *aSky = psSphereAlloc(); + + psPlane *bChip = psPlaneAlloc(); + psSphere *bSky = psSphereAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + aChip->x = x; + aChip->y = y; + bChip->x = x; + bChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (aSky, aTPA, fpa->toSky); + + pmAstromWCStoSky (bSky, wcs, bChip); + + while (aSky->r > 2*M_PI) + aSky->r -= 2*M_PI; + while (aSky->r < 0) + aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; + + ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD); + ok_float(aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, aSky->d*PM_DEG_RAD - bSky->d*PM_DEG_RAD); + } + } + psFree (aSky); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + + psFree (bSky); + psFree (bChip); + + psFree (wcs); + psFree (header2); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test2x() +{ + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + // build a DVO-style coordinate system + Coords coords; + strcpy (coords.ctype, "RA---TAN"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 20.0; + coords.crpix2 = 50.0; + coords.cdelt1 = 1.0/3600; + coords.cdelt2 = 1.0/3600; + coords.pc1_1 = 0.9; + coords.pc1_2 = 0.1; + coords.pc2_1 = -0.1; + coords.pc2_2 = 0.9; + coords.Npolyterms = 0; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + + psMetadata *header1 = WriteCoordsToHeader (&coords); + psMetadataConfigWrite (header1, "head1.md"); + + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false); + + ok (status, "converted WCS keywords to WCS astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + status = pmAstromWriteWCS (header2, fpa, chip); + psMetadataConfigWrite (header2, "head2.md"); + + pmAstromWCS *wcs = pmAstromWCSfromHeader(header2); + + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psSphere *aSky = psSphereAlloc(); + + psPlane *bChip = psPlaneAlloc(); + psSphere *bSky = psSphereAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + aChip->x = x; + aChip->y = y; + bChip->x = x; + bChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (aSky, aTPA, fpa->toSky); + + pmAstromWCStoSky (bSky, wcs, bChip); + + while (aSky->r > 2*M_PI) + aSky->r -= 2*M_PI; + while (aSky->r < 0) + aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; + + ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD); + ok_float(aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->d*PM_DEG_RAD, bSky->d*PM_DEG_RAD, aSky->d*PM_DEG_RAD - bSky->d*PM_DEG_RAD); + } + } + psFree (aSky); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + + psFree (bSky); + psFree (bChip); + + psFree (wcs); + psFree (header2); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test3x() +{ + diag("test pmAstromReadWCS"); + psMemId id = psMemGetId(); + + // build a DVO-style coordinate system + Coords coords; + strcpy (coords.ctype, "RA---TAN"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 20.0; + coords.crpix2 = 50.0; + coords.cdelt1 = 1.0/3600; + coords.cdelt2 = 1.0/3600; + coords.pc1_1 = 1.0; + coords.pc1_2 = 0.0; + coords.pc2_1 = 0.0; + coords.pc2_2 = 1.0; + coords.Npolyterms = 2; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + coords.polyterms[0][0] = 0.01; // L vs X^2 + coords.polyterms[2][1] = 0.01; // M vs Y^2 + + psMetadata *header1 = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false); + psMetadataConfigWrite (header1, "head1.md"); + + ok (status, "converted WCS keywords to WCS astrometry"); + skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + status = pmAstromWriteWCS (header2, fpa, chip); + psMetadataConfigWrite (header2, "head2.md"); + + pmAstromWCS *wcs = pmAstromWCSfromHeader(header2); + + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psSphere *aSky = psSphereAlloc(); + + psPlane *bChip = psPlaneAlloc(); + psSphere *bSky = psSphereAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + aChip->x = x; + aChip->y = y; + bChip->x = x; + bChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (aSky, aTPA, fpa->toSky); + + pmAstromWCStoSky (bSky, wcs, bChip); + + while (aSky->r > 2*M_PI) + aSky->r -= 2*M_PI; + while (aSky->r < 0) + aSky->r += 2*M_PI; + while (bSky->r > 2*M_PI) + bSky->r -= 2*M_PI; + while (bSky->r < 0) + bSky->r += 2*M_PI; ok_float(aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", aSky->r*PM_DEG_RAD, bSky->r*PM_DEG_RAD, aSky->r*PM_DEG_RAD - bSky->r*PM_DEG_RAD);