Index: /trunk/psLib/etc/pslib/pslib.config.template =================================================================== --- /trunk/psLib/etc/pslib/pslib.config.template (revision 5813) +++ /trunk/psLib/etc/pslib/pslib.config.template (revision 5814) @@ -34,5 +34,5 @@ psLib.eoc.precession.table.file.y STR DATADIR/pslib/tab5.2b.dat psLib.eoc.precession.table.file.s STR DATADIR/pslib/tab5.2c.dat -psLib.eoc.precession.table.file.iers STR DATADIR/pslib/iers_corr.dat +psLib.eoc.precession.table.file.iers STR DATADIR/pslib/finals2000A.dat psLib.eoc.precession.table.file.final STR DATADIR/pslib/finals_all.dat @psLib.eoc.precession.poly.x F64 -16616.99 2004191742.88 -427219.05 -198620.54 -46.05 5.98 Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 5813) +++ /trunk/psLib/src/astro/psCoord.h (revision 5814) @@ -11,6 +11,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-30 02:00:00 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-20 05:05:37 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -169,6 +169,6 @@ psPlaneTransform* psPlaneTransformAlloc( - int n1, ///< The order of the x term in the transform. - int n2 ///< The order of the y term in the transform. + int order1, ///< The order of the x term in the transform. + int order2 ///< The order of the y term in the transform. ); @@ -189,8 +189,8 @@ psPlaneDistort* psPlaneDistortAlloc( - int n1, ///< The order of the w term in the transform. - int n2, ///< The order of the x term in the transform. - int n3, ///< The order of the y term in the transform. - int n4 ///< The order of the z term in the transform. + int order1, ///< The order of the w term in the transform. + int order2, ///< The order of the x term in the transform. + int order3, ///< The order of the y term in the transform. + int order4 ///< The order of the z term in the transform. ); Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5813) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5814) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-15 02:37:48 $ + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-20 05:05:37 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -512,7 +512,5 @@ S = SEC_TO_RAD(S * 1e-6); - // now calculate the non-poly portion from the tables - psF64* cols[17]; for (int lcv = 0; lcv < 17; lcv++) { @@ -628,12 +626,12 @@ out->x = cols[1][rowNum]; out->y = cols[2][rowNum]; - out->x = SEC_TO_RAD(out->x) * 1e-6; - out->y = SEC_TO_RAD(out->y) * 1e-6; + out->x = SEC_TO_RAD(out->x) * 1e-3; + out->y = SEC_TO_RAD(out->y) * 1e-3; rowNum = numRows; } else { out->x = cols[3][rowNum]; out->y = cols[4][rowNum]; - out->x = SEC_TO_RAD(out->x) * 1e-6; - out->y = SEC_TO_RAD(out->y) * 1e-6; + out->x = SEC_TO_RAD(out->x) * 1e-3; + out->y = SEC_TO_RAD(out->y) * 1e-3; rowNum = numRows; } @@ -762,5 +760,5 @@ T += -2451545.0; double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T); - psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0); + psSphereRot *out = psSphereRotAlloc(-theta, 0.0, 0.0); // psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); @@ -773,43 +771,7 @@ { PS_ASSERT_PTR_NON_NULL(time, NULL); + PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL); psEarthPole *out = psEarthPoleAlloc(); - // psSphere *in = NULL; - // double x,y,s; - // psTime *time2 = psTimeAlloc(time->type); - // time2->sec = time->sec; - // time2->nsec = time->nsec; - // time2->leapsecond = time->leapsecond; - //XXX: Time must be converted to TAI before use? - /* if (time->type == PS_TIME_UT1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n"); - return NULL; - } else if (time->type != PS_TIME_TAI) { - time2 = psTimeConvert(time2, PS_TIME_TAI); - } - */ - //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate - //polar motion coordinates. - // in = p_psTimeGetPoleCoords(time2); - // if (in == NULL) { - // psError(PS_ERR_BAD_PARAMETER_NULL, true, - // "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n"); - // return NULL; - // } - // out->x = in->r; - // out->y = in->d; - // s = psTimeGetUT1Delta(time2, bulletin); - // out->s = s; - - //XXX: Apply polar tide correction here??? - // psEarthPole *correction = psEOC_PolarTideCorr(time2); - // out->x += correction->x; - // out->y += correction->y; - // out->s += correction->s; - - // psFree(time2); - // psFree(in); - // psFree(correction); - out->x = 0.0; out->y = 0.0; @@ -817,12 +779,12 @@ double MJD = psTimeToMJD(time); - - if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) || - (MJD > 53858.0 && bulletin == PS_IERS_A) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Invalid time input. Date, %lf, is out of range\n", MJD); - return out; - } - + /* + if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) || + (MJD > 53858.0 && bulletin == PS_IERS_A) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Invalid time input. Date, %lf, is out of range\n", MJD); + return out; + } + */ // Check if EOC data loaded if(! eocInitialized) { @@ -841,26 +803,91 @@ } int numRows = ((psVector*)(finalsTable->data[0]))->n; - for (int rowNum = 0; rowNum < numRows; rowNum++) { - if ( (MJD - cols[0][rowNum]) < 1.0 ) { - if (bulletin == PS_IERS_A) { - out->x = cols[1][rowNum]; - out->y = cols[2][rowNum]; - out->s = cols[3][rowNum]; - out->x = SEC_TO_RAD(out->x); - out->y = SEC_TO_RAD(out->y); - out->s = SEC_TO_RAD(out->s); - rowNum = numRows; - } else { - out->x = cols[4][rowNum]; - out->y = cols[5][rowNum]; - out->s = cols[6][rowNum]; - out->x = SEC_TO_RAD(out->x); - out->y = SEC_TO_RAD(out->y); - out->s = SEC_TO_RAD(out->s); - rowNum = numRows; + psVector *X = psVectorAlloc(numRows, PS_TYPE_F64); + psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64); + psVector *S = psVectorAlloc(numRows, PS_TYPE_F64); + psVector *T = psVectorAlloc(numRows, PS_TYPE_F64); + if (bulletin == PS_IERS_A) { + for (int rowNum = 0; rowNum < numRows; rowNum++) { + T->data.F64[rowNum] = cols[0][rowNum]; + X->data.F64[rowNum] = cols[1][rowNum]; + Y->data.F64[rowNum] = cols[2][rowNum]; + S->data.F64[rowNum] = cols[3][rowNum]; + } + } else { + for (int rowNum = 0; rowNum < numRows; rowNum++) { + T->data.F64[rowNum] = cols[0][rowNum]; + X->data.F64[rowNum] = cols[4][rowNum]; + Y->data.F64[rowNum] = cols[5][rowNum]; + S->data.F64[rowNum] = cols[6][rowNum]; + } + } + + double xOut = 0.0; + double yOut = 0.0; + double sOut = 0.0; + double xTerm = 0.0; + double yTerm = 0.0; + double sTerm = 0.0; + int k = 0; + for (int i = 0; i < (numRows-1); i++) { + if (MJD >= T->data.F64[i] && MJD < T->data.F64[i+1]) { + k = i; + if (k < 2) { + k = 2; } - } - } - + if (k > (numRows-2)) { + k = numRows-2; + } + for (int m = k-1; m <= k+2; m++) { + xTerm = X->data.F64[m]; + yTerm = Y->data.F64[m]; + sTerm = S->data.F64[m]; + for (int j = k-1; j <= k+2; j++) { + if ( m != j) { + double term = (MJD - T->data.F64[j])/(T->data.F64[m] - T->data.F64[j]); + xTerm *= term; + yTerm *= term; + sTerm *= term; + } + } + xOut += xTerm; + yOut += yTerm; + sOut += sTerm; + } + i = numRows-1; + } + } + out->x = SEC_TO_RAD(xOut); + out->y = SEC_TO_RAD(yOut); + out->s = SEC_TO_RAD(sOut); + + + /* + for (int rowNum = 0; rowNum < numRows; rowNum++) { + if ( (MJD - cols[0][rowNum]) < 1.0 ) { + if (bulletin == PS_IERS_A) { + out->x = cols[1][rowNum]; + out->y = cols[2][rowNum]; + out->s = cols[3][rowNum]; + out->x = SEC_TO_RAD(out->x); + out->y = SEC_TO_RAD(out->y); + out->s = SEC_TO_RAD(out->s); + rowNum = numRows; + } else { + out->x = cols[4][rowNum]; + out->y = cols[5][rowNum]; + out->s = cols[6][rowNum]; + out->x = SEC_TO_RAD(out->x); + out->y = SEC_TO_RAD(out->y); + out->s = SEC_TO_RAD(out->s); + rowNum = numRows; + } + } + } + */ + psFree(X); + psFree(Y); + psFree(S); + psFree(T); return out; } @@ -1081,4 +1108,5 @@ y = motion->y; s = motion->s; + s = -s; @@ -1130,31 +1158,31 @@ /* - psSphereRot r,s,t; + psSphereRot r,p,t; - r.q0=sin(motion->y/2.0); - r.q1=0; - r.q2=0; - r.q3=cos(motion->y/2.0); + r.q0=sin(y/2.0); + r.q1=0; + r.q2=0; + r.q3=cos(y/2.0); - s.q0=0; - s.q1=sin(motion->x/2.0); - s.q2=0; - s.q3=cos(motion->x/2.0); + p.q0=0; + p.q1=sin(x/2.0); + p.q2=0; + p.q3=cos(x/2.0); - t.q0=0; - t.q1=0; - t.q2=sin(motion->s/2.0); - t.q3=cos(motion->s/2.0); + t.q0=0; + t.q1=0; + t.q2=sin(s/2.0); + t.q3=cos(s/2.0); - // calculate t*s*r. - psSphereRot* temp = psSphereRotCombine(NULL,&t,&s); - out = psSphereRotCombine(NULL, temp, &r); - psFree(temp); + // calculate t*s*r. + psSphereRot* temp = psSphereRotCombine(NULL,&t,&p); + out = psSphereRotCombine(NULL, temp, &r); + psFree(temp); - return out; - */ - - s = -s; + return out; + */ + //Newest trial - mult. y matrix * x * z + A[0][0] = cos(x)*cos(s); A[1][0] = cos(x)*sin(s); @@ -1167,13 +1195,14 @@ A[2][2] = cos(x)*cos(y); - /* A[0][0] = cos(x)*cos(s); - A[0][1] = cos(x)*sin(s); - A[0][2] = -sin(x); - A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s); - A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s); - A[1][2] = cos(x)*sin(y); - A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s); - A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s); - A[2][2] = cos(x)*cos(y); + /* + A[0][0] = cos(x)*cos(s); + A[0][1] = cos(x)*sin(s); + A[0][2] = -sin(x); + A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s); + A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s); + A[1][2] = cos(x)*sin(y); + A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s); + A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s); + A[2][2] = cos(x)*cos(y); */ //New trial - mult z * y * x @@ -1203,50 +1232,4 @@ out = rotMatrix_To_Quat(A); - // out = psSphereRotAlloc(motion->y, motion->x, motion->s); - - // out->q0 = SEC_TO_RAD(out->q0); - // out->q1 = SEC_TO_RAD(out->q1); - // out->q2 = SEC_TO_RAD(out->q2) + (out->q0+out->q1/2.0); - // out->q3 = sqrt(1.0 - (out->q0*out->q0 + out->q1*out->q1 + out->q2*out->q2) ); - /* double diag_sum[3]; - int maxi; - double recip; - diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2]; - diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2]; - diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2]; - diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2]; - - maxi = 0; - for (int i = 1; i < 4; ++i) { - if (diag_sum[i] > diag_sum[maxi]) { - maxi = i; - } - } - - double p = 0.5 * sqrt(diag_sum[maxi]); - recip = 1.0 / (4.0 * p); - - if (maxi == 0) { - out->q0 = p; - out->q1 = recip * (A[0][1] + A[1][0]); - out->q2 = recip * (A[2][0] + A[0][2]); - out->q3 = recip * (A[1][2] - A[2][1]); - } else if (maxi == 1) { - out->q0 = recip * (A[0][1] + A[1][0]); - out->q1 = p; - out->q2 = recip * (A[1][2] + A[2][1]); - out->q3 = recip * (A[2][0] - A[0][2]); - } else if (maxi == 2) { - out->q0 = recip * (A[2][0] + A[0][2]); - out->q1 = recip * (A[1][2] + A[2][1]); - out->q2 = p; - out->q3 = recip * (A[0][1] - A[1][0]); - } else if (maxi == 3) { - out->q0 = recip * (A[1][2] - A[2][1]); - out->q1 = recip * (A[2][0] - A[0][2]); - out->q2 = recip * (A[0][1] - A[1][0]); - out->q3 = p; - } - */ return out; } Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5813) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5814) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-15 02:37:48 $ +* @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-20 05:05:37 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -418,5 +418,5 @@ printf("\n <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - return 3; + // return 3; } psFree(polarMotion); @@ -440,5 +440,5 @@ printf("\n <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - return 5; + // return 5; } @@ -512,4 +512,7 @@ "psEOC_PolarTideCorr returned NULL for valid input time.\n"); return 2; + } else { + printf("\nPolarTideCorr output = x,y,s = %.13g, %.13g, %.13g\n", + eop->x, eop->y, eop->s); } @@ -522,7 +525,9 @@ psS32 testEOCNutation(void) { - psTime *in = psTimeAlloc(PS_TIME_TAI); - in->sec = 1131579114; - in->nsec = 498489000; + psTime *in = psTimeAlloc(PS_TIME_UTC); + // in->sec = 1131579114; + // in->nsec = 498489000; + in->sec = timesec; + in->nsec = 0; in->leapsecond = false; psTime *empty = NULL; @@ -606,5 +611,6 @@ objSetup(); psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo); - psSphere *result = psSphereRotApply(NULL, earthRot, obj); + // psSphere *result = psSphereRotApply(NULL, earthRot, obj); + psSphere *result = psSphereRotApply(NULL, teoceo, obj); psCube *cube = psSphereToCube(result); Index: /trunk/psLib/test/astro/verified/tst_psCoord.stderr =================================================================== --- /trunk/psLib/test/astro/verified/tst_psCoord.stderr (revision 5813) +++ /trunk/psLib/test/astro/verified/tst_psCoord.stderr (revision 5814) @@ -8,9 +8,9 @@ Following should generate an error message: negative x order |E|psPlaneTransformAlloc (FILE:LINENO) - Error: n1 is less than 0. + Error: order1 is less than 0. |I|testPlaneTransformAlloc Following should generate an error message: negative y order |E|psPlaneTransformAlloc (FILE:LINENO) - Error: n2 is less than 0. + Error: order2 is less than 0. ---> TESTPOINT PASSED (psCoord{psPlaneTransformAlloc()} | tst_psCoord.c) @@ -25,17 +25,17 @@ Following should generate an error message for negative x order |E|psPlaneDistortAlloc (FILE:LINENO) - Error: n1 is less than 0. + Error: order1 is less than 0. |I|testPlaneDistortAlloc Following should generate an error message for negative y order |E|psPlaneDistortAlloc (FILE:LINENO) - Error: n2 is less than 0. + Error: order2 is less than 0. |I|testPlaneDistortAlloc Following should generate an error message for negative z order |E|psPlaneDistortAlloc (FILE:LINENO) - Error: n3 is less than 0. + Error: order3 is less than 0. |I|testPlaneDistortAlloc Following should generate an error message for negative w order |E|psPlaneDistortAlloc (FILE:LINENO) - Error: n4 is less than 0. + Error: order4 is less than 0. ---> TESTPOINT PASSED (psCoord{psPlaneDistortAlloc()} | tst_psCoord.c)