IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10775


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

polynomial fits

Location:
trunk/psModules
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryWCS.c

    r10712 r10775  
    77 *  @author EAM, IfA
    88 *
    9  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2006-12-14 06:26:15 $
     9 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2006-12-16 05:34:54 $
    1111 *
    1212 *  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    190190            wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
    191191            wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
    192             wcs->trans->y->coeff[1][0] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
     192            wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
    193193            return wcs;
    194194        }
    … …  
    207207                    if (i + j < 2)
    208208                        continue;
    209                     if (i + j > fitOrder)
     209                    if (i + j > fitOrder) {
     210                        wcs->trans->x->mask[i][j] = 1;
     211                        wcs->trans->y->mask[i][j] = 1;
    210212                        continue;
     213                    }
    211214                    sprintf (name, "PCA1X%1dY%1d", i, j);
    212215                    wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
    … …  
    487490{
    488491    // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
    489     if (!psPlaneDistortIsIdentity (fpa->toTPA))
     492    if (!psPlaneDistortIsDiagonal (fpa->toTPA))
    490493        psAbort ("psastro", "invalid TPA transformation");
    491494
    … …  
    505508    // crpix1,2 = X,Y(crval1,2)
    506509    // start with linear solution for Xo,Yo:
    507     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]);
    508     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;
    509     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;
     510    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]);
     511    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;
     512    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;
    510513
    511514    // iterate to actual solution: requires small non-linear terms
    … …  
    582585    double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD;
    583586    double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD;
     587    wcs->cdelt1 = cdelt1;
     588    wcs->cdelt2 = cdelt2;
    584589
    585590    // convert wcs->trans to a matrix which yields L,M in pixels
    … …  
    741746}
    742747
    743 // check that the given psPlaneDistort is the identity
    744 bool psPlaneDistortIsIdentity (psPlaneDistort *distort)
     748// check that the given psPlaneDistort is the identity * (Xs,Ys)
     749bool psPlaneDistortIsDiagonal (psPlaneDistort *distort)
    745750{
    746751
    … …  
    787792            if (i + j > order) {
    788793                // high-order cross terms must be masked (eg, x^3 y^2)
     794                status &= distort->x->mask[i][j][0][0];
     795            } else {
    789796                status &= !distort->x->mask[i][j][0][0];
    790             } else {
    791                 if (i + j != 1) {
     797                if ((i == 1) && (i + j == 1)) {
     798                    // linear, diagonal terms must be 1.0
     799                    status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON);
     800                } else {
    792801                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
    793                     status &= (distort->x->coeff[i][j][0][0] == 0.0);
    794                 } else {
    795                     // linear, diagonal terms must be 1.0
    796                     status &= (distort->x->coeff[i][j][0][0] == 1.0);
     802                    status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON);
    797803                }
    798804            }
    … …  
    805811            if (i + j > order) {
    806812                // high-order cross terms must be masked (eg, x^3 y^2)
     813                status &= distort->y->mask[i][j][0][0];
     814            } else {
    807815                status &= !distort->y->mask[i][j][0][0];
    808             } else {
    809                 if (i + j != 1) {
     816                if ((j == 1) && (i + j == 1)) {
     817                    // linear, diagonal terms must be 1.0
     818                    status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON);
     819                } else {
    810820                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
    811                     status &= (distort->y->coeff[i][j][0][0] == 0.0);
    812                 } else {
    813                     // linear, diagonal terms must be 1.0
    814                     status &= (distort->y->coeff[i][j][0][0] == 1.0);
     821                    status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON);
    815822                }
    816823            }

  • trunk/psModules/src/astrom/pmAstrometryWCS.h

    r10712 r10775  
    77*  @author EAM, IfA
    88*
    9 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2006-12-14 06:26:15 $
     9*  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2006-12-16 05:34:54 $
    1111*
    1212*  Copyright 2006 Institute for Astronomy, University of Hawaii
    … …  
    4242// move to pslib
    4343psPlaneDistort *psPlaneDistortIdentity (int order);
    44 bool psPlaneDistortIsIdentity (psPlaneDistort *distort);
     44bool psPlaneDistortIsDiagonal (psPlaneDistort *distort);
    4545
    4646# define PM_DEG_RAD 57.295779513082322

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

    r10711 r10775  
    2929    plan_tests(984);
    3030
    31     diag("pmAstrometryWCS tests compared with DVO coords routines");
    32     diag("this file tests pmAstromWCS <-> Header representations");
     31    // diag("pmAstrometryWCS tests compared with DVO coords routines");
     32    // diag("this file tests pmAstromWCS <-> Header representations");
    3333
    3434    test1in();
    … …  
    4444void test3in()
    4545{
    46     diag("test pmAstromWCStoHeader");
     46    // diag("test pmAstromWCStoHeader");
    4747    psMemId id = psMemGetId();
    4848
    … …  
    122122void test2in()
    123123{
    124     diag("test pmAstromWCStoHeader");
     124    // diag("test pmAstromWCStoHeader");
    125125    psMemId id = psMemGetId();
    126126
    … …  
    198198void test1in()
    199199{
    200     diag("test pmAstromWCStoHeader");
     200    // diag("test pmAstromWCStoHeader");
    201201    psMemId id = psMemGetId();
    202202
    … …  
    271271void test3 ()
    272272{
    273     diag("test pmAstromWCSfromHeader ");
     273    // diag("test pmAstromWCSfromHeader ");
    274274    psMemId id = psMemGetId();
    275275
    … …  
    332332void test2 ()
    333333{
    334     diag("test pmAstromWCSfromHeader");
     334    // diag("test pmAstromWCSfromHeader");
    335335    psMemId id = psMemGetId();
    336336
    … …  
    391391void test1()
    392392{
    393     diag("test pmAstromWCSfromHeader");
     393    // diag("test pmAstromWCSfromHeader");
    394394    psMemId id = psMemGetId();
    395395
    … …  
    450450void testA()
    451451{
    452     diag("test coord allocs");
     452    // diag("test coord allocs");
    453453    psMemId id = psMemGetId();
    454454
    … …  
    460460void testB()
    461461{
    462     diag("test coord allocs");
     462    // diag("test coord allocs");
    463463    psMemId id = psMemGetId();
    464464
    … …  
    525525    plan_tests(0);
    526526
    527     diag("pmAstrometryWCS tests compared with DVO coords routines : SKIPPED (libdvo not available)");
     527    // diag("pmAstrometryWCS tests compared with DVO coords routines : SKIPPED (libdvo not available)");
    528528
    529529    return exit_status();

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

    r10711 r10775  
    2424int main (void)
    2525{
    26     plan_tests(984);
     26    plan_tests(1472);
    2727
    2828    diag("pmAstromReadWCS tests compared with DVO coords routines");
    … …  
    474474{
    475475    diag("test the inversion of the non-linear polynomial for toFPA -> fromFPA in pmAstromReadWCS");
     476    diag("note that the tolerance for these tests are rather loose");
     477    diag("a 2nd order polynomial is not a great approximate to 1 over a 2nd order polynomial");
     478    diag("unless the non-linear terms are quite small");
    476479    psMemId id = psMemGetId();
    477480
    … …  
    533536            psPlaneTransformApply (bChip, chip->fromFPA, bFPA);
    534537
    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);
     538            // calculate appropriate tol values as f(x,y)
     539            ok_float_tol(aChip->x, bChip->x, 1.0, "coordinate match: %f vs %f (delta = %f)", aChip->x, bChip->x, aChip->x - bChip->x);
     540            ok_float_tol(aChip->y, bChip->y, 1.0, "coordinate match: %f vs %f (delta = %f)", aChip->y, bChip->y, aChip->y - bChip->y);
     541
     542            ok_float(aFPA->x, bFPA->x, "coordinate match: %f vs %f (delta = %f)", aFPA->x, bFPA->x, aFPA->x - bFPA->x);
     543            ok_float(aFPA->y, bFPA->y, "coordinate match: %f vs %f (delta = %f)", aFPA->y, bFPA->y, aFPA->y - bFPA->y);
    540544
    541545            // 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);
     546            ok_float_tol(aTPA->x, bTPA->x, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->x, bTPA->x, aTPA->x - bTPA->x);
     547            ok_float_tol(aTPA->y, bTPA->y, 1e-10, "coordinate match: %f vs %f (delta = %g)", aTPA->y, bTPA->y, aTPA->y - bTPA->y);
    544548        }
    545549    }

  • 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.