Index: /trunk/psModules/test/astrom/.cvsignore =================================================================== --- /trunk/psModules/test/astrom/.cvsignore (revision 10710) +++ /trunk/psModules/test/astrom/.cvsignore (revision 10711) @@ -8,2 +8,3 @@ tap_pmAstrometryWCS tap_pmAstrometryWCS_DVO2 +tap_pmAstrometryWCS_DVO3 Index: /trunk/psModules/test/astrom/Makefile.am =================================================================== --- /trunk/psModules/test/astrom/Makefile.am (revision 10710) +++ /trunk/psModules/test/astrom/Makefile.am (revision 10711) @@ -16,5 +16,6 @@ tap_pmAstrometryWCS \ tap_pmAstrometryWCS_DVO \ - tap_pmAstrometryWCS_DVO2 + tap_pmAstrometryWCS_DVO2 \ + tap_pmAstrometryWCS_DVO3 test: check Index: /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c =================================================================== --- /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c (revision 10710) +++ /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO.c (revision 10711) @@ -18,4 +18,8 @@ void test3(); // 2nd order term +void test1in(); // basic TAN projection, +void test2in(); // small rotation +void test3in(); // 2nd order term + void testA(); // alloc test void testB(); // alloc test @@ -23,7 +27,12 @@ int main (void) { - plan_tests(492); + plan_tests(984); diag("pmAstrometryWCS tests compared with DVO coords routines"); + diag("this file tests pmAstromWCS <-> Header representations"); + + test1in(); + test2in(); + test3in(); test1(); @@ -33,7 +42,7 @@ } -void test3 () -{ - diag("test pmAstromReadWCS"); +void test3in() +{ + diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -59,42 +68,59 @@ coords.polyterms[2][1] = 0.01; // M vs Y^2 - psMetadata *header = WriteCoordsToHeader (&coords); - pmAstromWCS *wcs = pmAstromWCSfromHeader (header); - - ok (wcs != NULL, "converted WCS keywords to WCS astrometry"); - skip_start (wcs == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); - - psSphere *sky = psSphereAlloc(); - psPlane *chip = psPlaneAlloc(); - - double rDVO, dDVO; + psMetadata *header1 = WriteCoordsToHeader(&coords); + pmAstromWCS *wcs1 = pmAstromWCSfromHeader(header1); + skip_start (wcs1 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + pmAstromWCStoHeader(header2, wcs1); + pmAstromWCS *wcs2 = pmAstromWCSfromHeader(header2); + + psMetadataConfigWrite (header1, "head1.md"); + psMetadataConfigWrite (header2, "head2.md"); + + ok (wcs2 != NULL, "converted WCS keywords to WCS astrometry"); + skip_start (wcs2 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psSphere *sky1 = psSphereAlloc(); + psSphere *sky2 = psSphereAlloc(); + psPlane *chip = psPlaneAlloc(); + for (double x = -2000; x <= +2000; x+= 500.0) { for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); chip->x = x; chip->y = y; - pmAstromWCStoSky (sky, wcs, chip); - while (sky->r > 2*M_PI) - sky->r -= 2*M_PI; - while (sky->r < 0) - sky->r += 2*M_PI; - - ok_float(rDVO, sky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, sky->r*PM_DEG_RAD, rDVO - sky->r*PM_DEG_RAD); - ok_float(dDVO, sky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, sky->d*PM_DEG_RAD, dDVO - sky->d*PM_DEG_RAD); + pmAstromWCStoSky (sky1, wcs1, chip); + pmAstromWCStoSky (sky2, wcs2, chip); + while (sky1->r > 2*M_PI) + sky1->r -= 2*M_PI; + while (sky1->r < 0) + sky1->r += 2*M_PI; + while (sky2->r > 2*M_PI) + sky2->r -= 2*M_PI; + while (sky2->r < 0) + sky2->r += 2*M_PI; + + ok_float(sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, sky1->r*PM_DEG_RAD - sky2->r*PM_DEG_RAD); + ok_float(sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, sky1->d*PM_DEG_RAD - sky2->d*PM_DEG_RAD); } } - psFree (sky); - psFree (chip); - - skip_end(); - psFree (wcs); - psFree (header); - - ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); -} - -void test2 () -{ - diag("test pmAstromReadWCS"); + psFree (sky1); + psFree (sky2); + psFree (chip); + + skip_end(); + psFree (wcs2); + psFree (header2); + + skip_end(); + psFree (wcs1); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test2in() +{ + diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -118,42 +144,59 @@ } - psMetadata *header = WriteCoordsToHeader (&coords); - pmAstromWCS *wcs = pmAstromWCSfromHeader (header); - - ok (wcs != NULL, "converted WCS keywords to WCS astrometry"); - skip_start (wcs == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); - - psSphere *sky = psSphereAlloc(); - psPlane *chip = psPlaneAlloc(); - - double rDVO, dDVO; + psMetadata *header1 = WriteCoordsToHeader(&coords); + pmAstromWCS *wcs1 = pmAstromWCSfromHeader(header1); + skip_start (wcs1 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + pmAstromWCStoHeader(header2, wcs1); + pmAstromWCS *wcs2 = pmAstromWCSfromHeader(header2); + + psMetadataConfigWrite (header1, "head1.md"); + psMetadataConfigWrite (header2, "head2.md"); + + ok (wcs2 != NULL, "converted WCS keywords to WCS astrometry"); + skip_start (wcs2 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psSphere *sky1 = psSphereAlloc(); + psSphere *sky2 = psSphereAlloc(); + psPlane *chip = psPlaneAlloc(); + for (double x = -2000; x <= +2000; x+= 500.0) { for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); chip->x = x; chip->y = y; - pmAstromWCStoSky (sky, wcs, chip); - while (sky->r > 2*M_PI) - sky->r -= 2*M_PI; - while (sky->r < 0) - sky->r += 2*M_PI; - - ok_float(rDVO, sky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, sky->r*PM_DEG_RAD, rDVO - sky->r*PM_DEG_RAD); - ok_float(dDVO, sky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, sky->d*PM_DEG_RAD, dDVO - sky->d*PM_DEG_RAD); + pmAstromWCStoSky (sky1, wcs1, chip); + pmAstromWCStoSky (sky2, wcs2, chip); + while (sky1->r > 2*M_PI) + sky1->r -= 2*M_PI; + while (sky1->r < 0) + sky1->r += 2*M_PI; + while (sky2->r > 2*M_PI) + sky2->r -= 2*M_PI; + while (sky2->r < 0) + sky2->r += 2*M_PI; + + ok_float(sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, sky1->r*PM_DEG_RAD - sky2->r*PM_DEG_RAD); + ok_float(sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, sky1->d*PM_DEG_RAD - sky2->d*PM_DEG_RAD); } } - psFree (sky); - psFree (chip); - - skip_end(); - psFree (wcs); - psFree (header); - - ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); -} - -void test1() -{ - diag("test pmAstromReadWCS"); + psFree (sky1); + psFree (sky2); + psFree (chip); + + skip_end(); + psFree (wcs2); + psFree (header2); + + skip_end(); + psFree (wcs1); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test1in() +{ + diag("test pmAstromWCStoHeader"); psMemId id = psMemGetId(); @@ -177,4 +220,197 @@ } + psMetadata *header1 = WriteCoordsToHeader(&coords); + pmAstromWCS *wcs1 = pmAstromWCSfromHeader(header1); + skip_start (wcs1 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psMetadata *header2 = psMetadataAlloc(); + pmAstromWCStoHeader(header2, wcs1); + pmAstromWCS *wcs2 = pmAstromWCSfromHeader(header2); + + ok (wcs2 != NULL, "converted WCS keywords to WCS astrometry"); + skip_start (wcs2 == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psSphere *sky1 = psSphereAlloc(); + psSphere *sky2 = psSphereAlloc(); + psPlane *chip = psPlaneAlloc(); + + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + chip->x = x; + chip->y = y; + pmAstromWCStoSky (sky1, wcs1, chip); + pmAstromWCStoSky (sky2, wcs2, chip); + while (sky1->r > 2*M_PI) + sky1->r -= 2*M_PI; + while (sky1->r < 0) + sky1->r += 2*M_PI; + while (sky2->r > 2*M_PI) + sky2->r -= 2*M_PI; + while (sky2->r < 0) + sky2->r += 2*M_PI; + + ok_float(sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->r*PM_DEG_RAD, sky2->r*PM_DEG_RAD, sky1->r*PM_DEG_RAD - sky2->r*PM_DEG_RAD); + ok_float(sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", sky1->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD, sky2->d*PM_DEG_RAD - sky2->d*PM_DEG_RAD); + } + } + psFree (sky1); + psFree (sky2); + psFree (chip); + + skip_end(); + psFree (wcs2); + psFree (header2); + + skip_end(); + psFree (wcs1); + psFree (header1); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test3 () +{ + diag("test pmAstromWCSfromHeader "); + 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 *header = WriteCoordsToHeader (&coords); + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + + ok (wcs != NULL, "converted WCS keywords to WCS astrometry"); + skip_start (wcs == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psSphere *sky = psSphereAlloc(); + psPlane *chip = psPlaneAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + chip->x = x; + chip->y = y; + pmAstromWCStoSky (sky, wcs, chip); + while (sky->r > 2*M_PI) + sky->r -= 2*M_PI; + while (sky->r < 0) + sky->r += 2*M_PI; + + ok_float(rDVO, sky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, sky->r*PM_DEG_RAD, rDVO - sky->r*PM_DEG_RAD); + ok_float(dDVO, sky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, sky->d*PM_DEG_RAD, dDVO - sky->d*PM_DEG_RAD); + } + } + psFree (sky); + psFree (chip); + + skip_end(); + psFree (wcs); + psFree (header); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test2 () +{ + diag("test pmAstromWCSfromHeader"); + 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 *header = WriteCoordsToHeader (&coords); + pmAstromWCS *wcs = pmAstromWCSfromHeader (header); + + ok (wcs != NULL, "converted WCS keywords to WCS astrometry"); + skip_start (wcs == NULL, 1, "*** WCS Conversion FAILS *** : skipping related tests"); + + psSphere *sky = psSphereAlloc(); + psPlane *chip = psPlaneAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + chip->x = x; + chip->y = y; + pmAstromWCStoSky (sky, wcs, chip); + while (sky->r > 2*M_PI) + sky->r -= 2*M_PI; + while (sky->r < 0) + sky->r += 2*M_PI; + + ok_float(rDVO, sky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, sky->r*PM_DEG_RAD, rDVO - sky->r*PM_DEG_RAD); + ok_float(dDVO, sky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, sky->d*PM_DEG_RAD, dDVO - sky->d*PM_DEG_RAD); + } + } + psFree (sky); + psFree (chip); + + skip_end(); + psFree (wcs); + psFree (header); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); +} + +void test1() +{ + diag("test pmAstromWCSfromHeader"); + 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 = 0; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + psMetadata *header = WriteCoordsToHeader (&coords); pmAstromWCS *wcs = pmAstromWCSfromHeader (header); Index: /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c =================================================================== --- /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c (revision 10710) +++ /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO2.c (revision 10711) @@ -20,4 +20,5 @@ void test2x(); // small rotation with central offset void test3x(); // 2nd order term with central offset +void test3inv(); // 2nd order term with central offset int main (void) @@ -33,5 +34,365 @@ test2x(); test3x(); + test3inv(); + return exit_status(); +} + +void test1() +{ + 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 = 0; + for (int i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0.0; + coords.polyterms[i][1] = 0.0; + } + + psMetadata *header = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header, 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"); + + psPlane *onChip = psPlaneAlloc(); + psPlane *onFPA = psPlaneAlloc(); + psPlane *onTPA = psPlaneAlloc(); + psSphere *onSky = psSphereAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + + onChip->x = x; + onChip->y = y; + + psPlaneTransformApply (onFPA, chip->toFPA, onChip); + psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); + psDeproject (onSky, onTPA, fpa->toSky); + + while (onSky->r > 2*M_PI) + onSky->r -= 2*M_PI; + while (onSky->r < 0) + onSky->r += 2*M_PI; + + ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); + ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + } + } + psFree (onSky); + psFree (onTPA); + psFree (onFPA); + psFree (onChip); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header); + + 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 *header = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header, 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"); + + psPlane *onChip = psPlaneAlloc(); + psPlane *onFPA = psPlaneAlloc(); + psPlane *onTPA = psPlaneAlloc(); + psSphere *onSky = psSphereAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + + onChip->x = x; + onChip->y = y; + + psPlaneTransformApply (onFPA, chip->toFPA, onChip); + psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); + psDeproject (onSky, onTPA, fpa->toSky); + + while (onSky->r > 2*M_PI) + onSky->r -= 2*M_PI; + while (onSky->r < 0) + onSky->r += 2*M_PI; + + ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); + ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + } + } + psFree (onSky); + psFree (onTPA); + psFree (onFPA); + psFree (onChip); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header); + + 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 *header = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header, 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"); + + psPlane *onChip = psPlaneAlloc(); + psPlane *onFPA = psPlaneAlloc(); + psPlane *onTPA = psPlaneAlloc(); + psSphere *onSky = psSphereAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + + onChip->x = x; + onChip->y = y; + + psPlaneTransformApply (onFPA, chip->toFPA, onChip); + psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); + psDeproject (onSky, onTPA, fpa->toSky); + + while (onSky->r > 2*M_PI) + onSky->r -= 2*M_PI; + while (onSky->r < 0) + onSky->r += 2*M_PI; + + ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); + ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + } + } + psFree (onSky); + psFree (onTPA); + psFree (onFPA); + psFree (onChip); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header); + + 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 = +50.0; + coords.crpix2 = -20.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 *header = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header, 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"); + + psPlane *onChip = psPlaneAlloc(); + psPlane *onFPA = psPlaneAlloc(); + psPlane *onTPA = psPlaneAlloc(); + psSphere *onSky = psSphereAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + + onChip->x = x; + onChip->y = y; + + psPlaneTransformApply (onFPA, chip->toFPA, onChip); + psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); + psDeproject (onSky, onTPA, fpa->toSky); + + while (onSky->r > 2*M_PI) + onSky->r -= 2*M_PI; + while (onSky->r < 0) + onSky->r += 2*M_PI; + + ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); + ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + } + } + psFree (onSky); + psFree (onTPA); + psFree (onFPA); + psFree (onChip); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header); + + 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 = +50.0; + coords.crpix2 = -20.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 *header = WriteCoordsToHeader (&coords); + pmFPA *fpa = pmFPAAlloc (NULL); + pmChip *chip = pmChipAlloc (fpa, NULL); + + bool status = pmAstromReadWCS (fpa, chip, header, 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"); + + psPlane *onChip = psPlaneAlloc(); + psPlane *onFPA = psPlaneAlloc(); + psPlane *onTPA = psPlaneAlloc(); + psSphere *onSky = psSphereAlloc(); + + double rDVO, dDVO; + for (double x = -2000; x <= +2000; x+= 500.0) { + for (double y = -2000; y <= +2000; y+= 500.0) { + XY_to_RD (&rDVO, &dDVO, x, y, &coords); + + onChip->x = x; + onChip->y = y; + + psPlaneTransformApply (onFPA, chip->toFPA, onChip); + psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); + psDeproject (onSky, onTPA, fpa->toSky); + + while (onSky->r > 2*M_PI) + onSky->r -= 2*M_PI; + while (onSky->r < 0) + onSky->r += 2*M_PI; + + ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); + ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + } + } + psFree (onSky); + psFree (onTPA); + psFree (onFPA); + psFree (onChip); + + skip_end(); + psFree (fpa); + psFree (chip); + psFree (header); + + ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); } @@ -110,78 +471,7 @@ } -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 = +50.0; - coords.crpix2 = -20.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 *header = WriteCoordsToHeader (&coords); - pmFPA *fpa = pmFPAAlloc (NULL); - pmChip *chip = pmChipAlloc (fpa, NULL); - - bool status = pmAstromReadWCS (fpa, chip, header, 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"); - - psPlane *onChip = psPlaneAlloc(); - psPlane *onFPA = psPlaneAlloc(); - psPlane *onTPA = psPlaneAlloc(); - psSphere *onSky = psSphereAlloc(); - - double rDVO, dDVO; - for (double x = -2000; x <= +2000; x+= 500.0) { - for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); - - onChip->x = x; - onChip->y = y; - - psPlaneTransformApply (onFPA, chip->toFPA, onChip); - psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); - psDeproject (onSky, onTPA, fpa->toSky); - - while (onSky->r > 2*M_PI) - onSky->r -= 2*M_PI; - while (onSky->r < 0) - onSky->r += 2*M_PI; - - ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); - ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); - } - } - psFree (onSky); - psFree (onTPA); - psFree (onFPA); - psFree (onChip); - - skip_end(); - psFree (fpa); - psFree (chip); - psFree (header); - - ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); -} - -void test1x() -{ - diag("test pmAstromReadWCS"); +void test3inv() +{ + diag("test the inversion of the non-linear polynomial for toFPA -> fromFPA in pmAstromReadWCS"); psMemId id = psMemGetId(); @@ -199,9 +489,12 @@ coords.pc2_1 = 0.0; coords.pc2_2 = 1.0; - coords.Npolyterms = 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 *header = WriteCoordsToHeader (&coords); @@ -213,20 +506,21 @@ skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests"); - psPlane *onChip = psPlaneAlloc(); - psPlane *onFPA = psPlaneAlloc(); - psPlane *onTPA = psPlaneAlloc(); + psPlane *aChip = psPlaneAlloc(); + psPlane *aFPA = psPlaneAlloc(); + psPlane *aTPA = psPlaneAlloc(); + psPlane *bChip = psPlaneAlloc(); + psPlane *bFPA = psPlaneAlloc(); + psPlane *bTPA = psPlaneAlloc(); psSphere *onSky = psSphereAlloc(); - double rDVO, dDVO; for (double x = -2000; x <= +2000; x+= 500.0) { for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); - - onChip->x = x; - onChip->y = y; - - psPlaneTransformApply (onFPA, chip->toFPA, onChip); - psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); - psDeproject (onSky, onTPA, fpa->toSky); + // convert up to sky + aChip->x = x; + aChip->y = y; + + psPlaneTransformApply (aFPA, chip->toFPA, aChip); + psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0); + psDeproject (onSky, aTPA, fpa->toSky); while (onSky->r > 2*M_PI) @@ -235,228 +529,26 @@ onSky->r += 2*M_PI; - ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); - ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); + psProject (bTPA, onSky, fpa->toSky); + psPlaneDistortApply (bFPA, fpa->fromTPA, bTPA, 0.0, 0.0); + 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); + + // 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); } } psFree (onSky); - psFree (onTPA); - psFree (onFPA); - psFree (onChip); - - skip_end(); - psFree (fpa); - psFree (chip); - psFree (header); - - 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 *header = WriteCoordsToHeader (&coords); - pmFPA *fpa = pmFPAAlloc (NULL); - pmChip *chip = pmChipAlloc (fpa, NULL); - - bool status = pmAstromReadWCS (fpa, chip, header, 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"); - - psPlane *onChip = psPlaneAlloc(); - psPlane *onFPA = psPlaneAlloc(); - psPlane *onTPA = psPlaneAlloc(); - psSphere *onSky = psSphereAlloc(); - - double rDVO, dDVO; - for (double x = -2000; x <= +2000; x+= 500.0) { - for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); - - onChip->x = x; - onChip->y = y; - - psPlaneTransformApply (onFPA, chip->toFPA, onChip); - psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); - psDeproject (onSky, onTPA, fpa->toSky); - - while (onSky->r > 2*M_PI) - onSky->r -= 2*M_PI; - while (onSky->r < 0) - onSky->r += 2*M_PI; - - ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); - ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); - } - } - psFree (onSky); - psFree (onTPA); - psFree (onFPA); - psFree (onChip); - - skip_end(); - psFree (fpa); - psFree (chip); - psFree (header); - - 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 *header = WriteCoordsToHeader (&coords); - pmFPA *fpa = pmFPAAlloc (NULL); - pmChip *chip = pmChipAlloc (fpa, NULL); - - bool status = pmAstromReadWCS (fpa, chip, header, 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"); - - psPlane *onChip = psPlaneAlloc(); - psPlane *onFPA = psPlaneAlloc(); - psPlane *onTPA = psPlaneAlloc(); - psSphere *onSky = psSphereAlloc(); - - double rDVO, dDVO; - for (double x = -2000; x <= +2000; x+= 500.0) { - for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); - - onChip->x = x; - onChip->y = y; - - psPlaneTransformApply (onFPA, chip->toFPA, onChip); - psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); - psDeproject (onSky, onTPA, fpa->toSky); - - while (onSky->r > 2*M_PI) - onSky->r -= 2*M_PI; - while (onSky->r < 0) - onSky->r += 2*M_PI; - - ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); - ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); - } - } - psFree (onSky); - psFree (onTPA); - psFree (onFPA); - psFree (onChip); - - skip_end(); - psFree (fpa); - psFree (chip); - psFree (header); - - ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks"); -} - -void test1() -{ - 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 = 0; - for (int i = 0; i < 7; i++) { - coords.polyterms[i][0] = 0.0; - coords.polyterms[i][1] = 0.0; - } - - psMetadata *header = WriteCoordsToHeader (&coords); - pmFPA *fpa = pmFPAAlloc (NULL); - pmChip *chip = pmChipAlloc (fpa, NULL); - - bool status = pmAstromReadWCS (fpa, chip, header, 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"); - - psPlane *onChip = psPlaneAlloc(); - psPlane *onFPA = psPlaneAlloc(); - psPlane *onTPA = psPlaneAlloc(); - psSphere *onSky = psSphereAlloc(); - - double rDVO, dDVO; - for (double x = -2000; x <= +2000; x+= 500.0) { - for (double y = -2000; y <= +2000; y+= 500.0) { - XY_to_RD (&rDVO, &dDVO, x, y, &coords); - - onChip->x = x; - onChip->y = y; - - psPlaneTransformApply (onFPA, chip->toFPA, onChip); - psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0); - psDeproject (onSky, onTPA, fpa->toSky); - - while (onSky->r > 2*M_PI) - onSky->r -= 2*M_PI; - while (onSky->r < 0) - onSky->r += 2*M_PI; - - ok_float(rDVO, onSky->r*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", rDVO, onSky->r*PM_DEG_RAD, rDVO - onSky->r*PM_DEG_RAD); - ok_float(dDVO, onSky->d*PM_DEG_RAD, "coordinate match: %f vs %f (delta = %f)", dDVO, onSky->d*PM_DEG_RAD, dDVO - onSky->d*PM_DEG_RAD); - } - } - psFree (onSky); - psFree (onTPA); - psFree (onFPA); - psFree (onChip); + psFree (aTPA); + psFree (aFPA); + psFree (aChip); + psFree (bTPA); + psFree (bFPA); + psFree (bChip); skip_end(); Index: /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c =================================================================== --- /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c (revision 10711) +++ /trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c (revision 10711) @@ -0,0 +1,190 @@ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include +#include +#include + +#include "tap.h" +#include "pstap.h" + +# if (HAVE_KAPA) + # include "dvo.h" + + psMetadata *WriteCoordsToHeader (Coords *coords); +void test1(); // basic TAN projection, +void test2(); // small rotation +void test3(); // 2nd order term +void test1x(); // basic TAN projection with central offset +void test2x(); // small rotation with central offset +void test3x(); // 2nd order term with central offset +void test3inv(); // 2nd order term with central offset + +int main (void) +{ + plan_tests(984); + + diag("pmAstromWriteWCS tests compared with DVO coords routines"); + + test1(); + // test2(); + // test3(); + // test1x(); + // test2x(); + // test3x(); + // test3inv(); + + return exit_status(); +} + +void test1() +{ + 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 = 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; + + 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"); +} + +psMetadata *WriteCoordsToHeader (Coords *coords) +{ + + char name[16]; + + // construct a header using coords as the input + psMetadata *header = psMetadataAlloc(); + + sprintf (name, "RA--%s", &coords[0].ctype[4]); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", name); + sprintf (name, "DEC-%s", &coords[0].ctype[4]); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", name); + + // center coords (R,D) + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", coords[0].crval1); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", coords[0].crval2); + + // center coords (X,Y) + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", coords[0].crpix1); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", coords[0].crpix2); + + // degrees per pixel + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", coords[0].cdelt1); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", coords[0].cdelt2); + + // rotation matrix + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", coords[0].pc1_1); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", coords[0].pc1_2); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", coords[0].pc2_1); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", coords[0].pc2_2); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", coords[0].polyterms[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", coords[0].polyterms[1][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", coords[0].polyterms[2][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", coords[0].polyterms[0][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", coords[0].polyterms[1][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", coords[0].polyterms[2][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", coords[0].polyterms[3][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", coords[0].polyterms[4][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", coords[0].polyterms[5][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", coords[0].polyterms[6][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", coords[0].polyterms[3][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", coords[0].polyterms[4][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", coords[0].polyterms[5][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", coords[0].polyterms[6][1]); + + psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", coords[0].Npolyterms); + + return header; +} + +# else + + int main (void) +{ + plan_tests(0); + + diag("pmAstrometryWCS tests compared with DVO coords routines : SKIPPED (libdvo not available)"); + + return exit_status(); +} + +# endif +