IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 15, 2006, 7:34:54 PM (20 years ago)
Author:
magnier
Message:

polynomial fits

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/test/astrom/tap_pmAstrometryWCS_DVO3.c

    r10711 r10775  
    2828    diag("pmAstromWriteWCS tests compared with DVO coords routines");
    2929
    30     test1();
     30    //    test1();
    3131    //    test2();
    3232    //    test3();
    3333    //    test1x();
    34     //    test2x();
     34    test2x();
    3535    //    test3x();
    3636    //    test3inv();
    … …  
    101101            while (aSky->r <      0)
    102102                aSky->r += 2*M_PI;
     103            while (bSky->r > 2*M_PI)
     104                bSky->r -= 2*M_PI;
     105            while (bSky->r <      0)
     106                bSky->r += 2*M_PI;
     107
     108            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);
     109            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);
     110        }
     111    }
     112    psFree (aSky);
     113    psFree (aTPA);
     114    psFree (aFPA);
     115    psFree (aChip);
     116
     117    psFree (bSky);
     118    psFree (bChip);
     119
     120    psFree (wcs);
     121    psFree (header2);
     122
     123    skip_end();
     124    psFree (fpa);
     125    psFree (chip);
     126    psFree (header1);
     127
     128    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     129}
     130
     131void test2()
     132{
     133    diag("test pmAstromReadWCS");
     134    psMemId id = psMemGetId();
     135
     136    // build a DVO-style coordinate system
     137    Coords coords;
     138    strcpy (coords.ctype, "RA---TAN");
     139    coords.crval1 = 0.0;
     140    coords.crval2 = 0.0;
     141    coords.crpix1 = 0.0;
     142    coords.crpix2 = 0.0;
     143    coords.cdelt1 = 1.0/3600;
     144    coords.cdelt2 = 1.0/3600;
     145    coords.pc1_1  = 0.9;
     146    coords.pc1_2  = 0.1;
     147    coords.pc2_1  = -0.1;
     148    coords.pc2_2  = 0.9;
     149    coords.Npolyterms = 0;
     150    for (int i = 0; i < 7; i++) {
     151        coords.polyterms[i][0] = 0.0;
     152        coords.polyterms[i][1] = 0.0;
     153    }
     154
     155    psMetadata *header1 = WriteCoordsToHeader (&coords);
     156    psMetadataConfigWrite (header1, "head1.md");
     157
     158    pmFPA *fpa = pmFPAAlloc (NULL);
     159    pmChip *chip = pmChipAlloc (fpa, NULL);
     160
     161    bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false);
     162
     163    ok (status, "converted WCS keywords to WCS astrometry");
     164    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     165
     166    psMetadata *header2 = psMetadataAlloc();
     167    status = pmAstromWriteWCS (header2, fpa, chip);
     168    psMetadataConfigWrite (header2, "head2.md");
     169
     170    pmAstromWCS *wcs = pmAstromWCSfromHeader(header2);
     171
     172    psPlane  *aChip = psPlaneAlloc();
     173    psPlane  *aFPA = psPlaneAlloc();
     174    psPlane  *aTPA = psPlaneAlloc();
     175    psSphere *aSky = psSphereAlloc();
     176
     177    psPlane  *bChip = psPlaneAlloc();
     178    psSphere *bSky = psSphereAlloc();
     179
     180    for (double x = -2000; x <= +2000; x+= 500.0) {
     181        for (double y = -2000; y <= +2000; y+= 500.0) {
     182            aChip->x = x;
     183            aChip->y = y;
     184            bChip->x = x;
     185            bChip->y = y;
     186
     187            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     188            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     189            psDeproject (aSky, aTPA, fpa->toSky);
     190
     191            pmAstromWCStoSky (bSky, wcs, bChip);
     192
     193            while (aSky->r > 2*M_PI)
     194                aSky->r -= 2*M_PI;
     195            while (aSky->r <      0)
     196                aSky->r += 2*M_PI;
     197            while (bSky->r > 2*M_PI)
     198                bSky->r -= 2*M_PI;
     199            while (bSky->r <      0)
     200                bSky->r += 2*M_PI;
     201
     202            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);
     203            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);
     204        }
     205    }
     206    psFree (aSky);
     207    psFree (aTPA);
     208    psFree (aFPA);
     209    psFree (aChip);
     210
     211    psFree (bSky);
     212    psFree (bChip);
     213
     214    psFree (wcs);
     215    psFree (header2);
     216
     217    skip_end();
     218    psFree (fpa);
     219    psFree (chip);
     220    psFree (header1);
     221
     222    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     223}
     224
     225void test3()
     226{
     227    diag("test pmAstromReadWCS");
     228    psMemId id = psMemGetId();
     229
     230    // build a DVO-style coordinate system
     231    Coords coords;
     232    strcpy (coords.ctype, "RA---TAN");
     233    coords.crval1 = 0.0;
     234    coords.crval2 = 0.0;
     235    coords.crpix1 = 0.0;
     236    coords.crpix2 = 0.0;
     237    coords.cdelt1 = 1.0/3600;
     238    coords.cdelt2 = 1.0/3600;
     239    coords.pc1_1  = 1.0;
     240    coords.pc1_2  = 0.0;
     241    coords.pc2_1  = 0.0;
     242    coords.pc2_2  = 1.0;
     243    coords.Npolyterms = 2;
     244    for (int i = 0; i < 7; i++) {
     245        coords.polyterms[i][0] = 0.0;
     246        coords.polyterms[i][1] = 0.0;
     247    }
     248    coords.polyterms[0][0] = 0.01; // L vs X^2
     249    coords.polyterms[2][1] = 0.01; // M vs Y^2
     250
     251    psMetadata *header1 = WriteCoordsToHeader (&coords);
     252    pmFPA *fpa = pmFPAAlloc (NULL);
     253    pmChip *chip = pmChipAlloc (fpa, NULL);
     254
     255    bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false);
     256    psMetadataConfigWrite (header1, "head1.md");
     257
     258    ok (status, "converted WCS keywords to WCS astrometry");
     259    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     260
     261    psMetadata *header2 = psMetadataAlloc();
     262    status = pmAstromWriteWCS (header2, fpa, chip);
     263    psMetadataConfigWrite (header2, "head2.md");
     264
     265    pmAstromWCS *wcs = pmAstromWCSfromHeader(header2);
     266
     267    psPlane  *aChip = psPlaneAlloc();
     268    psPlane  *aFPA = psPlaneAlloc();
     269    psPlane  *aTPA = psPlaneAlloc();
     270    psSphere *aSky = psSphereAlloc();
     271
     272    psPlane  *bChip = psPlaneAlloc();
     273    psSphere *bSky = psSphereAlloc();
     274
     275    for (double x = -2000; x <= +2000; x+= 500.0) {
     276        for (double y = -2000; y <= +2000; y+= 500.0) {
     277            aChip->x = x;
     278            aChip->y = y;
     279            bChip->x = x;
     280            bChip->y = y;
     281
     282            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     283            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     284            psDeproject (aSky, aTPA, fpa->toSky);
     285
     286            pmAstromWCStoSky (bSky, wcs, bChip);
     287
     288            while (aSky->r > 2*M_PI)
     289                aSky->r -= 2*M_PI;
     290            while (aSky->r <      0)
     291                aSky->r += 2*M_PI;
     292            while (bSky->r > 2*M_PI)
     293                bSky->r -= 2*M_PI;
     294            while (bSky->r <      0)
     295                bSky->r += 2*M_PI;
     296
     297            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);
     298            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);
     299        }
     300    }
     301    psFree (aSky);
     302    psFree (aTPA);
     303    psFree (aFPA);
     304    psFree (aChip);
     305
     306    psFree (bSky);
     307    psFree (bChip);
     308
     309    psFree (wcs);
     310    psFree (header2);
     311
     312    skip_end();
     313    psFree (fpa);
     314    psFree (chip);
     315    psFree (header1);
     316
     317    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     318}
     319
     320void test1x()
     321{
     322    diag("test pmAstromReadWCS");
     323    psMemId id = psMemGetId();
     324
     325    // build a DVO-style coordinate system
     326    Coords coords;
     327    strcpy (coords.ctype, "RA---TAN");
     328    coords.crval1 = 0.0;
     329    coords.crval2 = 0.0;
     330    coords.crpix1 = 20.0;
     331    coords.crpix2 = 50.0;
     332    coords.cdelt1 = 1.0/3600;
     333    coords.cdelt2 = 1.0/3600;
     334    coords.pc1_1  = 1.0;
     335    coords.pc1_2  = 0.0;
     336    coords.pc2_1  = 0.0;
     337    coords.pc2_2  = 1.0;
     338    coords.Npolyterms = 0;
     339    for (int i = 0; i < 7; i++) {
     340        coords.polyterms[i][0] = 0.0;
     341        coords.polyterms[i][1] = 0.0;
     342    }
     343
     344    psMetadata *header1 = WriteCoordsToHeader (&coords);
     345    pmFPA *fpa = pmFPAAlloc (NULL);
     346    pmChip *chip = pmChipAlloc (fpa, NULL);
     347
     348    bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false);
     349
     350    ok (status, "converted WCS keywords to WCS astrometry");
     351    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     352
     353    psMetadata *header2 = psMetadataAlloc();
     354    status = pmAstromWriteWCS (header2, fpa, chip);
     355    pmAstromWCS *wcs = pmAstromWCSfromHeader(header2);
     356
     357    psPlane  *aChip = psPlaneAlloc();
     358    psPlane  *aFPA = psPlaneAlloc();
     359    psPlane  *aTPA = psPlaneAlloc();
     360    psSphere *aSky = psSphereAlloc();
     361
     362    psPlane  *bChip = psPlaneAlloc();
     363    psSphere *bSky = psSphereAlloc();
     364
     365    for (double x = -2000; x <= +2000; x+= 500.0) {
     366        for (double y = -2000; y <= +2000; y+= 500.0) {
     367            aChip->x = x;
     368            aChip->y = y;
     369            bChip->x = x;
     370            bChip->y = y;
     371
     372            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     373            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     374            psDeproject (aSky, aTPA, fpa->toSky);
     375
     376            pmAstromWCStoSky (bSky, wcs, bChip);
     377
     378            while (aSky->r > 2*M_PI)
     379                aSky->r -= 2*M_PI;
     380            while (aSky->r <      0)
     381                aSky->r += 2*M_PI;
     382            while (bSky->r > 2*M_PI)
     383                bSky->r -= 2*M_PI;
     384            while (bSky->r <      0)
     385                bSky->r += 2*M_PI;
     386
     387            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);
     388            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);
     389        }
     390    }
     391    psFree (aSky);
     392    psFree (aTPA);
     393    psFree (aFPA);
     394    psFree (aChip);
     395
     396    psFree (bSky);
     397    psFree (bChip);
     398
     399    psFree (wcs);
     400    psFree (header2);
     401
     402    skip_end();
     403    psFree (fpa);
     404    psFree (chip);
     405    psFree (header1);
     406
     407    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     408}
     409
     410void test2x()
     411{
     412    diag("test pmAstromReadWCS");
     413    psMemId id = psMemGetId();
     414
     415    // build a DVO-style coordinate system
     416    Coords coords;
     417    strcpy (coords.ctype, "RA---TAN");
     418    coords.crval1 = 0.0;
     419    coords.crval2 = 0.0;
     420    coords.crpix1 = 20.0;
     421    coords.crpix2 = 50.0;
     422    coords.cdelt1 = 1.0/3600;
     423    coords.cdelt2 = 1.0/3600;
     424    coords.pc1_1  = 0.9;
     425    coords.pc1_2  = 0.1;
     426    coords.pc2_1  = -0.1;
     427    coords.pc2_2  = 0.9;
     428    coords.Npolyterms = 0;
     429    for (int i = 0; i < 7; i++) {
     430        coords.polyterms[i][0] = 0.0;
     431        coords.polyterms[i][1] = 0.0;
     432    }
     433
     434    psMetadata *header1 = WriteCoordsToHeader (&coords);
     435    psMetadataConfigWrite (header1, "head1.md");
     436
     437    pmFPA *fpa = pmFPAAlloc (NULL);
     438    pmChip *chip = pmChipAlloc (fpa, NULL);
     439
     440    bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false);
     441
     442    ok (status, "converted WCS keywords to WCS astrometry");
     443    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     444
     445    psMetadata *header2 = psMetadataAlloc();
     446    status = pmAstromWriteWCS (header2, fpa, chip);
     447    psMetadataConfigWrite (header2, "head2.md");
     448
     449    pmAstromWCS *wcs = pmAstromWCSfromHeader(header2);
     450
     451    psPlane  *aChip = psPlaneAlloc();
     452    psPlane  *aFPA = psPlaneAlloc();
     453    psPlane  *aTPA = psPlaneAlloc();
     454    psSphere *aSky = psSphereAlloc();
     455
     456    psPlane  *bChip = psPlaneAlloc();
     457    psSphere *bSky = psSphereAlloc();
     458
     459    for (double x = -2000; x <= +2000; x+= 500.0) {
     460        for (double y = -2000; y <= +2000; y+= 500.0) {
     461            aChip->x = x;
     462            aChip->y = y;
     463            bChip->x = x;
     464            bChip->y = y;
     465
     466            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     467            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     468            psDeproject (aSky, aTPA, fpa->toSky);
     469
     470            pmAstromWCStoSky (bSky, wcs, bChip);
     471
     472            while (aSky->r > 2*M_PI)
     473                aSky->r -= 2*M_PI;
     474            while (aSky->r <      0)
     475                aSky->r += 2*M_PI;
     476            while (bSky->r > 2*M_PI)
     477                bSky->r -= 2*M_PI;
     478            while (bSky->r <      0)
     479                bSky->r += 2*M_PI;
     480
     481            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);
     482            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);
     483        }
     484    }
     485    psFree (aSky);
     486    psFree (aTPA);
     487    psFree (aFPA);
     488    psFree (aChip);
     489
     490    psFree (bSky);
     491    psFree (bChip);
     492
     493    psFree (wcs);
     494    psFree (header2);
     495
     496    skip_end();
     497    psFree (fpa);
     498    psFree (chip);
     499    psFree (header1);
     500
     501    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     502}
     503
     504void test3x()
     505{
     506    diag("test pmAstromReadWCS");
     507    psMemId id = psMemGetId();
     508
     509    // build a DVO-style coordinate system
     510    Coords coords;
     511    strcpy (coords.ctype, "RA---TAN");
     512    coords.crval1 = 0.0;
     513    coords.crval2 = 0.0;
     514    coords.crpix1 = 20.0;
     515    coords.crpix2 = 50.0;
     516    coords.cdelt1 = 1.0/3600;
     517    coords.cdelt2 = 1.0/3600;
     518    coords.pc1_1  = 1.0;
     519    coords.pc1_2  = 0.0;
     520    coords.pc2_1  = 0.0;
     521    coords.pc2_2  = 1.0;
     522    coords.Npolyterms = 2;
     523    for (int i = 0; i < 7; i++) {
     524        coords.polyterms[i][0] = 0.0;
     525        coords.polyterms[i][1] = 0.0;
     526    }
     527    coords.polyterms[0][0] = 0.01; // L vs X^2
     528    coords.polyterms[2][1] = 0.01; // M vs Y^2
     529
     530    psMetadata *header1 = WriteCoordsToHeader (&coords);
     531    pmFPA *fpa = pmFPAAlloc (NULL);
     532    pmChip *chip = pmChipAlloc (fpa, NULL);
     533
     534    bool status = pmAstromReadWCS (fpa, chip, header1, PM_RAD_DEG*10.0/3600.0, false);
     535    psMetadataConfigWrite (header1, "head1.md");
     536
     537    ok (status, "converted WCS keywords to WCS astrometry");
     538    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     539
     540    psMetadata *header2 = psMetadataAlloc();
     541    status = pmAstromWriteWCS (header2, fpa, chip);
     542    psMetadataConfigWrite (header2, "head2.md");
     543
     544    pmAstromWCS *wcs = pmAstromWCSfromHeader(header2);
     545
     546    psPlane  *aChip = psPlaneAlloc();
     547    psPlane  *aFPA = psPlaneAlloc();
     548    psPlane  *aTPA = psPlaneAlloc();
     549    psSphere *aSky = psSphereAlloc();
     550
     551    psPlane  *bChip = psPlaneAlloc();
     552    psSphere *bSky = psSphereAlloc();
     553
     554    for (double x = -2000; x <= +2000; x+= 500.0) {
     555        for (double y = -2000; y <= +2000; y+= 500.0) {
     556            aChip->x = x;
     557            aChip->y = y;
     558            bChip->x = x;
     559            bChip->y = y;
     560
     561            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     562            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     563            psDeproject (aSky, aTPA, fpa->toSky);
     564
     565            pmAstromWCStoSky (bSky, wcs, bChip);
     566
     567            while (aSky->r > 2*M_PI)
     568                aSky->r -= 2*M_PI;
     569            while (aSky->r <      0)
     570                aSky->r += 2*M_PI;
     571            while (bSky->r > 2*M_PI)
     572                bSky->r -= 2*M_PI;
     573            while (bSky->r <      0)
     574                bSky->r += 2*M_PI;
    103575
    104576            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);
Note: See TracChangeset for help on using the changeset viewer.