IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9976


Ignore:
Timestamp:
Nov 14, 2006, 11:15:43 AM (20 years ago)
Author:
drobbin
Message:

Implemented testing for psEarthOrientation. Done with motion, working on corrections. Edited data file dir in Makefile.am and changed a couple of, if (MJD==NAN) to if (isnan(MJD))

Location:
trunk/psLib
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/Makefile.am

    r8914 r9976  
    11#Makefile for astronomy functions of psLib
    22#
    3 AM_CFLAGS=$(CFLAGS) -DPS_CONFIG_FILE_DEFAULT=\"$(sysconfdir)/pslib/pslib.config\"
     3AM_CFLAGS=$(CFLAGS) -DPS_CONFIG_FILE_DEFAULT=\"$(sysconfdir)/pslib.config\"
    44
    55noinst_LTLIBRARIES = libpslibastro.la
    66
    7 libpslibastro_la_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS) 
     7libpslibastro_la_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS)
    88libpslibastro_la_SOURCES = \
    99        psTime.c \

  • trunk/psLib/src/astro/psEarthOrientation.c

    r9665 r9976  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2006-10-20 00:17:21 $
     10 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2006-11-14 21:15:43 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    284284        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    285285                "Aberration speed should not be equal to 0.\n");
     286        return NULL;
     287    }
     288    if (fabs(speed) > 1.0) {
     289        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     290                "Aberration speed should not greater than the speed of light!.\n");
    286291        return NULL;
    287292    }
    … …  
    445450    // Convert psTime to MJD
    446451    double MJD = psTimeToMJD(time);
    447     if (MJD == NAN) {
     452    if (isnan(MJD)) {
    448453        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    449454                "Time conversion to MJD failed.  Invalid input time.\n");
    … …  
    594599    //Convert the input time to MJD.  If NAN is returned, return NULL for the function.
    595600    double MJD = psTimeToMJD(time);
    596     if (MJD == NAN) {
     601    if (isnan(MJD)) {
    597602        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    598603                "Time conversion to MJD failed.  Invalid input time.\n");
    … …  
    762767    //Check if tidal corrections should be included.
    763768    //If so, make sure values are positive and in the correct range.
    764     if (tidalCorr != NULL && tidalCorr->s != 0.0) {
     769    if (tidalCorr != NULL && fabs(tidalCorr->s) > FLT_EPSILON) {
    765770        int nsec = in->nsec + (int)(tidalCorr->s * 1e9);
    766771        if (nsec > 1e9) {
    … …  
    777782    //Calculate the Julian Date from the input time in UT1 format.
    778783    double T = psTimeToJD(in);
    779     if (T == NAN) {
     784    if (isnan(T)) {
    780785        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    781786                "Time conversion to JD failed.  Invalid input time.\n");
    … …  
    803808
    804809    double MJD = psTimeToMJD(time);
    805     if (MJD == NAN) {
     810    if ( isnan(MJD) ) {
    806811        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    807812                "Time conversion to MJD failed.  Invalid input time.\n");
    … …  
    858863                k = 2;
    859864            }
    860             if (k > (numRows-2)) {
    861                 k = numRows-2;
    862             }
     865            //            if (k > (numRows-2)) {
     866            //                k = numRows-2;
     867            //            }
    863868            for (int m = k-1; m <= k+2; m++) {
    864869                xTerm = X->data.F64[m];
    … …  
    935940    // Convert psTime to MJD
    936941    double MJD = psTimeToMJD(time);
    937     if (MJD == NAN) {
     942    if (isnan(MJD)) {
    938943        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    939944                "Time conversion to MJD failed.  Invalid input time.\n");
    … …  
    10191024    double MJD = psTimeToMJD(time);
    10201025    //    printf("\nMJD check = %.13g\n", MJD);
    1021     if (MJD == NAN) {
     1026    if (isnan(MJD)) {
    10221027        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    10231028                "Time conversion to MJD failed.  Invalid input time.\n");
    … …  
    12521257    // Calculate Julian centuries
    12531258    //If either input time is NULL, assume it to be J2000 -> from SDRS as of rev 18
     1259    psTime *from = NULL;
     1260    psTime *to = NULL;
    12541261    if (fromTime == NULL) {
    12551262        fromMJD = MJD_2000;
     1263        from = psTimeFromMJD(fromMJD);
    12561264    } else {
    12571265        fromMJD = psTimeToMJD(fromTime);
     1266        from = psTimeCopy(fromTime);
    12581267    }
    12591268    if (toTime == NULL) {
    12601269        toMJD = MJD_2000;
     1270        to = psTimeFromMJD(toMJD);
    12611271    } else {
    12621272        toMJD = psTimeToMJD(toTime);
    1263     }
    1264     psTime *from = NULL;
    1265     psTime *to = NULL;
     1273        to = psTimeCopy(toTime);
     1274    }
    12661275    if (fromMJD > toMJD) {
    12671276        psWarning("From time is later than to time in psSpherePrecess.\n");
    … …  
    12901299        // Create transform with proper constants
    12911300        psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP);
     1301        psFree(from);
     1302        psFree(to);
    12921303        return tmpST;
    12931304    } else if (mode == PS_PRECESS_IAU2000A) {
    … …  
    12961307        //rotation corresponding to fromTime to the toTime rotation.
    12971308
    1298         //Since the time inputs are allowed to be NULL, either convert the MJD time
    1299         //or copy to non-NULL time
    1300         if (fromTime == NULL) {
    1301             from = psTimeFromMJD(fromMJD);
    1302         } else {
    1303             from = psTimeCopy(fromTime);
    1304         }
    1305         if (toTime == NULL) {
    1306             to = psTimeFromMJD(toMJD);
    1307         } else {
    1308             to = psTimeCopy(toTime);
    1309         }
    13101309        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
    13111310        //Combine the quaternions to produce the output psSphereRot.
    … …  
    13291328        //for COMPLETE_A come from the IERS Bulletin A.
    13301329
    1331         //Since the time inputs are allowed to be NULL, either convert the MJD time
    1332         //or copy to non-NULL time
    1333         if (fromTime == NULL) {
    1334             from = psTimeFromMJD(fromMJD);
    1335         } else {
    1336             from = psTimeCopy(fromTime);
    1337         }
    1338         if (toTime == NULL) {
    1339             to = psTimeFromMJD(toMJD);
    1340         } else {
    1341             to = psTimeCopy(toTime);
    1342         }
    13431330        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
    13441331        //Add in the precession corrections from IERS bulletin A.
    … …  
    13711358        //for COMPLETE_B come from the IERS Bulletin B.
    13721359
    1373         //Since the time inputs are allowed to be NULL, either convert the MJD time
    1374         //or copy to non-NULL time
    1375         if (fromTime == NULL) {
    1376             from = psTimeFromMJD(fromMJD);
    1377         } else {
    1378             from = psTimeCopy(fromTime);
    1379         }
    1380         if (toTime == NULL) {
    1381             to = psTimeFromMJD(toMJD);
    1382         } else {
    1383             to = psTimeCopy(toTime);
    1384         }
    1385 
    13861360        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
    13871361        //Add in the precession corrections from IERS bulletin B.

  • trunk/psLib/test/astro/tap_psEarthOrientation_corrections.c

    r9948 r9976  
    1818#include "pstap.h"
    1919
    20 static void testSphereRotCreate(void);
    21 //static void testSphereRotConvert(void);
    22 //static void testSphereOffsets(void);
     20static void testEOCInit(void);
     21static void testAberration(void);
     22static void testGravDef(void);
     23static void testEOC_Corrs(void);
     24
     25#define MJD_1900  15021.0        // Modified Julian Day 1/1/1900 00:00:00
     26#define MJD_2100  88069.0        // Modified Julian Day 1/1/2100 00:00:00
    2327
    2428int main(void)
    2529{
    26     plan_tests(1);
     30    plan_tests(8);
    2731
    2832    diag("Tests for psEarthOrientation Correction Functions");
    2933
    30     testSphereRotCreate();
    31     //    testSphereRotConvert();
     34    // Initialize library internal structures
     35    psLibInit("pslib.config");
     36
     37    testEOCInit();
     38    testAberration();
     39    testGravDef();
     40    testEOC_Corrs();
    3241    //    testSphereOffsets();
     42
     43    // Cleanup library
     44    psLibFinalize();
    3345
    3446    done();
    3547}
    3648
    37 void testSphereRotCreate(void)
     49void testEOCInit(void)
    3850{
    39     diag("  >>>Test 1:  psSphereRot Creation Functions");
     51    diag("  >>>Test 1:  p_psEOCInit ");
     52    //Test for psEarthPoleAlloc
     53    //Return properly allocated psEarthPole
     54    /*    {
     55            skip_start(  ep == NULL, 3,
     56                         "Skipping 3 tests because psEarthPole is NULL!");
     57            ep = psEarthPoleAlloc();
     58            ok( ep != NULL,
     59                "psEarthPoleAlloc:               return properly allocated psEarthPole.");
     60            psFree(ep);
     61            ep = NULL;
     62        }
     63    */
     64
     65
     66    //Check for Memory leaks
     67    {
     68        checkMem();
     69    }
     70
     71}
     72
     73void testAberration(void)
     74{
     75    diag("  >>>Test 2:  psAberration ");
     76    psSphere *apparent = NULL;
     77    psSphere *empty = NULL;
     78    psCube *actualCube = psCubeAlloc();
     79    //values from After gravity deflection//
     80    actualCube->x = -0.35961949760293604;
     81    actualCube->y = 0.5555613950298085;
     82    actualCube->z = 0.7496835020836093;
     83    psSphere *actual = psCubeToSphere(actualCube);
     84    psCube *cubeDir = psCubeAlloc();
     85    cubeDir->x = 5148.713262821658;
     86    cubeDir->y = -26945.04752348012;
     87    cubeDir->z = -11682.787302030947;
     88    cubeDir->x += -357.6031690489248;
     89    cubeDir->y += 248.46429758174693;
     90    cubeDir->z += 0.09694774143797581;
     91    psSphere *direction = psCubeToSphere(cubeDir);
     92    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
     93    double c = 299792458.0; // Speed of light in vacuum (src:NIST) m/s
     94    speed /= c;
     95
     96    //Tests for psAberration
     97    //Return NULL for NULL actual sphere input
     98    {
     99        empty = psAberration(apparent, NULL, direction, speed);
     100        ok( empty == NULL,
     101            "psAberration:                   return NULL for NULL actual sphere input.");
     102    }
     103    //Return NULL for NULL direction sphere input
     104    {
     105        empty = psAberration(apparent, actual, NULL, speed);
     106        ok( empty == NULL,
     107            "psAberration:                   return NULL for NULL direction sphere input.");
     108    }
     109    //Return NULL for speed = 0.0
     110    {
     111        empty = psAberration(apparent, actual, direction, 0.0);
     112        ok( empty == NULL,
     113            "psAberration:                   return NULL for zero speed input.");
     114    }
     115    //Return NULL for speed > 1.0
     116    {
     117        empty = psAberration(apparent, actual, direction, 1.1);
     118        ok( empty == NULL,
     119            "psAberration:                   return NULL for impossible speed input (> c).");
     120    }
     121
     122    //Return correct values for valid inputs
     123    {
     124        double x, y, z;
     125        x = -0.35963388069046304;
     126        y = 0.5555192509816625;
     127        z = 0.7497078321908413;
     128        apparent = psAberration(apparent, actual, direction, speed);
     129        skip_start(  apparent == NULL, 3,
     130                     "Skipping 3 tests because psEarthPole is NULL!");
     131        psCube *outCube = psSphereToCube(apparent);
     132        ok_double_tol( outCube->x, x, 0.001,
     133                       "psAberration:                   return correct sphere for valid inputs.");
     134        ok_double_tol( outCube->y, y, 0.001,
     135                       "psAberration:                   return correct sphere for valid inputs.");
     136        ok_double_tol( outCube->z, z, 0.001,
     137                       "psAberration:                   return correct sphere for valid inputs.");
     138        psFree(outCube);
     139        skip_end();
     140    }
     141
     142    //Check for Memory leaks
     143    {
     144        psFree(actualCube);
     145        psFree(cubeDir);
     146        psFree(apparent);
     147        psFree(actual);
     148        psFree(direction);
     149        checkMem();
     150    }
     151
     152}
     153
     154void testGravDef(void)
     155{
     156    diag("  >>>Test 3:  psGravityDeflection ");
     157    //Test for psGravityDeflection
     158    //Return properly allocated psEarthPole
    40159    /*
    41     //Tests for psSphereRotAlloc
    42     //Return NULL for NAN input
    43     {
    44         psSphereRot *s1 = psSphereRotAlloc(ALPHA_P, NAN, PHI_P);
    45         ok( s1 == NULL,
    46             "psSphereRotAlloc:               return NULL for NAN input.");
    47     }
    48     //Tests for psMemCheckSphereRot
    49     //Make sure psMemCheckSphereRot works correctly - return false
    50     {
    51         int j = 2;
    52         ok( !psMemCheckSphereRot(&j),
    53              "psMemCheckSphereRot:            return false for non-SphereRot input.");
    54         skip_start(  s1 == NULL, 1,
    55                      "Skipping 1 tests because psSphereRot is NULL!");
    56         skip_end();
    57     }
     160        psSphere *apparent = NULL;
     161        psSphere *empty = NULL;
     162        psCube *actualCube = psCubeAlloc();
     163        actualCube->x = -0.3596195125758298;
     164        actualCube->y = 0.5555613903455866;
     165        actualCube->z = 0.7496834983724809;
     166        psSphere *actual = psCubeToSphere(actualCube);
     167        psCube *sunCube = psCubeAlloc();
     168        sunCube->x = 1.467797790127511e11;
     169        sunCube->y = 2.5880956908748722e10;
     170        sunCube->z = 1.1220046291457653e10;
     171        double sunLength = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z);
     172        sunCube->x /= sunLength;
     173        sunCube->y /= sunLength;
     174        sunCube->z /= sunLength;
     175        if (VERBOSE) {
     176        printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n",
     177        sunCube->x, sunCube->y, sunCube->z);
     178    }
     179        psSphere *sun = psCubeToSphere(sunCube);
     180        if (VERBOSE) {
     181        printf("sunSphere  = r, d = %.13g, %.13g\n", sun->r, sun->d);
     182    }
     183        psCube *outCube = psCubeAlloc();
     184     
     185        empty = psGravityDeflection(apparent, empty, sun);
     186        if (empty != NULL) {
     187        psError(PS_ERR_BAD_PARAMETER_NULL, true,
     188        "psGravityDeflection Failed to return NULL for NULL actual input sphere.\n");
     189        return 1;
     190    }
     191        empty = psGravityDeflection(apparent, actual, empty);
     192        if (empty != NULL) {
     193        psError(PS_ERR_BAD_PARAMETER_NULL, true,
     194        "psGravityDeflection Failed to return NULL for NULL sun input sphere.\n");
     195        return 2;
     196    }
     197     
     198        apparent = psGravityDeflection(NULL, actual, sun);
     199        //    apparent->r *= -1.0;
     200        //    apparent->d *= -1.0;
     201        psSphere *result2;
     202        //    psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     203        //    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     204        //    psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN);
     205        if (VERBOSE) {
     206        printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     207        actualCube->x, actualCube->y, actualCube->z);
     208    }
     209        //    psCube *outCube = psSphereToCube(result);
     210        //    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     211        //           outCube->x, outCube->y, outCube->z);
     212        psCube *outCube2 = psSphereToCube(apparent);
     213        if (VERBOSE) {
     214        printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g  -- \n",
     215        outCube2->x, outCube2->y, outCube2->z);
     216    }
     217        double x, y, z;
     218        x = -0.35961949760293604;
     219        y = 0.5555613950298085;
     220        z = 0.7496835020836093;
     221     
     222        if (VERBOSE) {
     223        printf(" -- expectCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
     224     
     225            //    if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
     226            //            fabs(z - outCube->z) > DBL_EPSILON ) {
     227            //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     228            //                "psGravityDeflection returned incorrect values.\n");
     229        printf("expect-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
     230        (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) );
     231        printf("expect-result=  x,y,z  = %.13g, %.13g, %.13g \n",
     232        (x - outCube2->x), (y - outCube2->y), (z - outCube2->z) );
     233        printf("result-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
     234        (outCube2->x - actualCube->x), (outCube2->y - actualCube->y),
     235        (outCube2->z - actualCube->z) );
     236    }
     237        //        return 1;
     238        //    }
     239        //    psFree(result2);
     240        outCube2->x = x;
     241        outCube2->y = y;
     242        outCube2->z = z;
     243        result2 = psCubeToSphere(outCube2);
     244        //    psFree(result);
     245        psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
     246        psFree(result2);
     247        result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     248        if (VERBOSE) {
     249        printf("The apparent output sphere = r,d = %.13g, %.13g\n", result2->r, result2->d);
     250        printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d);
     251    }
     252        psFree(result2);
     253     
     254        psFree(outCube);
     255        psFree(outCube2);
     256        psFree(sunCube);
     257        psFree(actualCube);
     258        psFree(result);
     259        psFree(actual);
     260        psFree(apparent);
     261        psFree(sun);
     262     
     263     
     264     
     265        {
     266            skip_start(  ep == NULL, 3,
     267                         "Skipping 3 tests because psEarthPole is NULL!");
     268            ep = psEarthPoleAlloc();
     269            ok( ep != NULL,
     270                "psEarthPoleAlloc:               return properly allocated psEarthPole.");
     271            psFree(ep);
     272            ep = NULL;
     273        }
    58274    */
     275
     276
    59277    //Check for Memory leaks
    60278    {
    … …  
    64282}
    65283
    66 
    67 
     284void testEOC_Corrs(void)
     285{
     286    diag("  >>>Test 4:  psEOC Correction Functions");
     287    //Tests for psEOC_PrecessionCorr
     288    /*
     289        psTime *empty = NULL;
     290        psTime *time2 = psTimeAlloc(PS_TIME_UTC);
     291        time2->sec = timesec;
     292        time2->nsec = 0;
     293        time2->leapsecond = false;
     294        //    time2 = psTimeConvert(time2, PS_TIME_TAI);
     295     
     296        //Tests for Precession Correction function//
     297        //Return NULL for NULL time input
     298        psEarthPole *pcorr = NULL;
     299        psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
     300        pcorr = psEOC_PrecessionCorr(empty, PS_IERS_A);
     301        if (pcorr != NULL) {
     302            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     303                    "psEOC_PrecessionCorr failed to return NULL for NULL time input.\n");
     304            return 5;
     305        }
     306     
     307        //Return NULL for Invalid IERS table
     308        psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
     309        pcorr = psEOC_PrecessionCorr(time2, 3);
     310        if (pcorr != NULL) {
     311            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     312                    "psEOC_PrecessionCorr failed to return NULL for incorrect IERS table.\n");
     313            return 6;
     314        }
     315        psFree(pcorr);
     316     
     317        //Check values from IERS table A
     318        pcorr = psEOC_PrecessionCorr(time2, PS_IERS_A);
     319        if ( pcorr == NULL ) {
     320            psError(PS_ERR_BAD_PARAMETER_NULL, false,
     321                    "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
     322            return 7;
     323        } else {
     324            if (VERBOSE) {
     325                printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g,  %.13g, %.13g\n",
     326                       pcorr->x, pcorr->y, pcorr->s);
     327            }
     328        }
     329        psFree(pcorr);
     330     
     331        //Check values from IERS table B
     332        pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B);
     333        if ( pcorr == NULL ) {
     334            psError(PS_ERR_BAD_PARAMETER_NULL, false,
     335                    "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
     336            return 9;
     337        } else {
     338            double xx, yy, ss;
     339            xx = 0.06295703125;
     340            yy = -0.0287618408203125;
     341            ss = 0.0;
     342            xx = SEC_TO_RAD(xx) * 1e-3;
     343            yy = SEC_TO_RAD(yy) * 1e-3;
     344            //        if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
     345            //                || fabs(pcorr->s - ss) > DBL_EPSILON) {
     346            //            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     347            //                    "   psEOC_PrecessionCorr return incorrect values.\n");
     348            if (VERBOSE) {
     349                printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     350                       pcorr->x, pcorr->y, pcorr->s);
     351                printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n", xx, yy, ss);
     352                printf("          A difference of:              %.13g,  %.13g, %.13g\n\n",
     353                       (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
     354            }
     355            //            return 10;
     356            //        }
     357        }
     358        //precess is the *actual* output from PrecessionModel + PrecessionCorr
     359        psEarthPole *precess = psEOC_PrecessionModel(time2);
     360        precess->x += pcorr->x;
     361        precess->y += pcorr->y;
     362        double xCorr, yCorr;
     363        xCorr = 3.05224300720406e-10;
     364        yCorr = -1.39441339235822e-10;
     365        //pcorr is the *expected* output from PrecessionModel//
     366        pcorr->x = 2.857175590089105e-4;
     367        pcorr->y = 2.3968739377734732e-5;
     368        pcorr->s = -1.3970066457904322e-8;
     369        pcorr->x += xCorr;
     370        pcorr->y += yCorr;
     371        psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess);
     372        psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv);
     373        double q0, q1, q2;
     374        q0 = -1.1984522406756289e-5;
     375        q1 = 1.4285893358610674e-4;
     376        q2 = 1.2191193518914336e-10;
     377        psSphereRot *pni = psSphereRot_CEOtoGCRS(pcorr);
     378        if (fabs(pni->q0-q0) > FLT_EPSILON || fabs(pni->q1-q1) > FLT_EPSILON ||
     379                fabs(pni->q2-q2) > FLT_EPSILON ) {
     380            printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n");
     381        }
     382        //    printf("  Output from CEOtoGCRS only  = %.13g,%.13g,%.13g,%.13g\n",
     383        //           pni->q0, pni->q1, pni->q2, pni->q3);
     384        if (VERBOSE) {
     385            printf("  Expected sphere rotation    = %.13g, %.13g, %.13g\n", q0,q1,q2);
     386            printf("  Result sphere rotation      = %.13g, %.13g, %.13g\n",
     387                   precessNutInv->q0, precessNutInv->q1, precessNutInv->q2);
     388            printf("     Difference         =        %.13g, %.13g, %.13g\n\n",
     389                   precessNutInv->q0-q0, precessNutInv->q1-q1, precessNutInv->q2-q2);
     390        }
     391        psCube *objC = psCubeAlloc();
     392        //    objC->x = -3.5963388069046304;
     393        //    objC->y = 0.5555192509816625;
     394        //    objC->z = 0.7497078321908413;
     395        //    objSetup();
     396     
     397        //This is the sphere rotation for the *expected* precession output//
     398        psSphereRot *pn = psSphereRotConjugate(NULL, pni);
     399     
     400        //    psSphere *sphere = psSphereAlloc();
     401        //    *sphere = *obj;
     402        //    psFree(obj);
     403     
     404        //create a psSphere for (from) the start position given in eoc_testing//
     405        objC->x = -0.35963388069046304;
     406        objC->y = 0.5555192509816625;
     407        objC->z = 0.7497078321908413;
     408        psSphere *sphere = psCubeToSphere(objC);
     409     
     410        psSphere *expect = psSphereRotApply(NULL, pn, sphere);
     411        //expected results below - stored in:  sphere  //
     412        double x,y,z;
     413        x = -0.3598480726985338;
     414        y = 0.5555012823608123;
     415        z = 0.7496183628158023;
     416        if (VERBOSE) {
     417            printf("\n<<Expected out       = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
     418        }
     419        //    psFree(objC);
     420        //    objC = psSphereToCube(expect);
     421        //printf("<<Expected out (CEO)  = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     422        //printf("     Difference     =           %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
     423        //    x = objC->x;
     424        //    y = objC->y;
     425        //    z = objC->z;
     426        psSphere *result = psSphereRotApply(NULL, precessNut, sphere);
     427        psFree(objC);
     428        objC = psSphereToCube(result);
     429        double xx = greatCircle(result, expect);
     430        if (VERBOSE) {
     431            printf("<<Resulting out      = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     432            printf("     Difference         =      %.13g, %.13g, %.13g\n\n",
     433                   objC->x-x, objC->y-y, objC->z-z);
     434            printf("GREAT CIRCLE DIFFERENCE = %.13g \n", xx);
     435        }
     436     
     437        psFree(precess);
     438        psFree(precessNut);
     439        psFree(precessNutInv);
     440        psFree(expect);
     441        psFree(objC);
     442     
     443        psFree(sphere);
     444        psFree(result);
     445        psFree(pn);
     446        psFree(pni);
     447        psFree(pcorr);
     448        psFree(time2);
     449        if (!p_psEOCFinalize() ) {
     450            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
     451            return 12;
     452        }
     453     
     454    */
     455
     456    //Tests for psEOC_PolarTideCorr
     457    /*
     458        psTime *in = psTimeAlloc(PS_TIME_UTC);
     459        in->sec = timesec;
     460        in->nsec = 0;
     461        in->leapsecond = false;
     462        psTime *empty = NULL;
     463        psEarthPole *eop = NULL;
     464     
     465        //Return NULL for NULL input time
     466        psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
     467        eop = psEOC_PolarTideCorr(empty);
     468        if (eop != NULL) {
     469        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     470        "psEOC_PolarTideCorr failed to return NULL for NULL input time.\n");
     471        return 1;
     472    }
     473     
     474        eop = psEOC_PolarTideCorr(in);
     475        if (eop == NULL) {
     476        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     477        "psEOC_PolarTideCorr returned NULL for valid input time.\n");
     478        return 2;
     479    } else {
     480        if (VERBOSE) {
     481        printf("\nPolarTideCorr output = x,y,s = %.13g, %.13g, %.13g\n",
     482        eop->x, eop->y, eop->s);
     483    }
     484    }
     485     
     486        psFree(in);
     487        psFree(eop);
     488     
     489    */
     490
     491
     492
     493    //Tests for psEOC_NutationCorr
     494    /*
     495        psTime *in = psTimeAlloc(PS_TIME_UTC);
     496        in->sec = timesec;
     497        in->nsec = 0;
     498        in->leapsecond = false;
     499        psTime *empty = NULL;
     500        psEarthPole *nute = NULL;
     501     
     502        //Return NULL for NULL input time.
     503        psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
     504        nute = psEOC_NutationCorr(empty);
     505        if (nute != NULL) {
     506        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     507        "psEOC_NutationCorr failed to return NULL for NULL input time.\n");
     508        return 1;
     509    }
     510        //Return NULL for UT1 time input
     511        *//*    psTime *UT1time = psTimeAlloc(PS_TIME_UT1);
     512        nute = psEOC_NutationCorr(UT1time);
     513        if (nute != NULL) {
     514        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     515        "psEOC_NutationCorr failed to return NULL for UT1 input time.\n");
     516        return 2;
     517    }
     518        psFree(UT1time);
     519        *//*
     520        //Check return values from valid nutation time input
     521        nute = psEOC_NutationCorr(in);
     522        if ( nute == NULL ) {
     523            psError(PS_ERR_BAD_PARAMETER_NULL, false,
     524                    "psEOC_NutationCorr returned NULL for valid input.\n");
     525            return 3;
     526        } else {
     527            if (VERBOSE) {
     528                printf("Nutation Correction output = x,y,s = %.13g, %.13g, %.13g\n\n",
     529                       nute->x, nute->y, nute->s);
     530            }
     531        }
     532        psFree(nute);
     533        psFree(in);
     534     
     535        if (!p_psEOCFinalize() ) {
     536            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
     537            return 12;
     538        }
     539     
     540    */
     541    //Return properly allocated psEarthPole
     542    /*
     543     
     544     
     545        {
     546            skip_start(  ep == NULL, 3,
     547                         "Skipping 3 tests because psEarthPole is NULL!");
     548            ep = psEarthPoleAlloc();
     549            ok( ep != NULL,
     550                "psEarthPoleAlloc:               return properly allocated psEarthPole.");
     551            psFree(ep);
     552            ep = NULL;
     553        }
     554    */
     555
     556
     557    //Check for Memory leaks
     558    {
     559        checkMem();
     560    }
     561
     562}
     563
     564

  • trunk/psLib/test/astro/tap_psEarthOrientation_motion.c

    r9948 r9976  
    44 * Description:  Tests for psEarthPoleAlloc, psEOC_PrecessionModel, psSpherePrecess
    55 *                         psSphereRot_CEOtoGCRS, psSphereRot_TEOtoCEO,
    6  *                         psEOC_GetPolarMotion, psSphereRot_ITRStoTEO,
     6 *                         psEOC_GetPolarMotion, psSphereRot_ITRStoTEO
    77 *
    88 * Author: dRob <David.Robbins@mhpcc.hpc.mil>, (C) 2006
    … …  
    1818#include "pstap.h"
    1919
    20 static void testSphereRotCreate(void);
    21 //static void testSphereRotConvert(void);
    22 //static void testSphereOffsets(void);
     20static void testPrecessionModel(void);
     21static void testPolarMotion(void);
     22static void testSphereRots(void);
     23static void testSpherePrecess(void);
     24
     25#define MJD_1900  15021.0        // Modified Julian Day 1/1/1900 00:00:00
     26#define MJD_2100  88069.0        // Modified Julian Day 1/1/2100 00:00:00
    2327
    2428int main(void)
    2529{
    26     plan_tests(1);
     30    plan_tests(48);
    2731
    2832    diag("Tests for psEarthOrientation Motion Functions");
    2933
    30     testSphereRotCreate();
    31     //    testSphereRotConvert();
    32     //    testSphereOffsets();
     34    // Initialize library internal structures
     35    psLibInit("pslib.config");
     36
     37    testPrecessionModel();
     38    testPolarMotion();
     39    testSphereRots();
     40    testSpherePrecess();
     41
     42    // Cleanup library
     43    psLibFinalize();
    3344
    3445    done();
    3546}
    3647
    37 void testSphereRotCreate(void)
     48void testPrecessionModel(void)
    3849{
    39     diag("  >>>Test 1:  psSphereRot Creation Functions");
     50    diag("  >>>Test 1:  psEOC_PrecessionModel");
     51
     52    psEarthPole *ep = NULL;
     53    psTime *testTime = NULL;
     54
     55    //Test for psEarthPoleAlloc
     56    //Return properly allocated psEarthPole
     57    {
     58        ep = psEarthPoleAlloc();
     59        ok( ep != NULL,
     60            "psEarthPoleAlloc:               return properly allocated psEarthPole.");
     61        psFree(ep);
     62        ep = NULL;
     63    }
     64    //Tests for psEOC_PrecessionModel
     65    //Return NULL for NULL time input.
     66    {
     67        ep = psEOC_PrecessionModel(NULL);
     68        ok( ep == NULL,
     69            "psEOC_PrecessionModel:          return NULL for NULL time input.");
     70    }
     71    //Return NULL for UT1 time input.
     72    {
     73        testTime = psTimeAlloc(PS_TIME_UT1);
     74        ep = psEOC_PrecessionModel(testTime);
     75        ok( ep == NULL,
     76            "psEOC_PrecessionModel:          return NULL for UT1 time input.");
     77    }
     78    //Return NULL for invalid time input.
     79    {
     80        psFree(testTime);
     81        testTime = psTimeAlloc(PS_TIME_UTC);
     82        testTime->nsec = 2e9;
     83        ep = psEOC_PrecessionModel(testTime);
     84        ok( ep == NULL,
     85            "psEOC_PrecessionModel:          return NULL for invalid time input.");
     86        testTime->nsec = 0;
     87    }
     88    //Return NULL for failed eoc init - missing file
     89    //    printf("\n filename = '%s' \n", p_psGetConfigFileName());
     90    {
     91        testTime->sec = 1049160600;
     92        testTime->nsec = 0;
     93        testTime->leapsecond = false;
     94        skip_start(rename("../../etc/pslib.config", "../../etc/pslib_config.bak"),
     95                   1, "Skipping 1 tests because file rename failed!");
     96        ep = psEOC_PrecessionModel(testTime);
     97        ok( ep == NULL,
     98            "psEOC_PrecessionModel:          return NULL for failed eoc init.");
     99        rename("../../etc/pslib_config.bak", "../../etc/pslib.config");
     100        skip_end();
     101    }
     102    //Return valid EarthPole for valid time input
     103    {
     104        double x, y, s;
     105        x = 2.857175590089105e-4;
     106        y = 2.3968739377734732e-5;
     107        s = -1.3970066457904322e-8;
     108        ep = psEOC_PrecessionModel(testTime);
     109        skip_start(  ep == NULL, 3,
     110                     "Skipping 3 tests because psEarthPole is NULL!");
     111        ok_double_tol(ep->x, x, 0.1,
     112                      "psEOC_PrecessionModel:          return valid EarthPole for valid inputs (x).");
     113        ok_double_tol(ep->y, y, 0.1,
     114                      "psEOC_PrecessionModel:          return valid EarthPole for valid inputs (y).");
     115        ok_double_tol(ep->s, s, 0.1,
     116                      "psEOC_PrecessionModel:          return valid EarthPole for valid inputs (s).");
     117        skip_end();
     118    }
     119
     120
    40121    /*
    41     //Tests for psSphereRotAlloc
    42     //Return NULL for NAN input
    43     {
    44     psSphereRot *s1 = psSphereRotAlloc(ALPHA_P, NAN, PHI_P);
    45     ok( s1 == NULL,
    46     "psSphereRotAlloc:               return NULL for NAN input.");
    47     }
    48122    //Tests for psMemCheckSphereRot
    49123    //Make sure psMemCheckSphereRot works correctly - return false
    … …  
    59133    //Check for Memory leaks
    60134    {
     135        psFree(testTime);
     136        psFree(ep);
    61137        checkMem();
    62138    }
    … …  
    64140}
    65141
     142void testPolarMotion(void)
     143{
     144    diag("  >>>Test 2:  psEOC_GetPolarMotion ");
     145
     146    psTime *in = psTimeAlloc(PS_TIME_UTC);
     147    in->sec = 1049160600;
     148    in->nsec = 0;
     149    in->leapsecond = false;
     150    psEarthPole *polarMotion = NULL;
     151
     152    //Tests for psEOC_GetPolarMotion
     153    //Return NULL for NULL time input.
     154    {
     155        polarMotion = psEOC_GetPolarMotion(NULL, PS_IERS_B);
     156        ok( polarMotion == NULL,
     157            "psEOC_GetPolarMotion:           return NULL for NULL time input.");
     158    }
     159    //Return NULL for invalid IERS table input
     160    {
     161        polarMotion = psEOC_GetPolarMotion(NULL, PS_IERS_B+1);
     162        ok( polarMotion == NULL,
     163            "psEOC_GetPolarMotion:          return NULL for invalid IERS table input.");
     164    }
     165    //Return NULL for invalid time input.
     166    {
     167        in->nsec = 2e9;
     168        polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
     169        ok( polarMotion == NULL,
     170            "psEOC_GetPolarMotion:          return NULL for invalid time input.");
     171        in->nsec = 0;
     172    }
     173    //Return NULL for failed eoc init - missing file
     174    {
     175        p_psEOCFinalize();
     176        skip_start(rename("../../etc/pslib.config", "../../etc/pslib_config.bak"),
     177                   1, "Skipping 1 tests because file rename failed!");
     178        polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
     179        ok( polarMotion == NULL,
     180            "psEOC_GetPolarMotion:          return NULL for failed eoc init.");
     181        rename("../../etc/pslib_config.bak", "../../etc/pslib.config");
     182        skip_end();
     183    }
     184    //Return valid EarthPole for valid time input - IERS B
     185    {
     186        double x, y, s;
     187        x = -6.454389659777e-07;
     188        y = 2.112606414597e-06;
     189        s = 0.0;
     190        polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
     191        skip_start(  polarMotion == NULL, 3,
     192                     "Skipping 3 tests because psEarthPole is NULL!");
     193        ok_double_tol(polarMotion->x, x, 0.1,
     194                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     195                      "(x) - IERS B.");
     196        ok_double_tol(polarMotion->y, y, 0.1,
     197                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     198                      "(y) - IERS B.");
     199        ok_double_tol(polarMotion->s, s, 0.1,
     200                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     201                      "(s) - IERS B.");
     202        skip_end();
     203        psFree(polarMotion);
     204    }
     205    //Return valid EarthPole for valid time input - IERS A
     206    {
     207        double x, y, s;
     208        x = -6.45381397904e-07;
     209        y = 2.112819726698e-06;
     210        s = 0.0;
     211        polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A);
     212        skip_start(  polarMotion == NULL, 3,
     213                     "Skipping 3 tests because psEarthPole is NULL!");
     214        ok_double_tol(polarMotion->x, x, 0.1,
     215                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     216                      "(x) - IERS A.");
     217        ok_double_tol(polarMotion->y, y, 0.1,
     218                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     219                      "(y) - IERS A.");
     220        ok_double_tol(polarMotion->s, s, 0.1,
     221                      "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs "
     222                      "(s) - IERS A.");
     223        skip_end();
     224        psFree(polarMotion);
     225    }
     226    //Return valid EarthPole for valid time input - IERS A
     227    {
     228        psTime *firstTime = psTimeFromMJD(41684.50);
     229        polarMotion = psEOC_GetPolarMotion(firstTime, PS_IERS_B);
     230        ok( polarMotion != NULL,
     231            "psEOC_GetPolarMotion:          return valid EarthPole for valid inputs.");
     232        psFree(firstTime);
     233    }
     234
     235
     236    //Check for Memory leaks
     237    {
     238        psFree(in);
     239        psFree(polarMotion);
     240        checkMem();
     241    }
     242
     243}
     244
     245static psSphere *objSetup(void)
     246{
     247    psSphere *obj = NULL;
     248    psCube *tempCube = psCubeAlloc();
     249    tempCube->x = -0.3598480726985338;
     250    tempCube->y = 0.5555012823608123;
     251    tempCube->z = 0.7496183628158023;
     252    obj = psCubeToSphere(tempCube);
     253    psFree(tempCube);
     254    return obj;
     255}
     256
     257void testSphereRots(void)
     258{
     259    diag("  >>>Test 3:  psSphereRot Functions");
     260    psSphereRot *out = NULL;
     261    psEarthPole *in = NULL;
     262    double q0,q1,q2,q3;
     263    q0 = -1.1984522406756289e-5;
     264    q1 = 1.4285893358610674e-4;
     265    q2 = 1.2191193518914336e-10;
     266    q3 = -0.9999999897238481;
     267
     268    //Tests for psSphereRot_CEOtoGCRS
     269    //Return NULL for NULL earthPole input
     270    {
     271        out = psSphereRot_CEOtoGCRS(in);
     272        ok( out == NULL,
     273            "psSphereRot_CEOtoGCRS:         return NULL for NULL earthPole input.");
     274        in = psEarthPoleAlloc();
     275    }
     276    in->x = 2.857175590089105e-4;
     277    in->y = 2.3968739377734732e-5;
     278    in->s = -1.3970066457904322e-8;
     279    //Return correct psSphereRot for valid input
     280    {
     281        out = psSphereRot_CEOtoGCRS(in);
     282        skip_start(  out == NULL, 4,
     283                     "Skipping 4 tests because psSphereRot output is NULL!");
     284        ok_double_tol( out->q0, q0, 0.0001,
     285                       "psSphereRot_CEOtoGCRS:         return correct psSphereRot for valid input (q0).");
     286        ok_double_tol( out->q1, q1, 0.0001,
     287                       "psSphereRot_CEOtoGCRS:         return correct psSphereRot for valid input (q1).");
     288        ok_double_tol( out->q2, q2, 0.0001,
     289                       "psSphereRot_CEOtoGCRS:         return correct psSphereRot for valid input (q2).");
     290        ok_double_tol( out->q3, -q3, 0.0001,
     291                       "psSphereRot_CEOtoGCRS:         return correct psSphereRot for valid input (q3).");
     292        skip_end();
     293    }
     294
     295    //Tests for psSphereRot_TEOtoCEO
     296    psFree(out);
     297    psTime *time = psTimeAlloc(PS_TIME_UT1);
     298    psTime *time2 = psTimeAlloc(PS_TIME_UTC);
     299    time->sec = 1049160600-1;
     300    time->nsec = 657017200;
     301    time->leapsecond = false;
     302    psEarthPole *temp = psEarthPoleAlloc();
     303    temp->s = -2.0;
     304    //Return NULL for NULL time input
     305    {
     306        out = psSphereRot_TEOtoCEO(NULL, NULL);
     307        ok( out == NULL,
     308            "psSphereRot_TEOtoCEO:          return NULL for NULL time input.");
     309    }
     310    //Return NULL for invalid time input - large nsec, UTC time type
     311    {
     312        time2->nsec = 3e9;
     313        out = psSphereRot_TEOtoCEO(time2, NULL);
     314        ok( out == NULL,
     315            "psSphereRot_TEOtoCEO:          return NULL for invalid time input.");
     316    }
     317    //Return NULL for invalid time input - UT1 time type, large nsec
     318    {
     319        psFree(time2);
     320        time2 = psTimeAlloc(PS_TIME_UT1);
     321        time2->nsec = 1e9;
     322        temp->s = 1.0;
     323        out = psSphereRot_TEOtoCEO(time2, temp);
     324        ok( out == NULL,
     325            "psSphereRot_TEOtoCEO:          return NULL for invalid time input.");
     326    }
     327    //Return NULL for invalid earthPole input
     328    {
     329        time2->nsec = 0;
     330        psEarthPole *temp = psEarthPoleAlloc();
     331        temp->s = -2.0;
     332        out = psSphereRot_TEOtoCEO(time2, temp);
     333        ok( out == NULL,
     334            "psSphereRot_TEOtoCEO:          return NULL for invalid earthPole input.");
     335    }
     336    //Return correct psSphereRot for valid inputs
     337    double x, y, z;
     338    x = 0.01698625430807123;
     339    y = -0.6616523084626379;
     340    z = 0.7496183628158023;
     341    {
     342        psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
     343        out = psSphereRot_TEOtoCEO(time, polarTideCorr);
     344        skip_start(  out == NULL, 3,
     345                     "Skipping 3 tests because psSphereRot output is NULL!");
     346        psSphere *obj = objSetup();
     347        psSphereRot *earthRot = psSphereRotConjugate(NULL, out);
     348        psSphere *result = psSphereRotApply(NULL, earthRot, obj);
     349        psCube *cube = psSphereToCube(result);
     350        ok_double_tol( cube->x, x, 0.0001,
     351                       "psSphereRot_TEOtoCEO:          return NULL for NULL time input. (x)");
     352        ok_double_tol( cube->y, y, 0.0001,
     353                       "psSphereRot_TEOtoCEO:          return NULL for NULL time input. (y)");
     354        ok_double_tol( cube->z, z, 0.0001,
     355                       "psSphereRot_TEOtoCEO:          return NULL for NULL time input. (z)");
     356        psFree(earthRot);
     357        psFree(result);
     358        psFree(cube);
     359        psFree(obj);
     360        skip_end();
     361        psFree(polarTideCorr);
     362    }
     363
     364    //Tests for psSphereRot_ITRStoTEO
     365    //Return NULL for NULL earthPole input
     366    {
     367        psFree(out);
     368        out = psSphereRot_ITRStoTEO(NULL);
     369        ok( out == NULL,
     370            "psSphereRot_ITRStoTEO:         return NULL for NULL earthPole input.");
     371        in = psEarthPoleAlloc();
     372    }
     373    //Return correct psSphereRot for valid input
     374    {
     375        q0 = -1.0567571848664005e-6;
     376        q1 = 3.218036562931509e-7;
     377        q2 = -3.3580195807204483e-12;
     378        q3 = -0.9999999999993899;
     379        in->x = SEC_TO_RAD(-0.13275353774074533);
     380        in->y = SEC_TO_RAD(0.4359436319739848);
     381        in->s = SEC_TO_RAD(-4.2376965863576153e-10);
     382        out = psSphereRot_ITRStoTEO(in);
     383        skip_start(  out == NULL, 4,
     384                     "Skipping 4 tests because psSphereRot output is NULL!");
     385        ok_double_tol( out->q0, q0, 0.0001,
     386                       "psSphereRot_ITRStoTEO:         return correct psSphereRot for valid input (q0).");
     387        ok_double_tol( out->q1, q1, 0.0001,
     388                       "psSphereRot_ITRStoTEO:         return correct psSphereRot for valid input (q1).");
     389        ok_double_tol( out->q2, q2, 0.0001,
     390                       "psSphereRot_ITRStoTEO:         return correct psSphereRot for valid input (q2).");
     391        ok_double_tol( out->q3, q3, 0.0001,
     392                       "psSphereRot_ITRStoTEO:         return correct psSphereRot for valid input (q3).");
     393        skip_end();
     394    }
     395
     396    //Check for Memory leaks
     397    {
     398        psFree(temp);
     399        psFree(time2);
     400        psFree(time);
     401        psFree(out);
     402        psFree(in);
     403        checkMem();
     404    }
     405}
     406
     407#define SPHERE_PRECESS_TP1_R            0.0               //    0.0       degrees
     408#define SPHERE_PRECESS_TP1_D            0.0               //    0.0       degrees
     409#define SPHERE_PRECESS_TP1_EXPECT_R     6.238453          //  357.437     degrees
     410#define SPHERE_PRECESS_TP1_EXPECT_D    -0.019426          //   -1.113     degrees
     411#define SPHERE_PRECESS_TP2_R            0.0               //    0.0       degrees
     412#define SPHERE_PRECESS_TP2_D            1.570796          //   90.0       degrees
     413#define SPHERE_PRECESS_TP2_EXPECT_R     6.260828          //  358.719     degrees
     414#define SPHERE_PRECESS_TP2_EXPECT_D     1.551353          //   88.886     degrees
     415#define SPHERE_PRECESS_TP3_R            3.141593          //  180.0       degrees
     416#define SPHERE_PRECESS_TP3_D            0.523599          //   30.0       degrees
     417#define SPHERE_PRECESS_TP3_EXPECT_R     3.096616          //  177.423     degrees
     418#define SPHERE_PRECESS_TP3_EXPECT_D     0.543024          //   31.113     degrees
     419#define ERROR_TOL   0.0001
     420
     421void testSpherePrecess(void)
     422{
     423    diag("  >>>Test 4:  psSpherePrecess");
     424    psSphereRot *rot = NULL;
     425    psTime *fromTime = NULL;
     426    psTime *toTime = NULL;
     427    psSphere* inputCoord  = psSphereAlloc();
     428    psSphere* outputCoord = NULL;
     429
     430    //Return NULL for NULL time inputs
     431    {
     432        rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
     433        ok( rot == NULL,
     434            "psSpherePrecess:               return NULL for NULL time inputs.");
     435    }
     436    //Return NULL for invalid mode input
     437    {
     438        fromTime    = psTimeFromMJD(MJD_2100);
     439        toTime      = psTimeFromMJD(MJD_1900);
     440        rot = psSpherePrecess(fromTime, toTime, -1);
     441        ok( rot == NULL,
     442            "psSpherePrecess:               return NULL for invalid mode input.");
     443    }
     444    //Return correct psSphereRot for valid inputs, mode = ROUGH
     445    {
     446        inputCoord->r = SPHERE_PRECESS_TP1_R;
     447        inputCoord->d = SPHERE_PRECESS_TP1_D;
     448        inputCoord->rErr = 0.0;
     449        inputCoord->dErr = 0.0;
     450
     451        rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
     452        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     453        if (outputCoord->r < -0.0001)
     454            outputCoord->r += 2.0 * M_PI;
     455        skip_start( rot == NULL || outputCoord == NULL, 2,
     456                    "Skipping 2 tests because psSphereRot output is NULL!");
     457        ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP1_EXPECT_R, ERROR_TOL,
     458                       "psSpherePrecess:               return correct psSphereRot for valid"
     459                       " inputs and PS_PRECESS_ROUGH mode. (r)");
     460        ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP1_EXPECT_D, ERROR_TOL,
     461                       "psSpherePrecess:               return correct psSphereRot for valid"
     462                       " inputs and PS_PRECESS_ROUGH mode. (d)");
     463        skip_end();
     464        psFree(outputCoord);
     465        psFree(rot);
     466    }
     467    //Return correct psSphereRot for valid inputs, mode = ROUGH
     468    {
     469        inputCoord->r = SPHERE_PRECESS_TP2_R;
     470        inputCoord->d = SPHERE_PRECESS_TP2_D;
     471        inputCoord->rErr = 0.0;
     472        inputCoord->dErr = 0.0;
     473
     474        rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
     475        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     476        if (outputCoord->r < -0.0001)
     477            outputCoord->r += 2.0 * M_PI;
     478        skip_start( rot == NULL || outputCoord == NULL, 2,
     479                    "Skipping 2 tests because psSphereRot output is NULL!");
     480        ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP2_EXPECT_R, ERROR_TOL,
     481                       "psSpherePrecess:               return correct psSphereRot for valid"
     482                       " inputs and PS_PRECESS_ROUGH mode. (r)");
     483        ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP2_EXPECT_D, ERROR_TOL,
     484                       "psSpherePrecess:               return correct psSphereRot for valid"
     485                       " inputs and PS_PRECESS_ROUGH mode. (d)");
     486        skip_end();
     487        psFree(outputCoord);
     488        psFree(rot);
     489    }
     490    //Return correct psSphereRot for valid inputs, mode = ROUGH
     491    {
     492        inputCoord->r = SPHERE_PRECESS_TP3_R;
     493        inputCoord->d = SPHERE_PRECESS_TP3_D;
     494        inputCoord->rErr = 0.0;
     495        inputCoord->dErr = 0.0;
     496
     497        rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
     498        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     499        if (outputCoord->r < -0.0001)
     500            outputCoord->r += 2.0 * M_PI;
     501        skip_start( rot == NULL || outputCoord == NULL, 2,
     502                    "Skipping 2 tests because psSphereRot output is NULL!");
     503        ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP3_EXPECT_R, ERROR_TOL,
     504                       "psSpherePrecess:               return correct psSphereRot for valid"
     505                       " inputs and PS_PRECESS_ROUGH mode. (r)");
     506        ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP3_EXPECT_D, ERROR_TOL,
     507                       "psSpherePrecess:               return correct psSphereRot for valid"
     508                       " inputs and PS_PRECESS_ROUGH mode. (d)");
     509        skip_end();
     510        psFree(outputCoord);
     511        psFree(rot);
     512    }
     513    //Return correct psSphereRot for valid inputs, mode = COMPLETE_A
     514    {
     515        rot = psSpherePrecess(NULL, toTime, PS_PRECESS_COMPLETE_A);
     516        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     517        ok( outputCoord != NULL && rot != NULL,
     518            "psSpherePrecess:               return correct psSphereRot for valid"
     519            " inputs and PS_PRECESS_COMPLETE_A mode.");
     520        psFree(outputCoord);
     521        psFree(rot);
     522    }
     523    //Return correct psSphereRot for valid inputs, mode = COMPLETE_B
     524    {
     525        rot = psSpherePrecess(NULL, toTime, PS_PRECESS_COMPLETE_B);
     526        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     527        ok( outputCoord != NULL && rot != NULL,
     528            "psSpherePrecess:               return correct psSphereRot for valid"
     529            " inputs and PS_PRECESS_COMPLETE_B mode.");
     530        psFree(outputCoord);
     531        psFree(rot);
     532    }
     533    //Return correct psSphereRot for valid inputs, mode = IAU2000A
     534    {
     535        rot = psSpherePrecess(NULL, toTime, PS_PRECESS_IAU2000A);
     536        outputCoord = psSphereRotApply(NULL, rot, inputCoord);
     537        ok( outputCoord != NULL && rot != NULL,
     538            "psSpherePrecess:               return correct psSphereRot for valid"
     539            " inputs and PS_PRECESS_IAU2000A mode.");
     540    }
     541
     542    //Check for Memory leaks
     543    {
     544        psFree(fromTime);
     545        psFree(toTime);
     546        psFree(outputCoord);
     547        psFree(rot);
     548        checkMem();
     549    }
     550
     551}
     552
     553
Note: See TracChangeset for help on using the changeset viewer.