Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5532) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5533) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-17 00:07:25 $ + * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-17 03:59:05 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -493,5 +493,6 @@ S -= X*Y/2.0; - psEarthPole* pole = psAlloc(sizeof(psEarthPole)); + // psEarthPole* pole = psAlloc(sizeof(psEarthPole)); + psEarthPole* pole = psEarthPoleAlloc(); pole->x = X; pole->y = Y; @@ -850,5 +851,116 @@ psEarthPole* psEOC_NutationCorr(psTime *time) { - return NULL; + // Check for null parameter + PS_ASSERT_PTR_NON_NULL(time, NULL); + if (time->type == PS_TIME_UT1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Invalid time input. Time cannot be of type UT1.\n"); + return NULL; + } + + // Convert psTime to MJD + double MJD = psTimeToMJD(time); + + // Calculate number of Julian centuries since 2000 + double t = ( MJD - MJD_2000 ) / JULIAN_CENTURY; + double t2 = t*t; + double t3 = t*t*t; + double t4 = t*t*t*t; + + //XXX: I think the t's should be inside of the SEC_TO_RAD conversion. + //Check this and for Precession Model as well! + double F[5]; + F[0] = DEG_TO_RAD(134.96340251) + + SEC_TO_RAD(1717915923.2178)*t + + SEC_TO_RAD(31.8792)*t2 + + SEC_TO_RAD(0.051635)*t3 - + SEC_TO_RAD(0.00024470)*t4; + + // Mean Anomaly of the Sun + F[1] = DEG_TO_RAD(357.52910918) + + SEC_TO_RAD(129596581.0481)*t - + SEC_TO_RAD(0.5532)*t2 + + SEC_TO_RAD(0.000136)*t3 - + SEC_TO_RAD(0.00001149)*t4; + + // L − Omega + F[2] = DEG_TO_RAD(93.27209062) + + SEC_TO_RAD(1739527262.8478)*t - + SEC_TO_RAD(12.7512)*t2 - + SEC_TO_RAD(0.001037)*t3 + + SEC_TO_RAD(0.00000417)*t4; + + // Mean Elongation of the Moon from the Sun + F[3] = DEG_TO_RAD(297.85019547) + + SEC_TO_RAD(1602961601.2090)*t - + SEC_TO_RAD(6.3706)*t2 + + SEC_TO_RAD(0.006593)*t3 - + SEC_TO_RAD(0.00003169)*t4; + + // Mean Longitude of the Ascending Node of the Moon + F[4] = DEG_TO_RAD(125.04455501) - + SEC_TO_RAD(6962890.5431)*t + + SEC_TO_RAD(7.4722)*t2 + + SEC_TO_RAD(0.007702)*t3 - + SEC_TO_RAD(0.0000593)*t4; + + //argument values taken from table 5.1 in IERS techical note No.32 + //http://maia.usno.navy.mil/conv2000/chapter5/tn32_c5.pdf, p38 + //Units are in micro-arcseconds here and must be converted to radians before using + double w_l[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-1.0), SEC_TO_RAD(1.0), 0.0, 0.0, + SEC_TO_RAD(-1.0), 0.0, 0.0, 0.0, SEC_TO_RAD(1.0)}; + double w_l_p[10] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}; + double w_F[10] = {SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), + SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0}; + double w_D[10] = {0.0,0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0, SEC_TO_RAD(2.0), 0.0,0.0,0.0}; + double w_Omega[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), + SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-1.0), 0.0}; + double xp_sin[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44), + SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76), + SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)}; + double xp_cos[10] = {SEC_TO_RAD(0.25), SEC_TO_RAD(1.32), SEC_TO_RAD(0.25), SEC_TO_RAD(1.23), + SEC_TO_RAD(6.52), SEC_TO_RAD(-0.48), SEC_TO_RAD(2.73), SEC_TO_RAD(-8.19), + SEC_TO_RAD(-1.11), SEC_TO_RAD(-0.43)}; + double yp_sin[10] = {SEC_TO_RAD(-0.25), SEC_TO_RAD(-1.32), SEC_TO_RAD(-0.25), + SEC_TO_RAD(-1.23), SEC_TO_RAD(-6.52), SEC_TO_RAD(0.48), SEC_TO_RAD(-2.73), + SEC_TO_RAD(8.19), SEC_TO_RAD(1.11), SEC_TO_RAD(0.43)}; + double yp_cos[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44), + SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76), + SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)}; + + double X = 0.0; + double Y = 0.0; + double arg = 0.0; + //This is from eqn 131 in the ADD - Note: pj_tj isn't included the first time. + //XXX: The xp_sin, yp_cos, etc. may need to be multiplied by pow(t,i) here? adding now... + double tj = 0.0; + + // calculate the polynomial portion first - the pj * t^j (poly coeff's) + // Check if EOC data loaded + if(! eocInitialized) { + eocInitialized = p_psEOCInit(); + if(!eocInitialized) { + // XXX: Move error message. + psError(PS_ERR_UNKNOWN, false, + "Could not initialize EOC tables -- check data files."); + return NULL; + } + } + X = psPolynomial1DEval(xPoly,t); + Y = psPolynomial1DEval(yPoly,t); + for (int i = 0; i < 10; i++) { + // tj = SEC_TO_RAD(pow(t, i)); + tj = pow(t, i); + arg = w_l[i]*F[0] + w_l_p[i]*F[1] + w_F[i]*F[2] + w_D[i]*F[3] + w_Omega[i]*F[4]; + X += xp_sin[i] * tj * sin(arg) + xp_cos[i] * tj * cos(arg); + Y += yp_sin[i] * tj * sin(arg) + yp_cos[i] * tj * cos(arg); + } + + psEarthPole *pole = psEarthPoleAlloc(); + pole->x = X; + pole->y = Y; + pole->s = SEC_TO_RAD(4.7e-5) * t; //XXX: This conv. should include t? + + return pole; } Index: /trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.h (revision 5532) +++ /trunk/psLib/src/astro/psEarthOrientation.h (revision 5533) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-16 02:40:13 $ +* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-17 03:59:05 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -45,4 +45,10 @@ psPrecessMethod; +/** Initialize EOC. + * + * Parses IERS table data and creates the polynomials used by certain EOC functions. + * + * @return bool: True if successful, otherwise false. + */ bool p_psEOCInit(void); Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5532) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5533) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-16 20:52:23 $ +* @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-17 03:59:05 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -14,6 +14,4 @@ #include "pslib_strict.h" -//static psS32 testEOCParallax(void); -static psS32 testEOCInit_Final(void); static psS32 testAberration(void); static psS32 testGravityDeflect(void); @@ -26,8 +24,6 @@ testDescription tests[] = { - {testEOCInit_Final, 665, "psEOCInit_Final", 0, false}, {testAberration, 666, "psAberration()", 0, false}, {testGravityDeflect, 667, "psGravityDeflect()", 0, false}, - // {testEOCParallax, 668, "psEOCParallax()", 0, false}, {testEOCPrecession, 669, "psEOCPrecession()", 0, false}, {testEOCPolar, 670, "psEOCPolar()", 0, false}, @@ -56,9 +52,4 @@ // return ( ! runTestSuite( stderr, "psEarthOrientation", tests, argc, argv ) ); -} - -psS32 testEOCInit_Final(void) -{ - return 0; } @@ -171,5 +162,4 @@ psFree(UT1time); //Check return values from valid precession input - pmodel = psEOC_PrecessionModel(time); if ( pmodel == NULL ) { @@ -288,4 +278,44 @@ psS32 testEOCNutation(void) { + psTime *in = psTimeAlloc(PS_TIME_TAI); + in->sec = 1131579114; + in->nsec = 498489000; + in->leapsecond = false; + psTime *empty = NULL; + psEarthPole *nute = NULL; + + //Return NULL for NULL input time. + 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_PrecessionModel(in); + if ( nute == NULL ) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_NutationCorr returned NULL for valid input.\n"); + return 3; + } else { + printf("Nutation Correction output = x,y,s = %.8g, %.8g, %.8g\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 0; Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5532) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5533) @@ -1,11 +1,2 @@ -/***************************** TESTPOINT ******************************************\ -* TestFile: tst_psEarthOrientation.c * -* TestPoint: psEarthOrientation{psEOCInit_Final} * -* TestType: Positive * -\**********************************************************************************/ - - ----> TESTPOINT PASSED (psEarthOrientation{psEOCInit_Final} | tst_psEarthOrientation.c) - /***************************** TESTPOINT ******************************************\ * TestFile: tst_psEarthOrientation.c * @@ -70,4 +61,8 @@ \**********************************************************************************/ +