Index: /trunk/psLib/src/astro/Makefile.am =================================================================== --- /trunk/psLib/src/astro/Makefile.am (revision 9975) +++ /trunk/psLib/src/astro/Makefile.am (revision 9976) @@ -1,9 +1,9 @@ #Makefile for astronomy functions of psLib # -AM_CFLAGS=$(CFLAGS) -DPS_CONFIG_FILE_DEFAULT=\"$(sysconfdir)/pslib/pslib.config\" +AM_CFLAGS=$(CFLAGS) -DPS_CONFIG_FILE_DEFAULT=\"$(sysconfdir)/pslib.config\" noinst_LTLIBRARIES = libpslibastro.la -libpslibastro_la_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS) +libpslibastro_la_CPPFLAGS = $(SRCINC) $(PSLIB_CFLAGS) libpslibastro_la_SOURCES = \ psTime.c \ Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 9975) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 9976) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-20 00:17:21 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-11-14 21:15:43 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -284,4 +284,9 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Aberration speed should not be equal to 0.\n"); + return NULL; + } + if (fabs(speed) > 1.0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Aberration speed should not greater than the speed of light!.\n"); return NULL; } @@ -445,5 +450,5 @@ // Convert psTime to MJD double MJD = psTimeToMJD(time); - if (MJD == NAN) { + if (isnan(MJD)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to MJD failed. Invalid input time.\n"); @@ -594,5 +599,5 @@ //Convert the input time to MJD. If NAN is returned, return NULL for the function. double MJD = psTimeToMJD(time); - if (MJD == NAN) { + if (isnan(MJD)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to MJD failed. Invalid input time.\n"); @@ -762,5 +767,5 @@ //Check if tidal corrections should be included. //If so, make sure values are positive and in the correct range. - if (tidalCorr != NULL && tidalCorr->s != 0.0) { + if (tidalCorr != NULL && fabs(tidalCorr->s) > FLT_EPSILON) { int nsec = in->nsec + (int)(tidalCorr->s * 1e9); if (nsec > 1e9) { @@ -777,5 +782,5 @@ //Calculate the Julian Date from the input time in UT1 format. double T = psTimeToJD(in); - if (T == NAN) { + if (isnan(T)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to JD failed. Invalid input time.\n"); @@ -803,5 +808,5 @@ double MJD = psTimeToMJD(time); - if (MJD == NAN) { + if ( isnan(MJD) ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to MJD failed. Invalid input time.\n"); @@ -858,7 +863,7 @@ k = 2; } - if (k > (numRows-2)) { - k = numRows-2; - } + // if (k > (numRows-2)) { + // k = numRows-2; + // } for (int m = k-1; m <= k+2; m++) { xTerm = X->data.F64[m]; @@ -935,5 +940,5 @@ // Convert psTime to MJD double MJD = psTimeToMJD(time); - if (MJD == NAN) { + if (isnan(MJD)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to MJD failed. Invalid input time.\n"); @@ -1019,5 +1024,5 @@ double MJD = psTimeToMJD(time); // printf("\nMJD check = %.13g\n", MJD); - if (MJD == NAN) { + if (isnan(MJD)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Time conversion to MJD failed. Invalid input time.\n"); @@ -1252,16 +1257,20 @@ // Calculate Julian centuries //If either input time is NULL, assume it to be J2000 -> from SDRS as of rev 18 + psTime *from = NULL; + psTime *to = NULL; if (fromTime == NULL) { fromMJD = MJD_2000; + from = psTimeFromMJD(fromMJD); } else { fromMJD = psTimeToMJD(fromTime); + from = psTimeCopy(fromTime); } if (toTime == NULL) { toMJD = MJD_2000; + to = psTimeFromMJD(toMJD); } else { toMJD = psTimeToMJD(toTime); - } - psTime *from = NULL; - psTime *to = NULL; + to = psTimeCopy(toTime); + } if (fromMJD > toMJD) { psWarning("From time is later than to time in psSpherePrecess.\n"); @@ -1290,4 +1299,6 @@ // Create transform with proper constants psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP); + psFree(from); + psFree(to); return tmpST; } else if (mode == PS_PRECESS_IAU2000A) { @@ -1296,16 +1307,4 @@ //rotation corresponding to fromTime to the toTime rotation. - //Since the time inputs are allowed to be NULL, either convert the MJD time - //or copy to non-NULL time - if (fromTime == NULL) { - from = psTimeFromMJD(fromMJD); - } else { - from = psTimeCopy(fromTime); - } - if (toTime == NULL) { - to = psTimeFromMJD(toMJD); - } else { - to = psTimeCopy(toTime); - } //Calculate the earthpoles and quaternions corresponding to each time (from, to). //Combine the quaternions to produce the output psSphereRot. @@ -1329,16 +1328,4 @@ //for COMPLETE_A come from the IERS Bulletin A. - //Since the time inputs are allowed to be NULL, either convert the MJD time - //or copy to non-NULL time - if (fromTime == NULL) { - from = psTimeFromMJD(fromMJD); - } else { - from = psTimeCopy(fromTime); - } - if (toTime == NULL) { - to = psTimeFromMJD(toMJD); - } else { - to = psTimeCopy(toTime); - } //Calculate the earthpoles and quaternions corresponding to each time (from, to). //Add in the precession corrections from IERS bulletin A. @@ -1371,17 +1358,4 @@ //for COMPLETE_B come from the IERS Bulletin B. - //Since the time inputs are allowed to be NULL, either convert the MJD time - //or copy to non-NULL time - if (fromTime == NULL) { - from = psTimeFromMJD(fromMJD); - } else { - from = psTimeCopy(fromTime); - } - if (toTime == NULL) { - to = psTimeFromMJD(toMJD); - } else { - to = psTimeCopy(toTime); - } - //Calculate the earthpoles and quaternions corresponding to each time (from, to). //Add in the precession corrections from IERS bulletin B. Index: /trunk/psLib/test/astro/tap_psEarthOrientation_corrections.c =================================================================== --- /trunk/psLib/test/astro/tap_psEarthOrientation_corrections.c (revision 9975) +++ /trunk/psLib/test/astro/tap_psEarthOrientation_corrections.c (revision 9976) @@ -18,43 +18,261 @@ #include "pstap.h" -static void testSphereRotCreate(void); -//static void testSphereRotConvert(void); -//static void testSphereOffsets(void); +static void testEOCInit(void); +static void testAberration(void); +static void testGravDef(void); +static void testEOC_Corrs(void); + +#define MJD_1900 15021.0 // Modified Julian Day 1/1/1900 00:00:00 +#define MJD_2100 88069.0 // Modified Julian Day 1/1/2100 00:00:00 int main(void) { - plan_tests(1); + plan_tests(8); diag("Tests for psEarthOrientation Correction Functions"); - testSphereRotCreate(); - // testSphereRotConvert(); + // Initialize library internal structures + psLibInit("pslib.config"); + + testEOCInit(); + testAberration(); + testGravDef(); + testEOC_Corrs(); // testSphereOffsets(); + + // Cleanup library + psLibFinalize(); done(); } -void testSphereRotCreate(void) +void testEOCInit(void) { - diag(" >>>Test 1: psSphereRot Creation Functions"); + diag(" >>>Test 1: p_psEOCInit "); + //Test for psEarthPoleAlloc + //Return properly allocated psEarthPole + /* { + skip_start( ep == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ep = psEarthPoleAlloc(); + ok( ep != NULL, + "psEarthPoleAlloc: return properly allocated psEarthPole."); + psFree(ep); + ep = NULL; + } + */ + + + //Check for Memory leaks + { + checkMem(); + } + +} + +void testAberration(void) +{ + diag(" >>>Test 2: psAberration "); + psSphere *apparent = NULL; + psSphere *empty = NULL; + psCube *actualCube = psCubeAlloc(); + //values from After gravity deflection// + actualCube->x = -0.35961949760293604; + actualCube->y = 0.5555613950298085; + actualCube->z = 0.7496835020836093; + psSphere *actual = psCubeToSphere(actualCube); + psCube *cubeDir = psCubeAlloc(); + cubeDir->x = 5148.713262821658; + cubeDir->y = -26945.04752348012; + cubeDir->z = -11682.787302030947; + cubeDir->x += -357.6031690489248; + cubeDir->y += 248.46429758174693; + cubeDir->z += 0.09694774143797581; + psSphere *direction = psCubeToSphere(cubeDir); + double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z); + double c = 299792458.0; // Speed of light in vacuum (src:NIST) m/s + speed /= c; + + //Tests for psAberration + //Return NULL for NULL actual sphere input + { + empty = psAberration(apparent, NULL, direction, speed); + ok( empty == NULL, + "psAberration: return NULL for NULL actual sphere input."); + } + //Return NULL for NULL direction sphere input + { + empty = psAberration(apparent, actual, NULL, speed); + ok( empty == NULL, + "psAberration: return NULL for NULL direction sphere input."); + } + //Return NULL for speed = 0.0 + { + empty = psAberration(apparent, actual, direction, 0.0); + ok( empty == NULL, + "psAberration: return NULL for zero speed input."); + } + //Return NULL for speed > 1.0 + { + empty = psAberration(apparent, actual, direction, 1.1); + ok( empty == NULL, + "psAberration: return NULL for impossible speed input (> c)."); + } + + //Return correct values for valid inputs + { + double x, y, z; + x = -0.35963388069046304; + y = 0.5555192509816625; + z = 0.7497078321908413; + apparent = psAberration(apparent, actual, direction, speed); + skip_start( apparent == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + psCube *outCube = psSphereToCube(apparent); + ok_double_tol( outCube->x, x, 0.001, + "psAberration: return correct sphere for valid inputs."); + ok_double_tol( outCube->y, y, 0.001, + "psAberration: return correct sphere for valid inputs."); + ok_double_tol( outCube->z, z, 0.001, + "psAberration: return correct sphere for valid inputs."); + psFree(outCube); + skip_end(); + } + + //Check for Memory leaks + { + psFree(actualCube); + psFree(cubeDir); + psFree(apparent); + psFree(actual); + psFree(direction); + checkMem(); + } + +} + +void testGravDef(void) +{ + diag(" >>>Test 3: psGravityDeflection "); + //Test for psGravityDeflection + //Return properly allocated psEarthPole /* - //Tests for psSphereRotAlloc - //Return NULL for NAN input - { - psSphereRot *s1 = psSphereRotAlloc(ALPHA_P, NAN, PHI_P); - ok( s1 == NULL, - "psSphereRotAlloc: return NULL for NAN input."); - } - //Tests for psMemCheckSphereRot - //Make sure psMemCheckSphereRot works correctly - return false - { - int j = 2; - ok( !psMemCheckSphereRot(&j), - "psMemCheckSphereRot: return false for non-SphereRot input."); - skip_start( s1 == NULL, 1, - "Skipping 1 tests because psSphereRot is NULL!"); - skip_end(); - } + psSphere *apparent = NULL; + psSphere *empty = NULL; + psCube *actualCube = psCubeAlloc(); + actualCube->x = -0.3596195125758298; + actualCube->y = 0.5555613903455866; + actualCube->z = 0.7496834983724809; + psSphere *actual = psCubeToSphere(actualCube); + psCube *sunCube = psCubeAlloc(); + sunCube->x = 1.467797790127511e11; + sunCube->y = 2.5880956908748722e10; + sunCube->z = 1.1220046291457653e10; + double sunLength = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z); + sunCube->x /= sunLength; + sunCube->y /= sunLength; + sunCube->z /= sunLength; + if (VERBOSE) { + printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n", + sunCube->x, sunCube->y, sunCube->z); + } + psSphere *sun = psCubeToSphere(sunCube); + if (VERBOSE) { + printf("sunSphere = r, d = %.13g, %.13g\n", sun->r, sun->d); + } + psCube *outCube = psCubeAlloc(); + + empty = psGravityDeflection(apparent, empty, sun); + if (empty != NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "psGravityDeflection Failed to return NULL for NULL actual input sphere.\n"); + return 1; + } + empty = psGravityDeflection(apparent, actual, empty); + if (empty != NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "psGravityDeflection Failed to return NULL for NULL sun input sphere.\n"); + return 2; + } + + apparent = psGravityDeflection(NULL, actual, sun); + // apparent->r *= -1.0; + // apparent->d *= -1.0; + psSphere *result2; + // psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + // psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + // psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN); + if (VERBOSE) { + printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g -- \n", + actualCube->x, actualCube->y, actualCube->z); + } + // psCube *outCube = psSphereToCube(result); + // printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", + // outCube->x, outCube->y, outCube->z); + psCube *outCube2 = psSphereToCube(apparent); + if (VERBOSE) { + printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g -- \n", + outCube2->x, outCube2->y, outCube2->z); + } + double x, y, z; + x = -0.35961949760293604; + y = 0.5555613950298085; + z = 0.7496835020836093; + + if (VERBOSE) { + printf(" -- expectCube = x,y,z = %.13g, %.13g, %.13g -- \n\n", x, y, z); + + // if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON || + // fabs(z - outCube->z) > DBL_EPSILON ) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // "psGravityDeflection returned incorrect values.\n"); + printf("expect-actual= x,y,z = %.13g, %.13g, %.13g \n", + (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) ); + printf("expect-result= x,y,z = %.13g, %.13g, %.13g \n", + (x - outCube2->x), (y - outCube2->y), (z - outCube2->z) ); + printf("result-actual= x,y,z = %.13g, %.13g, %.13g \n", + (outCube2->x - actualCube->x), (outCube2->y - actualCube->y), + (outCube2->z - actualCube->z) ); + } + // return 1; + // } + // psFree(result2); + outCube2->x = x; + outCube2->y = y; + outCube2->z = z; + result2 = psCubeToSphere(outCube2); + // psFree(result); + psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN); + psFree(result2); + result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + if (VERBOSE) { + printf("The apparent output sphere = r,d = %.13g, %.13g\n", result2->r, result2->d); + printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d); + } + psFree(result2); + + psFree(outCube); + psFree(outCube2); + psFree(sunCube); + psFree(actualCube); + psFree(result); + psFree(actual); + psFree(apparent); + psFree(sun); + + + + { + skip_start( ep == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ep = psEarthPoleAlloc(); + ok( ep != NULL, + "psEarthPoleAlloc: return properly allocated psEarthPole."); + psFree(ep); + ep = NULL; + } */ + + //Check for Memory leaks { @@ -64,4 +282,283 @@ } - - +void testEOC_Corrs(void) +{ + diag(" >>>Test 4: psEOC Correction Functions"); + //Tests for psEOC_PrecessionCorr + /* + psTime *empty = NULL; + psTime *time2 = psTimeAlloc(PS_TIME_UTC); + time2->sec = timesec; + time2->nsec = 0; + time2->leapsecond = false; + // time2 = psTimeConvert(time2, PS_TIME_TAI); + + //Tests for Precession Correction function// + //Return NULL for NULL time input + psEarthPole *pcorr = NULL; + psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message"); + pcorr = psEOC_PrecessionCorr(empty, PS_IERS_A); + if (pcorr != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_PrecessionCorr failed to return NULL for NULL time input.\n"); + return 5; + } + + //Return NULL for Invalid IERS table + psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message"); + pcorr = psEOC_PrecessionCorr(time2, 3); + if (pcorr != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_PrecessionCorr failed to return NULL for incorrect IERS table.\n"); + return 6; + } + psFree(pcorr); + + //Check values from IERS table A + pcorr = psEOC_PrecessionCorr(time2, PS_IERS_A); + if ( pcorr == NULL ) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_PrecessionCorr returned NULL for valid inputs.\n"); + return 7; + } else { + if (VERBOSE) { + printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", + pcorr->x, pcorr->y, pcorr->s); + } + } + psFree(pcorr); + + //Check values from IERS table B + pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B); + if ( pcorr == NULL ) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_PrecessionCorr returned NULL for valid inputs.\n"); + return 9; + } else { + double xx, yy, ss; + xx = 0.06295703125; + yy = -0.0287618408203125; + ss = 0.0; + xx = SEC_TO_RAD(xx) * 1e-3; + yy = SEC_TO_RAD(yy) * 1e-3; + // if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON + // || fabs(pcorr->s - ss) > DBL_EPSILON) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // " psEOC_PrecessionCorr return incorrect values.\n"); + if (VERBOSE) { + printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + pcorr->x, pcorr->y, pcorr->s); + printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n", xx, yy, ss); + printf(" A difference of: %.13g, %.13g, %.13g\n\n", + (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) ); + } + // return 10; + // } + } + //precess is the *actual* output from PrecessionModel + PrecessionCorr + psEarthPole *precess = psEOC_PrecessionModel(time2); + precess->x += pcorr->x; + precess->y += pcorr->y; + double xCorr, yCorr; + xCorr = 3.05224300720406e-10; + yCorr = -1.39441339235822e-10; + //pcorr is the *expected* output from PrecessionModel// + pcorr->x = 2.857175590089105e-4; + pcorr->y = 2.3968739377734732e-5; + pcorr->s = -1.3970066457904322e-8; + pcorr->x += xCorr; + pcorr->y += yCorr; + psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess); + psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv); + double q0, q1, q2; + q0 = -1.1984522406756289e-5; + q1 = 1.4285893358610674e-4; + q2 = 1.2191193518914336e-10; + psSphereRot *pni = psSphereRot_CEOtoGCRS(pcorr); + if (fabs(pni->q0-q0) > FLT_EPSILON || fabs(pni->q1-q1) > FLT_EPSILON || + fabs(pni->q2-q2) > FLT_EPSILON ) { + printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n"); + } + // printf(" Output from CEOtoGCRS only = %.13g,%.13g,%.13g,%.13g\n", + // pni->q0, pni->q1, pni->q2, pni->q3); + if (VERBOSE) { + printf(" Expected sphere rotation = %.13g, %.13g, %.13g\n", q0,q1,q2); + printf(" Result sphere rotation = %.13g, %.13g, %.13g\n", + precessNutInv->q0, precessNutInv->q1, precessNutInv->q2); + printf(" Difference = %.13g, %.13g, %.13g\n\n", + precessNutInv->q0-q0, precessNutInv->q1-q1, precessNutInv->q2-q2); + } + psCube *objC = psCubeAlloc(); + // objC->x = -3.5963388069046304; + // objC->y = 0.5555192509816625; + // objC->z = 0.7497078321908413; + // objSetup(); + + //This is the sphere rotation for the *expected* precession output// + psSphereRot *pn = psSphereRotConjugate(NULL, pni); + + // psSphere *sphere = psSphereAlloc(); + // *sphere = *obj; + // psFree(obj); + + //create a psSphere for (from) the start position given in eoc_testing// + objC->x = -0.35963388069046304; + objC->y = 0.5555192509816625; + objC->z = 0.7497078321908413; + psSphere *sphere = psCubeToSphere(objC); + + psSphere *expect = psSphereRotApply(NULL, pn, sphere); + //expected results below - stored in: sphere // + double x,y,z; + x = -0.3598480726985338; + y = 0.5555012823608123; + z = 0.7496183628158023; + if (VERBOSE) { + printf("\n<x, objC->y, objC->z); + //printf(" Difference = %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z); + // x = objC->x; + // y = objC->y; + // z = objC->z; + psSphere *result = psSphereRotApply(NULL, precessNut, sphere); + psFree(objC); + objC = psSphereToCube(result); + double xx = greatCircle(result, expect); + if (VERBOSE) { + printf("<x, objC->y, objC->z); + printf(" Difference = %.13g, %.13g, %.13g\n\n", + objC->x-x, objC->y-y, objC->z-z); + printf("GREAT CIRCLE DIFFERENCE = %.13g \n", xx); + } + + psFree(precess); + psFree(precessNut); + psFree(precessNutInv); + psFree(expect); + psFree(objC); + + psFree(sphere); + psFree(result); + psFree(pn); + psFree(pni); + psFree(pcorr); + psFree(time2); + if (!p_psEOCFinalize() ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); + return 12; + } + + */ + + //Tests for psEOC_PolarTideCorr + /* + psTime *in = psTimeAlloc(PS_TIME_UTC); + in->sec = timesec; + in->nsec = 0; + in->leapsecond = false; + psTime *empty = NULL; + psEarthPole *eop = NULL; + + //Return NULL for NULL input time + psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message"); + eop = psEOC_PolarTideCorr(empty); + if (eop != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_PolarTideCorr failed to return NULL for NULL input time.\n"); + return 1; + } + + eop = psEOC_PolarTideCorr(in); + if (eop == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_PolarTideCorr returned NULL for valid input time.\n"); + return 2; + } else { + if (VERBOSE) { + printf("\nPolarTideCorr output = x,y,s = %.13g, %.13g, %.13g\n", + eop->x, eop->y, eop->s); + } + } + + psFree(in); + psFree(eop); + + */ + + + + //Tests for psEOC_NutationCorr + /* + psTime *in = psTimeAlloc(PS_TIME_UTC); + in->sec = timesec; + in->nsec = 0; + in->leapsecond = false; + psTime *empty = NULL; + psEarthPole *nute = NULL; + + //Return NULL for NULL input time. + psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message"); + nute = psEOC_NutationCorr(empty); + if (nute != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_NutationCorr failed to return NULL for NULL input time.\n"); + return 1; + } + //Return NULL for UT1 time input + *//* psTime *UT1time = psTimeAlloc(PS_TIME_UT1); + nute = psEOC_NutationCorr(UT1time); + if (nute != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_NutationCorr failed to return NULL for UT1 input time.\n"); + return 2; + } + psFree(UT1time); + *//* + //Check return values from valid nutation time input + nute = psEOC_NutationCorr(in); + if ( nute == NULL ) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_NutationCorr returned NULL for valid input.\n"); + return 3; + } else { + if (VERBOSE) { + printf("Nutation Correction output = x,y,s = %.13g, %.13g, %.13g\n\n", + nute->x, nute->y, nute->s); + } + } + psFree(nute); + psFree(in); + + if (!p_psEOCFinalize() ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); + return 12; + } + + */ + //Return properly allocated psEarthPole + /* + + + { + skip_start( ep == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ep = psEarthPoleAlloc(); + ok( ep != NULL, + "psEarthPoleAlloc: return properly allocated psEarthPole."); + psFree(ep); + ep = NULL; + } + */ + + + //Check for Memory leaks + { + checkMem(); + } + +} + + Index: /trunk/psLib/test/astro/tap_psEarthOrientation_motion.c =================================================================== --- /trunk/psLib/test/astro/tap_psEarthOrientation_motion.c (revision 9975) +++ /trunk/psLib/test/astro/tap_psEarthOrientation_motion.c (revision 9976) @@ -4,5 +4,5 @@ * Description: Tests for psEarthPoleAlloc, psEOC_PrecessionModel, psSpherePrecess * psSphereRot_CEOtoGCRS, psSphereRot_TEOtoCEO, - * psEOC_GetPolarMotion, psSphereRot_ITRStoTEO, + * psEOC_GetPolarMotion, psSphereRot_ITRStoTEO * * Author: dRob , (C) 2006 @@ -18,32 +18,106 @@ #include "pstap.h" -static void testSphereRotCreate(void); -//static void testSphereRotConvert(void); -//static void testSphereOffsets(void); +static void testPrecessionModel(void); +static void testPolarMotion(void); +static void testSphereRots(void); +static void testSpherePrecess(void); + +#define MJD_1900 15021.0 // Modified Julian Day 1/1/1900 00:00:00 +#define MJD_2100 88069.0 // Modified Julian Day 1/1/2100 00:00:00 int main(void) { - plan_tests(1); + plan_tests(48); diag("Tests for psEarthOrientation Motion Functions"); - testSphereRotCreate(); - // testSphereRotConvert(); - // testSphereOffsets(); + // Initialize library internal structures + psLibInit("pslib.config"); + + testPrecessionModel(); + testPolarMotion(); + testSphereRots(); + testSpherePrecess(); + + // Cleanup library + psLibFinalize(); done(); } -void testSphereRotCreate(void) +void testPrecessionModel(void) { - diag(" >>>Test 1: psSphereRot Creation Functions"); + diag(" >>>Test 1: psEOC_PrecessionModel"); + + psEarthPole *ep = NULL; + psTime *testTime = NULL; + + //Test for psEarthPoleAlloc + //Return properly allocated psEarthPole + { + ep = psEarthPoleAlloc(); + ok( ep != NULL, + "psEarthPoleAlloc: return properly allocated psEarthPole."); + psFree(ep); + ep = NULL; + } + //Tests for psEOC_PrecessionModel + //Return NULL for NULL time input. + { + ep = psEOC_PrecessionModel(NULL); + ok( ep == NULL, + "psEOC_PrecessionModel: return NULL for NULL time input."); + } + //Return NULL for UT1 time input. + { + testTime = psTimeAlloc(PS_TIME_UT1); + ep = psEOC_PrecessionModel(testTime); + ok( ep == NULL, + "psEOC_PrecessionModel: return NULL for UT1 time input."); + } + //Return NULL for invalid time input. + { + psFree(testTime); + testTime = psTimeAlloc(PS_TIME_UTC); + testTime->nsec = 2e9; + ep = psEOC_PrecessionModel(testTime); + ok( ep == NULL, + "psEOC_PrecessionModel: return NULL for invalid time input."); + testTime->nsec = 0; + } + //Return NULL for failed eoc init - missing file + // printf("\n filename = '%s' \n", p_psGetConfigFileName()); + { + testTime->sec = 1049160600; + testTime->nsec = 0; + testTime->leapsecond = false; + skip_start(rename("../../etc/pslib.config", "../../etc/pslib_config.bak"), + 1, "Skipping 1 tests because file rename failed!"); + ep = psEOC_PrecessionModel(testTime); + ok( ep == NULL, + "psEOC_PrecessionModel: return NULL for failed eoc init."); + rename("../../etc/pslib_config.bak", "../../etc/pslib.config"); + skip_end(); + } + //Return valid EarthPole for valid time input + { + double x, y, s; + x = 2.857175590089105e-4; + y = 2.3968739377734732e-5; + s = -1.3970066457904322e-8; + ep = psEOC_PrecessionModel(testTime); + skip_start( ep == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ok_double_tol(ep->x, x, 0.1, + "psEOC_PrecessionModel: return valid EarthPole for valid inputs (x)."); + ok_double_tol(ep->y, y, 0.1, + "psEOC_PrecessionModel: return valid EarthPole for valid inputs (y)."); + ok_double_tol(ep->s, s, 0.1, + "psEOC_PrecessionModel: return valid EarthPole for valid inputs (s)."); + skip_end(); + } + + /* - //Tests for psSphereRotAlloc - //Return NULL for NAN input - { - psSphereRot *s1 = psSphereRotAlloc(ALPHA_P, NAN, PHI_P); - ok( s1 == NULL, - "psSphereRotAlloc: return NULL for NAN input."); - } //Tests for psMemCheckSphereRot //Make sure psMemCheckSphereRot works correctly - return false @@ -59,4 +133,6 @@ //Check for Memory leaks { + psFree(testTime); + psFree(ep); checkMem(); } @@ -64,2 +140,414 @@ } +void testPolarMotion(void) +{ + diag(" >>>Test 2: psEOC_GetPolarMotion "); + + psTime *in = psTimeAlloc(PS_TIME_UTC); + in->sec = 1049160600; + in->nsec = 0; + in->leapsecond = false; + psEarthPole *polarMotion = NULL; + + //Tests for psEOC_GetPolarMotion + //Return NULL for NULL time input. + { + polarMotion = psEOC_GetPolarMotion(NULL, PS_IERS_B); + ok( polarMotion == NULL, + "psEOC_GetPolarMotion: return NULL for NULL time input."); + } + //Return NULL for invalid IERS table input + { + polarMotion = psEOC_GetPolarMotion(NULL, PS_IERS_B+1); + ok( polarMotion == NULL, + "psEOC_GetPolarMotion: return NULL for invalid IERS table input."); + } + //Return NULL for invalid time input. + { + in->nsec = 2e9; + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); + ok( polarMotion == NULL, + "psEOC_GetPolarMotion: return NULL for invalid time input."); + in->nsec = 0; + } + //Return NULL for failed eoc init - missing file + { + p_psEOCFinalize(); + skip_start(rename("../../etc/pslib.config", "../../etc/pslib_config.bak"), + 1, "Skipping 1 tests because file rename failed!"); + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); + ok( polarMotion == NULL, + "psEOC_GetPolarMotion: return NULL for failed eoc init."); + rename("../../etc/pslib_config.bak", "../../etc/pslib.config"); + skip_end(); + } + //Return valid EarthPole for valid time input - IERS B + { + double x, y, s; + x = -6.454389659777e-07; + y = 2.112606414597e-06; + s = 0.0; + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); + skip_start( polarMotion == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ok_double_tol(polarMotion->x, x, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(x) - IERS B."); + ok_double_tol(polarMotion->y, y, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(y) - IERS B."); + ok_double_tol(polarMotion->s, s, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(s) - IERS B."); + skip_end(); + psFree(polarMotion); + } + //Return valid EarthPole for valid time input - IERS A + { + double x, y, s; + x = -6.45381397904e-07; + y = 2.112819726698e-06; + s = 0.0; + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A); + skip_start( polarMotion == NULL, 3, + "Skipping 3 tests because psEarthPole is NULL!"); + ok_double_tol(polarMotion->x, x, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(x) - IERS A."); + ok_double_tol(polarMotion->y, y, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(y) - IERS A."); + ok_double_tol(polarMotion->s, s, 0.1, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs " + "(s) - IERS A."); + skip_end(); + psFree(polarMotion); + } + //Return valid EarthPole for valid time input - IERS A + { + psTime *firstTime = psTimeFromMJD(41684.50); + polarMotion = psEOC_GetPolarMotion(firstTime, PS_IERS_B); + ok( polarMotion != NULL, + "psEOC_GetPolarMotion: return valid EarthPole for valid inputs."); + psFree(firstTime); + } + + + //Check for Memory leaks + { + psFree(in); + psFree(polarMotion); + checkMem(); + } + +} + +static psSphere *objSetup(void) +{ + psSphere *obj = NULL; + psCube *tempCube = psCubeAlloc(); + tempCube->x = -0.3598480726985338; + tempCube->y = 0.5555012823608123; + tempCube->z = 0.7496183628158023; + obj = psCubeToSphere(tempCube); + psFree(tempCube); + return obj; +} + +void testSphereRots(void) +{ + diag(" >>>Test 3: psSphereRot Functions"); + psSphereRot *out = NULL; + psEarthPole *in = NULL; + double q0,q1,q2,q3; + q0 = -1.1984522406756289e-5; + q1 = 1.4285893358610674e-4; + q2 = 1.2191193518914336e-10; + q3 = -0.9999999897238481; + + //Tests for psSphereRot_CEOtoGCRS + //Return NULL for NULL earthPole input + { + out = psSphereRot_CEOtoGCRS(in); + ok( out == NULL, + "psSphereRot_CEOtoGCRS: return NULL for NULL earthPole input."); + in = psEarthPoleAlloc(); + } + in->x = 2.857175590089105e-4; + in->y = 2.3968739377734732e-5; + in->s = -1.3970066457904322e-8; + //Return correct psSphereRot for valid input + { + out = psSphereRot_CEOtoGCRS(in); + skip_start( out == NULL, 4, + "Skipping 4 tests because psSphereRot output is NULL!"); + ok_double_tol( out->q0, q0, 0.0001, + "psSphereRot_CEOtoGCRS: return correct psSphereRot for valid input (q0)."); + ok_double_tol( out->q1, q1, 0.0001, + "psSphereRot_CEOtoGCRS: return correct psSphereRot for valid input (q1)."); + ok_double_tol( out->q2, q2, 0.0001, + "psSphereRot_CEOtoGCRS: return correct psSphereRot for valid input (q2)."); + ok_double_tol( out->q3, -q3, 0.0001, + "psSphereRot_CEOtoGCRS: return correct psSphereRot for valid input (q3)."); + skip_end(); + } + + //Tests for psSphereRot_TEOtoCEO + psFree(out); + psTime *time = psTimeAlloc(PS_TIME_UT1); + psTime *time2 = psTimeAlloc(PS_TIME_UTC); + time->sec = 1049160600-1; + time->nsec = 657017200; + time->leapsecond = false; + psEarthPole *temp = psEarthPoleAlloc(); + temp->s = -2.0; + //Return NULL for NULL time input + { + out = psSphereRot_TEOtoCEO(NULL, NULL); + ok( out == NULL, + "psSphereRot_TEOtoCEO: return NULL for NULL time input."); + } + //Return NULL for invalid time input - large nsec, UTC time type + { + time2->nsec = 3e9; + out = psSphereRot_TEOtoCEO(time2, NULL); + ok( out == NULL, + "psSphereRot_TEOtoCEO: return NULL for invalid time input."); + } + //Return NULL for invalid time input - UT1 time type, large nsec + { + psFree(time2); + time2 = psTimeAlloc(PS_TIME_UT1); + time2->nsec = 1e9; + temp->s = 1.0; + out = psSphereRot_TEOtoCEO(time2, temp); + ok( out == NULL, + "psSphereRot_TEOtoCEO: return NULL for invalid time input."); + } + //Return NULL for invalid earthPole input + { + time2->nsec = 0; + psEarthPole *temp = psEarthPoleAlloc(); + temp->s = -2.0; + out = psSphereRot_TEOtoCEO(time2, temp); + ok( out == NULL, + "psSphereRot_TEOtoCEO: return NULL for invalid earthPole input."); + } + //Return correct psSphereRot for valid inputs + double x, y, z; + x = 0.01698625430807123; + y = -0.6616523084626379; + z = 0.7496183628158023; + { + psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time); + out = psSphereRot_TEOtoCEO(time, polarTideCorr); + skip_start( out == NULL, 3, + "Skipping 3 tests because psSphereRot output is NULL!"); + psSphere *obj = objSetup(); + psSphereRot *earthRot = psSphereRotConjugate(NULL, out); + psSphere *result = psSphereRotApply(NULL, earthRot, obj); + psCube *cube = psSphereToCube(result); + ok_double_tol( cube->x, x, 0.0001, + "psSphereRot_TEOtoCEO: return NULL for NULL time input. (x)"); + ok_double_tol( cube->y, y, 0.0001, + "psSphereRot_TEOtoCEO: return NULL for NULL time input. (y)"); + ok_double_tol( cube->z, z, 0.0001, + "psSphereRot_TEOtoCEO: return NULL for NULL time input. (z)"); + psFree(earthRot); + psFree(result); + psFree(cube); + psFree(obj); + skip_end(); + psFree(polarTideCorr); + } + + //Tests for psSphereRot_ITRStoTEO + //Return NULL for NULL earthPole input + { + psFree(out); + out = psSphereRot_ITRStoTEO(NULL); + ok( out == NULL, + "psSphereRot_ITRStoTEO: return NULL for NULL earthPole input."); + in = psEarthPoleAlloc(); + } + //Return correct psSphereRot for valid input + { + q0 = -1.0567571848664005e-6; + q1 = 3.218036562931509e-7; + q2 = -3.3580195807204483e-12; + q3 = -0.9999999999993899; + in->x = SEC_TO_RAD(-0.13275353774074533); + in->y = SEC_TO_RAD(0.4359436319739848); + in->s = SEC_TO_RAD(-4.2376965863576153e-10); + out = psSphereRot_ITRStoTEO(in); + skip_start( out == NULL, 4, + "Skipping 4 tests because psSphereRot output is NULL!"); + ok_double_tol( out->q0, q0, 0.0001, + "psSphereRot_ITRStoTEO: return correct psSphereRot for valid input (q0)."); + ok_double_tol( out->q1, q1, 0.0001, + "psSphereRot_ITRStoTEO: return correct psSphereRot for valid input (q1)."); + ok_double_tol( out->q2, q2, 0.0001, + "psSphereRot_ITRStoTEO: return correct psSphereRot for valid input (q2)."); + ok_double_tol( out->q3, q3, 0.0001, + "psSphereRot_ITRStoTEO: return correct psSphereRot for valid input (q3)."); + skip_end(); + } + + //Check for Memory leaks + { + psFree(temp); + psFree(time2); + psFree(time); + psFree(out); + psFree(in); + checkMem(); + } +} + +#define SPHERE_PRECESS_TP1_R 0.0 // 0.0 degrees +#define SPHERE_PRECESS_TP1_D 0.0 // 0.0 degrees +#define SPHERE_PRECESS_TP1_EXPECT_R 6.238453 // 357.437 degrees +#define SPHERE_PRECESS_TP1_EXPECT_D -0.019426 // -1.113 degrees +#define SPHERE_PRECESS_TP2_R 0.0 // 0.0 degrees +#define SPHERE_PRECESS_TP2_D 1.570796 // 90.0 degrees +#define SPHERE_PRECESS_TP2_EXPECT_R 6.260828 // 358.719 degrees +#define SPHERE_PRECESS_TP2_EXPECT_D 1.551353 // 88.886 degrees +#define SPHERE_PRECESS_TP3_R 3.141593 // 180.0 degrees +#define SPHERE_PRECESS_TP3_D 0.523599 // 30.0 degrees +#define SPHERE_PRECESS_TP3_EXPECT_R 3.096616 // 177.423 degrees +#define SPHERE_PRECESS_TP3_EXPECT_D 0.543024 // 31.113 degrees +#define ERROR_TOL 0.0001 + +void testSpherePrecess(void) +{ + diag(" >>>Test 4: psSpherePrecess"); + psSphereRot *rot = NULL; + psTime *fromTime = NULL; + psTime *toTime = NULL; + psSphere* inputCoord = psSphereAlloc(); + psSphere* outputCoord = NULL; + + //Return NULL for NULL time inputs + { + rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH); + ok( rot == NULL, + "psSpherePrecess: return NULL for NULL time inputs."); + } + //Return NULL for invalid mode input + { + fromTime = psTimeFromMJD(MJD_2100); + toTime = psTimeFromMJD(MJD_1900); + rot = psSpherePrecess(fromTime, toTime, -1); + ok( rot == NULL, + "psSpherePrecess: return NULL for invalid mode input."); + } + //Return correct psSphereRot for valid inputs, mode = ROUGH + { + inputCoord->r = SPHERE_PRECESS_TP1_R; + inputCoord->d = SPHERE_PRECESS_TP1_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + if (outputCoord->r < -0.0001) + outputCoord->r += 2.0 * M_PI; + skip_start( rot == NULL || outputCoord == NULL, 2, + "Skipping 2 tests because psSphereRot output is NULL!"); + ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP1_EXPECT_R, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (r)"); + ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP1_EXPECT_D, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (d)"); + skip_end(); + psFree(outputCoord); + psFree(rot); + } + //Return correct psSphereRot for valid inputs, mode = ROUGH + { + inputCoord->r = SPHERE_PRECESS_TP2_R; + inputCoord->d = SPHERE_PRECESS_TP2_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + if (outputCoord->r < -0.0001) + outputCoord->r += 2.0 * M_PI; + skip_start( rot == NULL || outputCoord == NULL, 2, + "Skipping 2 tests because psSphereRot output is NULL!"); + ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP2_EXPECT_R, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (r)"); + ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP2_EXPECT_D, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (d)"); + skip_end(); + psFree(outputCoord); + psFree(rot); + } + //Return correct psSphereRot for valid inputs, mode = ROUGH + { + inputCoord->r = SPHERE_PRECESS_TP3_R; + inputCoord->d = SPHERE_PRECESS_TP3_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + if (outputCoord->r < -0.0001) + outputCoord->r += 2.0 * M_PI; + skip_start( rot == NULL || outputCoord == NULL, 2, + "Skipping 2 tests because psSphereRot output is NULL!"); + ok_double_tol( outputCoord->r, SPHERE_PRECESS_TP3_EXPECT_R, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (r)"); + ok_double_tol( outputCoord->d, SPHERE_PRECESS_TP3_EXPECT_D, ERROR_TOL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_ROUGH mode. (d)"); + skip_end(); + psFree(outputCoord); + psFree(rot); + } + //Return correct psSphereRot for valid inputs, mode = COMPLETE_A + { + rot = psSpherePrecess(NULL, toTime, PS_PRECESS_COMPLETE_A); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + ok( outputCoord != NULL && rot != NULL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_COMPLETE_A mode."); + psFree(outputCoord); + psFree(rot); + } + //Return correct psSphereRot for valid inputs, mode = COMPLETE_B + { + rot = psSpherePrecess(NULL, toTime, PS_PRECESS_COMPLETE_B); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + ok( outputCoord != NULL && rot != NULL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_COMPLETE_B mode."); + psFree(outputCoord); + psFree(rot); + } + //Return correct psSphereRot for valid inputs, mode = IAU2000A + { + rot = psSpherePrecess(NULL, toTime, PS_PRECESS_IAU2000A); + outputCoord = psSphereRotApply(NULL, rot, inputCoord); + ok( outputCoord != NULL && rot != NULL, + "psSpherePrecess: return correct psSphereRot for valid" + " inputs and PS_PRECESS_IAU2000A mode."); + } + + //Check for Memory leaks + { + psFree(fromTime); + psFree(toTime); + psFree(outputCoord); + psFree(rot); + checkMem(); + } + +} + +