IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 13, 2004, 1:58:20 PM (22 years ago)
Author:
gusciora
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/test/astronomy/tst_psCoord.c

    r2051 r2096  
    11/** @file  tst_psImageManip.c
    22*
    3 *  @brief Contains the tests for psImageManip.[ch]
     3*  @brief The code will ...
    44*
    55*
    6 *  @author Robert DeSonia, MHPCC
     6*  @author GLG, MHPCC
    77*
    8 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9 *  @date $Date: 2004-10-12 03:28:47 $
     8*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     9*  @date $Date: 2004-10-13 23:58:20 $
    1010*
    1111*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1818static int test2( void );
    1919static int test3( void );
    20 
    21 // test descrription consists of:
    22 //    * test function
    23 //
    24 //    * testpoint number (item number in gforge testpoint log)
    25 //
    26 //    * name of the testpoint
    27 //
    28 //    * expected return value (negative numbers are signal types, e.g.,
    29 //      -6 for aborts).  Usually either -6 (psAbort called) or 0.
    30 //
    31 //    * boolean to signify if the test function is already in list, i.e., if
    32 //      the test function covers multiple testpoints.  If true, the line is
    33 //      ignored by runTestSuite if the user didn't explicitly specify the
    34 //      testpoint with the -n or -t option, as to not run a test function
    35 //      multiple times.
    36 //
     20static int test4( void );
     21static int test5( void );
     22static int test6( void );
     23static int test7( void );
     24static int test8( void );
     25static int test20( void );
     26static int test21( void );
     27static int test40( void );
     28static int test41( void );
     29
    3730testDescription tests[] = {
    38                               {
    39                                   test1, 1111, "psFunctionBar", 0, false
    40                               },
    41                               {
    42                                   test2, 1113, "psFunctionFoo", 0, false
    43                               },
    44                               {
    45                                   test3, 1114, "psFunctionFunk", 0, false
    46                               },
    47 
    48                               // testpoint #1112 handled by same test function as #1111
    49                               {
    50                                   test1, 1112, "psFunctionBang", 0, true
    51                               },
    52 
    53                               // A null terminates the testDescription list
    54                               {
    55                                   NULL
    56                               }
     31                              {test1, 0000, "psSphereTransformAlloc()", 0, false},
     32                              {test2, 0000, "psPlaneTransformApply()", 0, false},
     33                              {test3, 0000, "psPlaneDistortApply()", 0, false},
     34                              {test4, 0000, "psPSphereTransformApply()", 0, false},
     35                              {test5, 0000, "psSphereTransformICRStoEcliptic()", 0, false},
     36                              {test6, 0000, "psSphereTransformEcliptictoICRS()", 0, false},
     37                              {test7, 0000, "psSphereTransformICRStoGalatic()", 0, false},
     38                              {test8, 0000, "psSphereTransformGalatictoICRS()", 0, false},
     39                              {test20, 0000, "psProject()", 0, false},
     40                              {test21, 0000, "psDeProject()", 0, false},
     41                              {test40, 0000, "psSphereGetOffset()", 0, false},
     42                              {test41, 0000, "psSphereSetOffset()", 0, false},
     43                              {NULL}
    5744                          };
    5845
     46#define ARCSEC_TO_RAD (M_PI / (180.0 * 60.0 * 60.0))
     47#define ARCMIN_TO_RAD (M_PI / (180.0 * 60.0))
     48#define DEG_TO_RAD (M_PI / 180.0)
     49#define RAD_TO_DEG (180.0 / M_PI)
     50#define RAD_TO_ARCMIN ((180.0 * 60.0) / M_PI)
     51#define RAD_TO_ARCSEC ((180.0 * 60.0 * 60.0) / M_PI)
     52
     53
    5954int main( int argc, char* argv[] )
    6055{
    … …  
    6459}
    6560
     61#define MY_TINY 0.0001
     62#define PS_COMPARE_TINY_THEN_PRINT_ERROR(ACTUAL, EXPECT, TESTSTATUS) \
     63if (MY_TINY < fabs(EXPECT - ACTUAL)) { \
     64    printf("ERROR: %s is %f, should be %f\n", #ACTUAL, ACTUAL, EXPECT); \
     65    testStatus = TESTSTATUS; \
     66}
     67#define MY_TINY 0.0001
     68#define PS_COMPARE_EPSILON_THEN_PRINT_ERROR(ACTUAL, EXPECT, TESTSTATUS) \
     69if (MY_TINY < fabs(EXPECT - ACTUAL)) { \
     70    printf("ERROR: %s is %f, should be %f\n", #ACTUAL, ACTUAL, EXPECT); \
     71    testStatus = TESTSTATUS; \
     72}
     73
     74#define PS_COT(X) (1.0 / atan(X))
     75float p_psArg(float x,
     76              float y)
     77{
     78    if (y!= 0.0) {
     79        return(atan(y/x));
     80    } else {
     81        return(0.5 * M_PI);
     82    }
     83    return (0.0);
     84}
     85
     86#define NPLAT 1.0
     87#define X0_BIG 2.0
     88#define X0_SMALL 2.0
     89
    6690int test1( void )
    6791{
    68     // no need to check for memory leaks, as the runTestSuite does that for you.
    69 
    70     // here is where one implements the tests for generally a single testpoint.
    71 
    72     return 0;  // the value that indicates success is part of the testDescription
    73 }
    74 
     92    float tmpF32;
     93    int testStatus = 0;
     94    psSphereTransform *myST = psSphereTransformAlloc(NPLAT, X0_BIG, X0_SMALL);
     95
     96    tmpF32 = sin(NPLAT);
     97    if (FLT_EPSILON < fabs(tmpF32 - myST->sinPhi)) {
     98        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, tmpF32);
     99        testStatus = 1;
     100    }
     101    tmpF32 = cos(NPLAT);
     102    if (FLT_EPSILON < fabs(tmpF32 - myST->cosPhi)) {
     103        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, tmpF32);
     104        testStatus = 2;
     105    }
     106
     107    if (FLT_EPSILON < fabs(X0_BIG - myST->Xo)) {
     108        printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG);
     109        testStatus = 3;
     110    }
     111
     112    if (FLT_EPSILON < fabs(X0_SMALL - myST->xo)) {
     113        printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL);
     114        testStatus = 4;
     115    }
     116
     117    psFree(myST);
     118    return(testStatus);
     119}
     120
     121#define N 10
     122// We do a simple identity transformation on a few x,y pairs.
    75123int test2( void )
    76124{
    77     // no need to check for memory leaks, as the runTestSuite does that for you.
    78 
    79     // here is where one implements the tests for generally a single testpoint.
    80 
    81     return 0;  // the value that indicates success is part of the testDescription
    82 }
    83 
     125    int i;
     126    int testStatus = 0;
     127    psPlane in;
     128    psPlane out;
     129    psPlaneTransform pt;
     130
     131    pt.x = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD);
     132    pt.y = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD);
     133    pt.x->coeff[1][0] = 1.0;
     134    pt.y->coeff[0][1] = 1.0;
     135    for (i=0;i<N;i++) {
     136        in.x = (float) i;
     137        in.y = (float) (i + 5.0);
     138        in.xErr = 0.0;
     139        in.yErr = 0.0;
     140        psPlaneTransformApply(&out, &pt, &in);
     141
     142        if (FLT_EPSILON < fabs(out.x - in.x)) {
     143            printf("ERROR: out.x is %f, should be %f\n", out.x, in.x);
     144            testStatus = 3;
     145        }
     146
     147        if (FLT_EPSILON < fabs(out.y - in.y)) {
     148            printf("ERROR: out.y is %f, should be %f\n", out.y, in.y);
     149            testStatus = 4;
     150        }
     151        //printf("psPlaneTransformApply(): (%f, %f) -> (%f, %f)\n", in.x, in.y, out.x, out.y);
     152    }
     153
     154    psFree(pt.x);
     155    psFree(pt.y);
     156    return(testStatus);
     157}
     158
     159#define COLOR 1.0
     160#define MAGNITUDE 1.0
     161// We do a simple identity transformation on a few x,y pairs.  For x and y,
     162// we add in the COLOR and MAGNITUDE.
    84163int test3( void )
    85164{
    86     // no need to check for memory leaks, as the runTestSuite does that for you.
    87 
    88     // here is where one implements the tests for generally a single testpoint.
    89 
    90     return 0;  // the value that indicates success is part of the testDescription
    91 }
    92 
     165    int i;
     166    int testStatus = 0;
     167    psPlane in;
     168    psPlane out;
     169    psPlaneDistort pt;
     170
     171    pt.x = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD);
     172    pt.y = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD);
     173    pt.x->coeff[1][0][0][0] = 1.0;
     174    pt.x->coeff[0][0][1][0] = 1.0;
     175    pt.x->coeff[0][0][0][1] = 1.0;
     176    pt.y->coeff[0][1][0][0] = 1.0;
     177    pt.y->coeff[0][0][1][0] = 1.0;
     178    pt.y->coeff[0][0][0][1] = 1.0;
     179    for (i=0;i<N;i++) {
     180        in.x = (float) i;
     181        in.y = (float) (i + 5.0);
     182        in.xErr = 0.0;
     183        in.yErr = 0.0;
     184        psPlaneDistortApply(&out, &pt, &in, COLOR, MAGNITUDE);
     185
     186        if (FLT_EPSILON < fabs(out.x - COLOR - MAGNITUDE - in.x)) {
     187            printf("ERROR: out.x is %f, should be %f\n", out.x, in.x + COLOR + MAGNITUDE);
     188            testStatus = 3;
     189        }
     190
     191        if (FLT_EPSILON < fabs(out.y - COLOR - MAGNITUDE - in.y)) {
     192            printf("ERROR: out.y is %f, should be %f\n", out.y, in.y + COLOR + MAGNITUDE);
     193            testStatus = 3;
     194        }
     195        //printf("psPlaneTransformApply(): (%f, %f) -> (%f, %f)\n", in.x, in.y, out.x, out.y);
     196    }
     197
     198    psFree(pt.x);
     199    psFree(pt.y);
     200    return(testStatus);
     201}
     202
     203#define DEG_INC 30.0
     204// We do a simple identity transformation on a few RA, DEC pairs.
     205int test4( void )
     206{
     207    int i;
     208    int testStatus = 0;
     209    psSphere in;
     210    psSphere out;
     211    psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0);
     212
     213    in.r = 45.0 * DEG_TO_RAD;
     214    in.d = 30.0 * DEG_TO_RAD;
     215    in.rErr = 0.0;
     216    in.dErr = 0.0;
     217
     218    for (float r=0.0;r<180.0;r+=DEG_INC) {
     219        for (float d=0.0;d<90.0;d+=DEG_INC) {
     220            in.r = r * DEG_TO_RAD;
     221            in.d = d * DEG_TO_RAD;
     222            in.rErr = 0.0;
     223            in.dErr = 0.0;
     224
     225            psSphereTransformApply(&out, myST, &in);
     226
     227            if (FLT_EPSILON < fabs(out.r - in.r)) {
     228                printf("ERROR: out.r is %f, should be %f\n", out.r, in.r);
     229                testStatus = 4;
     230            }
     231
     232            if (FLT_EPSILON < fabs(out.d - in.d)) {
     233                printf("ERROR: out.d is %f, should be %f\n", out.d, in.d);
     234                testStatus = 5;
     235            }
     236            // printf("psSphereTransformApply (%f, %f) -> (%f, %f)\n", out.r, out.d, in.r, in.d);
     237        }
     238    }
     239
     240    psFree(myST);
     241    return(testStatus);
     242}
     243
     244int test5( void )
     245{
     246    int testStatus = 0;
     247    psTime* now = psTimeGetTime(PS_TIME_UTC);
     248    struct tm *tm_time = psTimeToTM(now);
     249
     250    // XXX: This test code is simply a copy of the original source code.
     251    double year = (double)(1900 + tm_time->tm_year);
     252    double T = year / 100.0;
     253    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
     254    psSphereTransform *myST = psSphereTransformICRStoEcliptic(*now);
     255
     256    if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) {
     257        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi));
     258        testStatus = 1;
     259    }
     260
     261    if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) {
     262        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi));
     263        testStatus = 2;
     264    }
     265
     266    if (FLT_EPSILON < fabs(0.0 - myST->Xo)) {
     267        printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG);
     268        testStatus = 3;
     269    }
     270
     271    if (FLT_EPSILON < fabs(0.0 - myST->xo)) {
     272        printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL);
     273        testStatus = 4;
     274    }
     275
     276    psFree(myST);
     277    psFree(tm_time);
     278    psFree(now);
     279    return(testStatus);
     280}
     281
     282
     283int test6( void )
     284{
     285    int testStatus = 0;
     286    psTime* now = psTimeGetTime(PS_TIME_UTC);
     287    struct tm *tm_time = psTimeToTM(now);
     288
     289    // XXX: This test code is simply a copy of the original source code.
     290    double year = (double)(1900 + tm_time->tm_year);
     291    double T = year / 100.0;
     292    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
     293    psSphereTransform *myST = psSphereTransformEcliptictoICRS(*now);
     294
     295    if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) {
     296        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi));
     297        testStatus = 1;
     298    }
     299
     300    if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) {
     301        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi));
     302        testStatus = 2;
     303    }
     304
     305    if (FLT_EPSILON < fabs(0.0 - myST->Xo)) {
     306        printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG);
     307        testStatus = 3;
     308    }
     309
     310    if (FLT_EPSILON < fabs(0.0 - myST->xo)) {
     311        printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL);
     312        testStatus = 4;
     313    }
     314
     315    psFree(myST);
     316    psFree(tm_time);
     317    psFree(now);
     318    return(testStatus);
     319}
     320
     321
     322int test7( void )
     323{
     324    int testStatus = 0;
     325    // XXX: This test code is simply a copy of the original source code.
     326    double phi = 62.6;
     327    double Xo = 282.25;
     328    double xo = 33.0;
     329    psSphereTransform *myST = psSphereTransformICRStoGalatic();
     330
     331    if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) {
     332        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi));
     333        testStatus = 1;
     334    }
     335
     336    if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) {
     337        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi));
     338        testStatus = 2;
     339    }
     340
     341    if (FLT_EPSILON < fabs(Xo - myST->Xo)) {
     342        printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, Xo);
     343        testStatus = 3;
     344    }
     345
     346    if (FLT_EPSILON < fabs(xo - myST->xo)) {
     347        printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, xo);
     348        testStatus = 4;
     349    }
     350
     351    psFree(myST);
     352    return(testStatus);
     353}
     354
     355int test8( void )
     356{
     357    int testStatus = 0;
     358    // XXX: This test code is simply a copy of the original source code.
     359    double phi = -62.6;
     360    double Xo = 33.0;
     361    double xo = 282.25;
     362    psSphereTransform *myST = psSphereTransformGalatictoICRS();
     363
     364    if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) {
     365        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi));
     366        testStatus = 1;
     367    }
     368
     369    if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) {
     370        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi));
     371        testStatus = 2;
     372    }
     373
     374    if (FLT_EPSILON < fabs(Xo - myST->Xo)) {
     375        printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, Xo);
     376        testStatus = 3;
     377    }
     378
     379    if (FLT_EPSILON < fabs(xo - myST->xo)) {
     380        printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, xo);
     381        testStatus = 4;
     382    }
     383
     384    psFree(myST);
     385    return(testStatus);
     386}
     387
     388int test20( void )
     389{
     390    int i;
     391    int testStatus = 0;
     392    psSphere in;
     393    double R;
     394    double expectX;
     395    double expectY;
     396    psPlane *out;
     397    psProjection myProjection;
     398    myProjection.R = 20.0 * DEG_TO_RAD;
     399    myProjection.D = 10.0 * DEG_TO_RAD;
     400    myProjection.Xs = 1.0;
     401    myProjection.Ys = 1.0;
     402
     403    for (float r=0.0;r<180.0;r+=DEG_INC) {
     404        for (float d=0.0;d<90.0;d+=DEG_INC) {
     405            in.r = r * DEG_TO_RAD;
     406            in.d = d * DEG_TO_RAD;
     407            in.rErr = 0.0;
     408            in.dErr = 0.0;
     409
     410            /******************************************************************
     411             Tangent Plane Projection
     412             *****************************************************************/
     413            myProjection.type = PS_PROJ_TAN;
     414            out = psProject(&in, &myProjection);
     415
     416            R = PS_COT(in.r) * (180.0 / M_PI);
     417            expectX = R * sin(in.d);
     418            expectY = R * cos(in.d);
     419            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 1);
     420            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 2);
     421            psFree(out);
     422
     423            /******************************************************************
     424             PS_PROJ_SIN Projection
     425             *****************************************************************/
     426            myProjection.type = PS_PROJ_SIN;
     427            out = psProject(&in, &myProjection);
     428
     429            R = cos(in.r) * (180.0 / M_PI);
     430            expectX = R * sin(in.d);
     431            expectY = R * cos(in.d);
     432            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 3);
     433            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 4);
     434            psFree(out);
     435
     436            /******************************************************************
     437             PS_PROJ_CAR Projection
     438             *****************************************************************/
     439            myProjection.type = PS_PROJ_CAR;
     440            out = psProject(&in, &myProjection);
     441
     442            expectX = in.d;
     443            expectY = in.r;
     444            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 5);
     445            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 6);
     446            psFree(out);
     447
     448            /******************************************************************
     449             PS_PROJ_MER Projection
     450             *****************************************************************/
     451            myProjection.type = PS_PROJ_MER;
     452            out = psProject(&in, &myProjection);
     453
     454            expectX = in.d;
     455            expectY = log(tan(45.0 + (0.5 * in.r))) * 180.0 / M_PI;
     456
     457            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 7);
     458            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 8);
     459            psFree(out);
     460
     461            /******************************************************************
     462             PS_PROJ_AIT Projection
     463             *****************************************************************/
     464            myProjection.type = PS_PROJ_AIT;
     465            out = psProject(&in, &myProjection);
     466
     467            double alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(in.r) *
     468                                  cos(0.5 * in.d) * 0.5)));
     469            expectX = 2.0 * alpha * cos(in.r) * sin(0.5 * in.d);
     470            expectY = alpha * sin(in.d);
     471
     472            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 7);
     473            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 8);
     474            psFree(out);
     475
     476
     477        }
     478    }
     479    return(testStatus);
     480}
     481
     482#define SPACE_INC 10.0
     483int test21( void )
     484{
     485    int testStatus = 0;
     486    psPlane in;
     487    double R;
     488    double expectD;
     489    double expectR;
     490    psSphere *out;
     491    psProjection myProjection;
     492    myProjection.R = 20.0 * DEG_TO_RAD;
     493    myProjection.D = 10.0 * DEG_TO_RAD;
     494    myProjection.Xs = 1.0;
     495    myProjection.Ys = 1.0;
     496
     497    in.xErr = 0.0;
     498    in.yErr = 0.0;
     499    for (in.x=0.0; in.x<100.0; in.x+=SPACE_INC) {
     500        for (in.y=0.0;in.y<100.0;in.y+=SPACE_INC) {
     501
     502            /******************************************************************
     503             Tangent Plane Projection
     504             *****************************************************************/
     505            myProjection.type = PS_PROJ_TAN;
     506            out = psDeproject(&in, &myProjection);
     507
     508            R = sqrt((in.x * in.x) + (in.y * in.y));
     509            expectD = p_psArg(-in.y, in.x);
     510            expectR = atan(180.0 / (R * M_PI));
     511
     512            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 1);
     513            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 2);
     514            psFree(out);
     515
     516            /******************************************************************
     517             PS_PROJ_SIN Projection
     518             *****************************************************************/
     519            myProjection.type = PS_PROJ_SIN;
     520            out = psDeproject(&in, &myProjection);
     521
     522            R = sqrt((in.x * in.x) + (in.y * in.y));
     523            expectD = p_psArg(-in.y, in.x);
     524            expectR = acos((R * M_PI) / 180.0);
     525
     526            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
     527            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
     528            psFree(out);
     529
     530            /******************************************************************
     531             PS_PROJ_CAR Projection
     532             *****************************************************************/
     533            myProjection.type = PS_PROJ_CAR;
     534            out = psDeproject(&in, &myProjection);
     535
     536            expectD = in.x;
     537            expectR = in.y;
     538
     539            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
     540            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
     541            psFree(out);
     542
     543            /******************************************************************
     544             PS_PROJ_MER Projection
     545             *****************************************************************/
     546            myProjection.type = PS_PROJ_MER;
     547            out = psDeproject(&in, &myProjection);
     548
     549            expectD = in.x;
     550            expectR = (2.0 * atan(exp((in.y * M_PI / 180.0)))) - 180.0;
     551
     552            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
     553            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
     554            psFree(out);
     555
     556
     557            /******************************************************************
     558             PS_PROJ_AIT Projection
     559             *****************************************************************/
     560            myProjection.type = PS_PROJ_AIT;
     561            out = psDeproject(&in, &myProjection);
     562
     563            float chu1 = (in.x * M_PI) / 720.0;
     564            chu1 *= chu1;
     565            float chu2 = (in.y * M_PI) / 360.0;
     566            chu2 *= chu2;
     567            float chu = sqrt(1.0 - chu1 - chu2);
     568            expectD = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (in.x * chu * M_PI) / 360.0);
     569            expectR = asin((in.y * chu * M_PI) / 180.0);
     570
     571            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
     572            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
     573            psFree(out);
     574        }
     575    }
     576    return(testStatus);
     577}
     578
     579int test40( void )
     580{
     581    int testStatus = 0;
     582    psSphere position1;
     583    psSphere position2;
     584    psSphere *offset = NULL;
     585
     586    position1.r = 90.0 * DEG_TO_RAD;
     587    position1.d = 45.0 * DEG_TO_RAD;
     588    position1.rErr = 0.0;
     589    position1.dErr = 0.0;
     590
     591    for (float r=0.0;r<180.0;r+=DEG_INC) {
     592        for (float d=0.0;d<90.0;d+=DEG_INC) {
     593            position2.r = r * DEG_TO_RAD;
     594            position2.d = d * DEG_TO_RAD;
     595            position2.rErr = 0.0;
     596            position2.dErr = 0.0;
     597
     598            offset = psSphereGetOffset(&position1, &position2,
     599                                       PS_SPHERICAL, PS_RADIAN);
     600            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1);
     601            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1);
     602            psFree(offset);
     603
     604            offset = psSphereGetOffset(&position1, &position2,
     605                                       PS_SPHERICAL, PS_DEGREE);
     606            offset->r = offset->r * DEG_TO_RAD;
     607            offset->d = offset->d * DEG_TO_RAD;
     608            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
     609            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
     610            psFree(offset);
     611
     612            offset = psSphereGetOffset(&position1, &position2,
     613                                       PS_SPHERICAL, PS_ARCMIN);
     614            offset->r = offset->r * ARCMIN_TO_RAD;
     615            offset->d = offset->d * ARCMIN_TO_RAD;
     616            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
     617            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
     618            psFree(offset);
     619
     620            offset = psSphereGetOffset(&position1, &position2,
     621                                       PS_SPHERICAL, PS_ARCSEC);
     622            offset->r = offset->r * ARCSEC_TO_RAD;
     623            offset->d = offset->d * ARCSEC_TO_RAD;
     624            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
     625            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
     626            psFree(offset);
     627
     628            /* XXX: This code does not work correctly.
     629                        offset = psSphereGetOffset(&position1, &position2,
     630                                                   PS_LINEAR, 0);
     631                        printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n",
     632                                position1.r, position1.d, position2.r, position2.d,
     633                                offset->r, offset->d);
     634                        PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1);
     635                        PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1);
     636                        psFree(offset);
     637            */
     638        }
     639    }
     640
     641    return(testStatus);
     642}
     643
     644int test41( void )
     645{
     646    int testStatus = 0;
     647    psSphere position1;
     648    psSphere *position2;
     649    psSphere offset;
     650    psSphere tmpOffset;
     651
     652    position1.r = 90.0 * DEG_TO_RAD;
     653    position1.d = 45.0 * DEG_TO_RAD;
     654    position1.rErr = 0.0;
     655    position1.dErr = 0.0;
     656
     657    for (float r=0.0;r<180.0;r+=DEG_INC) {
     658        for (float d=0.0;d<90.0;d+=DEG_INC) {
     659            offset.r = r * DEG_TO_RAD;
     660            offset.d = d * DEG_TO_RAD;
     661            offset.rErr = 0.0;
     662            offset.dErr = 0.0;
     663
     664            position2 = psSphereSetOffset(&position1, &offset,
     665                                          PS_SPHERICAL, PS_RADIAN);
     666            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
     667            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
     668            psFree(position2);
     669
     670            tmpOffset.r = offset.r * RAD_TO_DEG;
     671            tmpOffset.d = offset.d * RAD_TO_DEG;
     672            tmpOffset.rErr = 0.0;
     673            tmpOffset.dErr = 0.0;
     674            position2 = psSphereSetOffset(&position1, &tmpOffset,
     675                                          PS_SPHERICAL, PS_DEGREE);
     676            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
     677            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
     678            psFree(position2);
     679
     680            tmpOffset.r = offset.r * RAD_TO_ARCMIN;
     681            tmpOffset.d = offset.d * RAD_TO_ARCMIN;
     682            tmpOffset.rErr = 0.0;
     683            tmpOffset.dErr = 0.0;
     684            position2 = psSphereSetOffset(&position1, &tmpOffset,
     685                                          PS_SPHERICAL, PS_ARCMIN);
     686            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
     687            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
     688            psFree(position2);
     689
     690            tmpOffset.r = offset.r * RAD_TO_ARCSEC;
     691            tmpOffset.d = offset.d * RAD_TO_ARCSEC;
     692            tmpOffset.rErr = 0.0;
     693            tmpOffset.dErr = 0.0;
     694            position2 = psSphereSetOffset(&position1, &tmpOffset,
     695                                          PS_SPHERICAL, PS_ARCSEC);
     696            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
     697            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
     698            psFree(position2);
     699
     700            /* XXX: This code does not work correctly.
     701                        position2 = psSphereSetOffset(&position1, &offset,
     702                                                       PS_LINEAR, 0);
     703                        printf("--------------- (%f, %f) and (%f, %f) is (%f, %f) ---------------\n",
     704                                position1.r, position1.d, offset.r, offset.d,
     705                                position2->r, position2->d);
     706                        PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
     707                        PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
     708                        psFree(position2);
     709            */
     710        }
     711    }
     712
     713    return(testStatus);
     714}
Note: See TracChangeset for help on using the changeset viewer.