IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2006, 8:25:29 PM (20 years ago)
Author:
magnier
Message:

more tests of mosaic astrometry

Location:
trunk/psModules/test/astrom
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/test/astrom

    • Property svn:ignore

      •  

        old new  
        88tap_pmAstrometryWCS
        99tap_pmAstrometryWCS_DVO2
         10tap_pmAstrometryWCS_DVO3

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

    r10688 r10711  
    2020void test2x(); // small rotation with central offset
    2121void test3x(); // 2nd order term with central offset
     22void test3inv(); // 2nd order term with central offset
    2223
    2324int main (void)
    … …  
    3334    test2x();
    3435    test3x();
     36    test3inv();
     37
    3538    return exit_status();
     39}
     40
     41void test1()
     42{
     43    diag("test pmAstromReadWCS");
     44    psMemId id = psMemGetId();
     45
     46    // build a DVO-style coordinate system
     47    Coords coords;
     48    strcpy (coords.ctype, "RA---TAN");
     49    coords.crval1 = 0.0;
     50    coords.crval2 = 0.0;
     51    coords.crpix1 = 0.0;
     52    coords.crpix2 = 0.0;
     53    coords.cdelt1 = 1.0/3600;
     54    coords.cdelt2 = 1.0/3600;
     55    coords.pc1_1  = 1.0;
     56    coords.pc1_2  = 0.0;
     57    coords.pc2_1  = 0.0;
     58    coords.pc2_2  = 1.0;
     59    coords.Npolyterms = 0;
     60    for (int i = 0; i < 7; i++) {
     61        coords.polyterms[i][0] = 0.0;
     62        coords.polyterms[i][1] = 0.0;
     63    }
     64
     65    psMetadata *header = WriteCoordsToHeader (&coords);
     66    pmFPA *fpa = pmFPAAlloc (NULL);
     67    pmChip *chip = pmChipAlloc (fpa, NULL);
     68
     69    bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
     70    ok (status, "converted WCS keywords to WCS astrometry");
     71    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     72
     73    psPlane *onChip = psPlaneAlloc();
     74    psPlane *onFPA = psPlaneAlloc();
     75    psPlane *onTPA = psPlaneAlloc();
     76    psSphere *onSky = psSphereAlloc();
     77
     78    double rDVO, dDVO;
     79    for (double x = -2000; x <= +2000; x+= 500.0) {
     80        for (double y = -2000; y <= +2000; y+= 500.0) {
     81            XY_to_RD (&rDVO, &dDVO, x, y, &coords);
     82
     83            onChip->x = x;
     84            onChip->y = y;
     85
     86            psPlaneTransformApply (onFPA, chip->toFPA, onChip);
     87            psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
     88            psDeproject (onSky, onTPA, fpa->toSky);
     89
     90            while (onSky->r > 2*M_PI)
     91                onSky->r -= 2*M_PI;
     92            while (onSky->r <      0)
     93                onSky->r += 2*M_PI;
     94
     95            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);
     96            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);
     97        }
     98    }
     99    psFree (onSky);
     100    psFree (onTPA);
     101    psFree (onFPA);
     102    psFree (onChip);
     103
     104    skip_end();
     105    psFree (fpa);
     106    psFree (chip);
     107    psFree (header);
     108
     109    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     110}
     111
     112void test2()
     113{
     114    diag("test pmAstromReadWCS");
     115    psMemId id = psMemGetId();
     116
     117    // build a DVO-style coordinate system
     118    Coords coords;
     119    strcpy (coords.ctype, "RA---TAN");
     120    coords.crval1 = 0.0;
     121    coords.crval2 = 0.0;
     122    coords.crpix1 = 0.0;
     123    coords.crpix2 = 0.0;
     124    coords.cdelt1 = 1.0/3600;
     125    coords.cdelt2 = 1.0/3600;
     126    coords.pc1_1  = 0.9;
     127    coords.pc1_2  = 0.1;
     128    coords.pc2_1  =-0.1;
     129    coords.pc2_2  = 0.9;
     130    coords.Npolyterms = 0;
     131    for (int i = 0; i < 7; i++) {
     132        coords.polyterms[i][0] = 0.0;
     133        coords.polyterms[i][1] = 0.0;
     134    }
     135
     136    psMetadata *header = WriteCoordsToHeader (&coords);
     137    pmFPA *fpa = pmFPAAlloc (NULL);
     138    pmChip *chip = pmChipAlloc (fpa, NULL);
     139
     140    bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
     141    ok (status, "converted WCS keywords to WCS astrometry");
     142    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     143
     144    psPlane *onChip = psPlaneAlloc();
     145    psPlane *onFPA = psPlaneAlloc();
     146    psPlane *onTPA = psPlaneAlloc();
     147    psSphere *onSky = psSphereAlloc();
     148
     149    double rDVO, dDVO;
     150    for (double x = -2000; x <= +2000; x+= 500.0) {
     151        for (double y = -2000; y <= +2000; y+= 500.0) {
     152            XY_to_RD (&rDVO, &dDVO, x, y, &coords);
     153
     154            onChip->x = x;
     155            onChip->y = y;
     156
     157            psPlaneTransformApply (onFPA, chip->toFPA, onChip);
     158            psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
     159            psDeproject (onSky, onTPA, fpa->toSky);
     160
     161            while (onSky->r > 2*M_PI)
     162                onSky->r -= 2*M_PI;
     163            while (onSky->r <      0)
     164                onSky->r += 2*M_PI;
     165
     166            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);
     167            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);
     168        }
     169    }
     170    psFree (onSky);
     171    psFree (onTPA);
     172    psFree (onFPA);
     173    psFree (onChip);
     174
     175    skip_end();
     176    psFree (fpa);
     177    psFree (chip);
     178    psFree (header);
     179
     180    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     181}
     182
     183void test3()
     184{
     185    diag("test pmAstromReadWCS");
     186    psMemId id = psMemGetId();
     187
     188    // build a DVO-style coordinate system
     189    Coords coords;
     190    strcpy (coords.ctype, "RA---TAN");
     191    coords.crval1 = 0.0;
     192    coords.crval2 = 0.0;
     193    coords.crpix1 = 0.0;
     194    coords.crpix2 = 0.0;
     195    coords.cdelt1 = 1.0/3600;
     196    coords.cdelt2 = 1.0/3600;
     197    coords.pc1_1  = 1.0;
     198    coords.pc1_2  = 0.0;
     199    coords.pc2_1  = 0.0;
     200    coords.pc2_2  = 1.0;
     201    coords.Npolyterms = 2;
     202    for (int i = 0; i < 7; i++) {
     203        coords.polyterms[i][0] = 0.0;
     204        coords.polyterms[i][1] = 0.0;
     205    }
     206    coords.polyterms[0][0] = 0.01; // L vs X^2
     207    coords.polyterms[2][1] = 0.01; // M vs Y^2
     208
     209
     210    psMetadata *header = WriteCoordsToHeader (&coords);
     211    pmFPA *fpa = pmFPAAlloc (NULL);
     212    pmChip *chip = pmChipAlloc (fpa, NULL);
     213
     214    bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
     215    ok (status, "converted WCS keywords to WCS astrometry");
     216    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     217
     218    psPlane *onChip = psPlaneAlloc();
     219    psPlane *onFPA = psPlaneAlloc();
     220    psPlane *onTPA = psPlaneAlloc();
     221    psSphere *onSky = psSphereAlloc();
     222
     223    double rDVO, dDVO;
     224    for (double x = -2000; x <= +2000; x+= 500.0) {
     225        for (double y = -2000; y <= +2000; y+= 500.0) {
     226            XY_to_RD (&rDVO, &dDVO, x, y, &coords);
     227
     228            onChip->x = x;
     229            onChip->y = y;
     230
     231            psPlaneTransformApply (onFPA, chip->toFPA, onChip);
     232            psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
     233            psDeproject (onSky, onTPA, fpa->toSky);
     234
     235            while (onSky->r > 2*M_PI)
     236                onSky->r -= 2*M_PI;
     237            while (onSky->r <      0)
     238                onSky->r += 2*M_PI;
     239
     240            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);
     241            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);
     242        }
     243    }
     244    psFree (onSky);
     245    psFree (onTPA);
     246    psFree (onFPA);
     247    psFree (onChip);
     248
     249    skip_end();
     250    psFree (fpa);
     251    psFree (chip);
     252    psFree (header);
     253
     254    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     255}
     256
     257void test1x()
     258{
     259    diag("test pmAstromReadWCS");
     260    psMemId id = psMemGetId();
     261
     262    // build a DVO-style coordinate system
     263    Coords coords;
     264    strcpy (coords.ctype, "RA---TAN");
     265    coords.crval1 = 0.0;
     266    coords.crval2 = 0.0;
     267    coords.crpix1 = +50.0;
     268    coords.crpix2 = -20.0;
     269    coords.cdelt1 = 1.0/3600;
     270    coords.cdelt2 = 1.0/3600;
     271    coords.pc1_1  = 1.0;
     272    coords.pc1_2  = 0.0;
     273    coords.pc2_1  = 0.0;
     274    coords.pc2_2  = 1.0;
     275    coords.Npolyterms = 0;
     276    for (int i = 0; i < 7; i++) {
     277        coords.polyterms[i][0] = 0.0;
     278        coords.polyterms[i][1] = 0.0;
     279    }
     280
     281    psMetadata *header = WriteCoordsToHeader (&coords);
     282    pmFPA *fpa = pmFPAAlloc (NULL);
     283    pmChip *chip = pmChipAlloc (fpa, NULL);
     284
     285    bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
     286    ok (status, "converted WCS keywords to WCS astrometry");
     287    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     288
     289    psPlane *onChip = psPlaneAlloc();
     290    psPlane *onFPA = psPlaneAlloc();
     291    psPlane *onTPA = psPlaneAlloc();
     292    psSphere *onSky = psSphereAlloc();
     293
     294    double rDVO, dDVO;
     295    for (double x = -2000; x <= +2000; x+= 500.0) {
     296        for (double y = -2000; y <= +2000; y+= 500.0) {
     297            XY_to_RD (&rDVO, &dDVO, x, y, &coords);
     298
     299            onChip->x = x;
     300            onChip->y = y;
     301
     302            psPlaneTransformApply (onFPA, chip->toFPA, onChip);
     303            psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
     304            psDeproject (onSky, onTPA, fpa->toSky);
     305
     306            while (onSky->r > 2*M_PI)
     307                onSky->r -= 2*M_PI;
     308            while (onSky->r <      0)
     309                onSky->r += 2*M_PI;
     310
     311            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);
     312            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);
     313        }
     314    }
     315    psFree (onSky);
     316    psFree (onTPA);
     317    psFree (onFPA);
     318    psFree (onChip);
     319
     320    skip_end();
     321    psFree (fpa);
     322    psFree (chip);
     323    psFree (header);
     324
     325    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
     326}
     327
     328void test2x()
     329{
     330    diag("test pmAstromReadWCS");
     331    psMemId id = psMemGetId();
     332
     333    // build a DVO-style coordinate system
     334    Coords coords;
     335    strcpy (coords.ctype, "RA---TAN");
     336    coords.crval1 = 0.0;
     337    coords.crval2 = 0.0;
     338    coords.crpix1 = +50.0;
     339    coords.crpix2 = -20.0;
     340    coords.cdelt1 = 1.0/3600;
     341    coords.cdelt2 = 1.0/3600;
     342    coords.pc1_1  = 0.9;
     343    coords.pc1_2  = 0.1;
     344    coords.pc2_1  =-0.1;
     345    coords.pc2_2  = 0.9;
     346    coords.Npolyterms = 0;
     347    for (int i = 0; i < 7; i++) {
     348        coords.polyterms[i][0] = 0.0;
     349        coords.polyterms[i][1] = 0.0;
     350    }
     351
     352    psMetadata *header = WriteCoordsToHeader (&coords);
     353    pmFPA *fpa = pmFPAAlloc (NULL);
     354    pmChip *chip = pmChipAlloc (fpa, NULL);
     355
     356    bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
     357    ok (status, "converted WCS keywords to WCS astrometry");
     358    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
     359
     360    psPlane *onChip = psPlaneAlloc();
     361    psPlane *onFPA = psPlaneAlloc();
     362    psPlane *onTPA = psPlaneAlloc();
     363    psSphere *onSky = psSphereAlloc();
     364
     365    double rDVO, dDVO;
     366    for (double x = -2000; x <= +2000; x+= 500.0) {
     367        for (double y = -2000; y <= +2000; y+= 500.0) {
     368            XY_to_RD (&rDVO, &dDVO, x, y, &coords);
     369
     370            onChip->x = x;
     371            onChip->y = y;
     372
     373            psPlaneTransformApply (onFPA, chip->toFPA, onChip);
     374            psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
     375            psDeproject (onSky, onTPA, fpa->toSky);
     376
     377            while (onSky->r > 2*M_PI)
     378                onSky->r -= 2*M_PI;
     379            while (onSky->r <      0)
     380                onSky->r += 2*M_PI;
     381
     382            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);
     383            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);
     384        }
     385    }
     386    psFree (onSky);
     387    psFree (onTPA);
     388    psFree (onFPA);
     389    psFree (onChip);
     390
     391    skip_end();
     392    psFree (fpa);
     393    psFree (chip);
     394    psFree (header);
     395
     396    ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
    36397}
    37398
    … …  
    110471}
    111472
    112 void test2x()
    113 {
    114     diag("test pmAstromReadWCS");
    115     psMemId id = psMemGetId();
    116 
    117     // build a DVO-style coordinate system
    118     Coords coords;
    119     strcpy (coords.ctype, "RA---TAN");
    120     coords.crval1 = 0.0;
    121     coords.crval2 = 0.0;
    122     coords.crpix1 = +50.0;
    123     coords.crpix2 = -20.0;
    124     coords.cdelt1 = 1.0/3600;
    125     coords.cdelt2 = 1.0/3600;
    126     coords.pc1_1  = 0.9;
    127     coords.pc1_2  = 0.1;
    128     coords.pc2_1  =-0.1;
    129     coords.pc2_2  = 0.9;
    130     coords.Npolyterms = 0;
    131     for (int i = 0; i < 7; i++) {
    132         coords.polyterms[i][0] = 0.0;
    133         coords.polyterms[i][1] = 0.0;
    134     }
    135 
    136     psMetadata *header = WriteCoordsToHeader (&coords);
    137     pmFPA *fpa = pmFPAAlloc (NULL);
    138     pmChip *chip = pmChipAlloc (fpa, NULL);
    139 
    140     bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
    141     ok (status, "converted WCS keywords to WCS astrometry");
    142     skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
    143 
    144     psPlane *onChip = psPlaneAlloc();
    145     psPlane *onFPA = psPlaneAlloc();
    146     psPlane *onTPA = psPlaneAlloc();
    147     psSphere *onSky = psSphereAlloc();
    148 
    149     double rDVO, dDVO;
    150     for (double x = -2000; x <= +2000; x+= 500.0) {
    151         for (double y = -2000; y <= +2000; y+= 500.0) {
    152             XY_to_RD (&rDVO, &dDVO, x, y, &coords);
    153 
    154             onChip->x = x;
    155             onChip->y = y;
    156 
    157             psPlaneTransformApply (onFPA, chip->toFPA, onChip);
    158             psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
    159             psDeproject (onSky, onTPA, fpa->toSky);
    160 
    161             while (onSky->r > 2*M_PI)
    162                 onSky->r -= 2*M_PI;
    163             while (onSky->r <      0)
    164                 onSky->r += 2*M_PI;
    165 
    166             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);
    167             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);
    168         }
    169     }
    170     psFree (onSky);
    171     psFree (onTPA);
    172     psFree (onFPA);
    173     psFree (onChip);
    174 
    175     skip_end();
    176     psFree (fpa);
    177     psFree (chip);
    178     psFree (header);
    179 
    180     ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
    181 }
    182 
    183 void test1x()
    184 {
    185     diag("test pmAstromReadWCS");
     473void test3inv()
     474{
     475    diag("test the inversion of the non-linear polynomial for toFPA -> fromFPA in pmAstromReadWCS");
    186476    psMemId id = psMemGetId();
    187477
    … …  
    199489    coords.pc2_1  = 0.0;
    200490    coords.pc2_2  = 1.0;
    201     coords.Npolyterms = 0;
     491    coords.Npolyterms = 2;
    202492    for (int i = 0; i < 7; i++) {
    203493        coords.polyterms[i][0] = 0.0;
    204494        coords.polyterms[i][1] = 0.0;
    205495    }
     496    coords.polyterms[0][0] = 0.01; // L vs X^2
     497    coords.polyterms[2][1] = 0.01; // M vs Y^2
     498
    206499
    207500    psMetadata *header = WriteCoordsToHeader (&coords);
    … …  
    213506    skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
    214507
    215     psPlane *onChip = psPlaneAlloc();
    216     psPlane *onFPA = psPlaneAlloc();
    217     psPlane *onTPA = psPlaneAlloc();
     508    psPlane *aChip  = psPlaneAlloc();
     509    psPlane *aFPA   = psPlaneAlloc();
     510    psPlane *aTPA   = psPlaneAlloc();
     511    psPlane *bChip  = psPlaneAlloc();
     512    psPlane *bFPA   = psPlaneAlloc();
     513    psPlane *bTPA   = psPlaneAlloc();
    218514    psSphere *onSky = psSphereAlloc();
    219515
    220     double rDVO, dDVO;
    221516    for (double x = -2000; x <= +2000; x+= 500.0) {
    222517        for (double y = -2000; y <= +2000; y+= 500.0) {
    223             XY_to_RD (&rDVO, &dDVO, x, y, &coords);
    224 
    225             onChip->x = x;
    226             onChip->y = y;
    227 
    228             psPlaneTransformApply (onFPA, chip->toFPA, onChip);
    229             psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
    230             psDeproject (onSky, onTPA, fpa->toSky);
     518            // convert up to sky
     519            aChip->x = x;
     520            aChip->y = y;
     521
     522            psPlaneTransformApply (aFPA, chip->toFPA, aChip);
     523            psPlaneDistortApply (aTPA, fpa->toTPA, aFPA, 0.0, 0.0);
     524            psDeproject (onSky, aTPA, fpa->toSky);
    231525
    232526            while (onSky->r > 2*M_PI)
    … …  
    235529                onSky->r += 2*M_PI;
    236530
    237             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);
    238             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);
     531            psProject (bTPA, onSky, fpa->toSky);
     532            psPlaneDistortApply (bFPA, fpa->fromTPA, bTPA, 0.0, 0.0);
     533            psPlaneTransformApply (bChip, chip->fromFPA, bFPA);
     534
     535            ok_float(aChip->x, bChip->x, "coordinate match: %f vs %f (delta = %f)", aChip->x, bChip->x, aChip->x - bChip->x);
     536            ok_float(aChip->y, bChip->y, "coordinate match: %f vs %f (delta = %f)", aChip->y, bChip->y, aChip->y - bChip->y);
     537
     538            // ok_float(aFPA->x, bFPA->x, "coordinate match: %f vs %f (delta = %f)", aFPA->x, bFPA->x, aFPA->x - bFPA->x);
     539            // ok_float(aFPA->y, bFPA->y, "coordinate match: %f vs %f (delta = %f)", aFPA->y, bFPA->y, aFPA->y - bFPA->y);
     540
     541            // in this example, TPA coordinates are 10 arcsec/mm; the tol. below represent 1nano-arcsec
     542            // ok_float_tol(aTPA->x, bTPA->x, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->x, bTPA->x, aTPA->x - bTPA->x);
     543            // ok_float_tol(aTPA->y, bTPA->y, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->y, bTPA->y, aTPA->y - bTPA->y);
    239544        }
    240545    }
    241546    psFree (onSky);
    242     psFree (onTPA);
    243     psFree (onFPA);
    244     psFree (onChip);
    245 
    246     skip_end();
    247     psFree (fpa);
    248     psFree (chip);
    249     psFree (header);
    250 
    251     ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
    252 }
    253 
    254 void test3()
    255 {
    256     diag("test pmAstromReadWCS");
    257     psMemId id = psMemGetId();
    258 
    259     // build a DVO-style coordinate system
    260     Coords coords;
    261     strcpy (coords.ctype, "RA---TAN");
    262     coords.crval1 = 0.0;
    263     coords.crval2 = 0.0;
    264     coords.crpix1 = 0.0;
    265     coords.crpix2 = 0.0;
    266     coords.cdelt1 = 1.0/3600;
    267     coords.cdelt2 = 1.0/3600;
    268     coords.pc1_1  = 1.0;
    269     coords.pc1_2  = 0.0;
    270     coords.pc2_1  = 0.0;
    271     coords.pc2_2  = 1.0;
    272     coords.Npolyterms = 2;
    273     for (int i = 0; i < 7; i++) {
    274         coords.polyterms[i][0] = 0.0;
    275         coords.polyterms[i][1] = 0.0;
    276     }
    277     coords.polyterms[0][0] = 0.01; // L vs X^2
    278     coords.polyterms[2][1] = 0.01; // M vs Y^2
    279 
    280 
    281     psMetadata *header = WriteCoordsToHeader (&coords);
    282     pmFPA *fpa = pmFPAAlloc (NULL);
    283     pmChip *chip = pmChipAlloc (fpa, NULL);
    284 
    285     bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
    286     ok (status, "converted WCS keywords to WCS astrometry");
    287     skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
    288 
    289     psPlane *onChip = psPlaneAlloc();
    290     psPlane *onFPA = psPlaneAlloc();
    291     psPlane *onTPA = psPlaneAlloc();
    292     psSphere *onSky = psSphereAlloc();
    293 
    294     double rDVO, dDVO;
    295     for (double x = -2000; x <= +2000; x+= 500.0) {
    296         for (double y = -2000; y <= +2000; y+= 500.0) {
    297             XY_to_RD (&rDVO, &dDVO, x, y, &coords);
    298 
    299             onChip->x = x;
    300             onChip->y = y;
    301 
    302             psPlaneTransformApply (onFPA, chip->toFPA, onChip);
    303             psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
    304             psDeproject (onSky, onTPA, fpa->toSky);
    305 
    306             while (onSky->r > 2*M_PI)
    307                 onSky->r -= 2*M_PI;
    308             while (onSky->r <      0)
    309                 onSky->r += 2*M_PI;
    310 
    311             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);
    312             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);
    313         }
    314     }
    315     psFree (onSky);
    316     psFree (onTPA);
    317     psFree (onFPA);
    318     psFree (onChip);
    319 
    320     skip_end();
    321     psFree (fpa);
    322     psFree (chip);
    323     psFree (header);
    324 
    325     ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
    326 }
    327 
    328 void test2()
    329 {
    330     diag("test pmAstromReadWCS");
    331     psMemId id = psMemGetId();
    332 
    333     // build a DVO-style coordinate system
    334     Coords coords;
    335     strcpy (coords.ctype, "RA---TAN");
    336     coords.crval1 = 0.0;
    337     coords.crval2 = 0.0;
    338     coords.crpix1 = 0.0;
    339     coords.crpix2 = 0.0;
    340     coords.cdelt1 = 1.0/3600;
    341     coords.cdelt2 = 1.0/3600;
    342     coords.pc1_1  = 0.9;
    343     coords.pc1_2  = 0.1;
    344     coords.pc2_1  =-0.1;
    345     coords.pc2_2  = 0.9;
    346     coords.Npolyterms = 0;
    347     for (int i = 0; i < 7; i++) {
    348         coords.polyterms[i][0] = 0.0;
    349         coords.polyterms[i][1] = 0.0;
    350     }
    351 
    352     psMetadata *header = WriteCoordsToHeader (&coords);
    353     pmFPA *fpa = pmFPAAlloc (NULL);
    354     pmChip *chip = pmChipAlloc (fpa, NULL);
    355 
    356     bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
    357     ok (status, "converted WCS keywords to WCS astrometry");
    358     skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
    359 
    360     psPlane *onChip = psPlaneAlloc();
    361     psPlane *onFPA = psPlaneAlloc();
    362     psPlane *onTPA = psPlaneAlloc();
    363     psSphere *onSky = psSphereAlloc();
    364 
    365     double rDVO, dDVO;
    366     for (double x = -2000; x <= +2000; x+= 500.0) {
    367         for (double y = -2000; y <= +2000; y+= 500.0) {
    368             XY_to_RD (&rDVO, &dDVO, x, y, &coords);
    369 
    370             onChip->x = x;
    371             onChip->y = y;
    372 
    373             psPlaneTransformApply (onFPA, chip->toFPA, onChip);
    374             psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
    375             psDeproject (onSky, onTPA, fpa->toSky);
    376 
    377             while (onSky->r > 2*M_PI)
    378                 onSky->r -= 2*M_PI;
    379             while (onSky->r <      0)
    380                 onSky->r += 2*M_PI;
    381 
    382             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);
    383             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);
    384         }
    385     }
    386     psFree (onSky);
    387     psFree (onTPA);
    388     psFree (onFPA);
    389     psFree (onChip);
    390 
    391     skip_end();
    392     psFree (fpa);
    393     psFree (chip);
    394     psFree (header);
    395 
    396     ok(!psMemCheckLeaks (id, NULL, stderr, false), "no memory leaks");
    397 }
    398 
    399 void test1()
    400 {
    401     diag("test pmAstromReadWCS");
    402     psMemId id = psMemGetId();
    403 
    404     // build a DVO-style coordinate system
    405     Coords coords;
    406     strcpy (coords.ctype, "RA---TAN");
    407     coords.crval1 = 0.0;
    408     coords.crval2 = 0.0;
    409     coords.crpix1 = 0.0;
    410     coords.crpix2 = 0.0;
    411     coords.cdelt1 = 1.0/3600;
    412     coords.cdelt2 = 1.0/3600;
    413     coords.pc1_1  = 1.0;
    414     coords.pc1_2  = 0.0;
    415     coords.pc2_1  = 0.0;
    416     coords.pc2_2  = 1.0;
    417     coords.Npolyterms = 0;
    418     for (int i = 0; i < 7; i++) {
    419         coords.polyterms[i][0] = 0.0;
    420         coords.polyterms[i][1] = 0.0;
    421     }
    422 
    423     psMetadata *header = WriteCoordsToHeader (&coords);
    424     pmFPA *fpa = pmFPAAlloc (NULL);
    425     pmChip *chip = pmChipAlloc (fpa, NULL);
    426 
    427     bool status = pmAstromReadWCS (fpa, chip, header, PM_RAD_DEG*10.0/3600.0, false);
    428     ok (status, "converted WCS keywords to WCS astrometry");
    429     skip_start (!status, 1, "*** WCS Conversion FAILS *** : skipping related tests");
    430 
    431     psPlane *onChip = psPlaneAlloc();
    432     psPlane *onFPA = psPlaneAlloc();
    433     psPlane *onTPA = psPlaneAlloc();
    434     psSphere *onSky = psSphereAlloc();
    435 
    436     double rDVO, dDVO;
    437     for (double x = -2000; x <= +2000; x+= 500.0) {
    438         for (double y = -2000; y <= +2000; y+= 500.0) {
    439             XY_to_RD (&rDVO, &dDVO, x, y, &coords);
    440 
    441             onChip->x = x;
    442             onChip->y = y;
    443 
    444             psPlaneTransformApply (onFPA, chip->toFPA, onChip);
    445             psPlaneDistortApply (onTPA, fpa->toTPA, onFPA, 0.0, 0.0);
    446             psDeproject (onSky, onTPA, fpa->toSky);
    447 
    448             while (onSky->r > 2*M_PI)
    449                 onSky->r -= 2*M_PI;
    450             while (onSky->r <      0)
    451                 onSky->r += 2*M_PI;
    452 
    453             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);
    454             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);
    455         }
    456     }
    457     psFree (onSky);
    458     psFree (onTPA);
    459     psFree (onFPA);
    460     psFree (onChip);
     547    psFree (aTPA);
     548    psFree (aFPA);
     549    psFree (aChip);
     550    psFree (bTPA);
     551    psFree (bFPA);
     552    psFree (bChip);
    461553
    462554    skip_end();
Note: See TracChangeset for help on using the changeset viewer.