Index: /trunk/psLib/psTableParse.c =================================================================== --- /trunk/psLib/psTableParse.c (revision 6029) +++ /trunk/psLib/psTableParse.c (revision 6030) @@ -7,6 +7,6 @@ * @author Dave Robbins, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-12 21:18:48 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-18 00:41:29 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -44,4 +44,12 @@ input = fopen(argv[1], "r"); output = fopen(argv[2], "w"); + if (input == NULL) { + printf("\nERROR. Could not open input file.\n"); + return 0; + } + if (output == NULL) { + printf("\nERROR. Could not open or create output file.\n"); + return 0; + } char data2[200]; printHeader(argv[1]); Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 6029) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 6030) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-12 20:40:12 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-18 00:41:29 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -39,5 +39,5 @@ #define PS_AU 149597890000.0 /* meters */ // Modified Julian Day 1/1/2000 00:00:00 -#define MJD_2000 51544.0 +#define MJD_2000 51544.5 // Days in Julian century #define JULIAN_CENTURY 36525.0 @@ -299,11 +299,9 @@ double a = 0.0; - //mu = apparent * direction; I believe this should be: - // mu = apparent->r * direction->r + apparent->d * direction->d; psCube* directionVector = psSphereToCube(direction); if (directionVector == NULL) printf("directionVector is null\n"); psCube* actualVector = psSphereToCube(actual); - if (directionVector == NULL) + if (actualVector == NULL) printf("actualVector is null\n"); // mu = acos(directionVector->x*actualVector->x + @@ -311,13 +309,9 @@ directionVector->y*actualVector->y + directionVector->z*actualVector->z); - - // rp = apparent - mu * direction; - // rp->r = actual->r - mu * direction->r; - // rp->d = actual->d - mu * direction->d; rp->x = actualVector->x - mu * directionVector->x; rp->y = actualVector->y - mu * directionVector->y; rp->z = actualVector->z - mu * directionVector->z; - mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu)); + mu_p = mu - speed * ((mu * mu - 1.0) / (1.0 - speed * mu)); //Not sure if this is right. ADD gets a scalar from division of rp (a vector) @@ -345,7 +339,7 @@ //XXX: Must be a sign error somewhere above? Magnitude of change is correct but wrong way... //This will fix the problem but is somewhat of a hack. - r_p->x = actualVector->x - (r_p->x - actualVector->x); - r_p->y = actualVector->y - (r_p->y - actualVector->y); - r_p->z = actualVector->z - (r_p->z - actualVector->z); + // r_p->x = actualVector->x - (r_p->x - actualVector->x); + // r_p->y = actualVector->y - (r_p->y - actualVector->y); + // r_p->z = actualVector->z - (r_p->z - actualVector->z); psSphere *r_pSphere = psCubeToSphere(r_p); @@ -370,4 +364,5 @@ PS_ASSERT_PTR_NON_NULL(sun, NULL); + psSphere *temp = psSphereAlloc(); // calculating the apparent angle from the actual angle and the sun position @@ -378,11 +373,21 @@ // N.B., assuming the psSphereToCube function returns a unit vector. // double theta = acos(sunVector->x*actualVector->x + - double theta = (sunVector->x*actualVector->x + - sunVector->y*actualVector->y + - sunVector->z*actualVector->z); - printf(" Theta = %lf\n", theta); + double dotProd = (sunVector->x*actualVector->x + + sunVector->y*actualVector->y + sunVector->z*actualVector->z); + double theta, sunMag, actMag; + sunMag = sqrt(sunVector->x*sunVector->x + sunVector->y*sunVector->y + + sunVector->z*sunVector->z); + actMag = sqrt(actualVector->x*actualVector->x + actualVector->y*actualVector->y + + actualVector->z*actualVector->z); + dotProd = dotProd / (sunMag * actMag); + theta = acos(dotProd); + + // theta = acos( cos(sun->d)*cos(actual->d)*cos(sun->r-actual->r) + sin(sun->r)*sin(actual->r) ); + + printf(" Theta = %.13g\n", theta); // theta = acos(-theta); // printf("\n Theta = %lf\n", theta); double r0 = PS_AU * tan(theta); + printf(" r0 = %.19e\n", r0); double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0); @@ -390,5 +395,5 @@ double limit = SEC_TO_RAD(1.75); printf(" deflection = %.13g\n", deflection); - printf(" limit = %lf\n", limit); + //printf(" limit = %lf\n", limit); if (deflection > limit) { // deflection = limit; @@ -400,13 +405,16 @@ } - if (apparent == NULL) { - apparent = psSphereAlloc(); - } else { - apparent->r = 0.0; - apparent->d = 0.0; - apparent->rErr = 0.0; - apparent->dErr = 0.0; - } - + /* if (apparent == NULL) { + apparent = psSphereAlloc(); + } else { + apparent->r = 0.0; + apparent->d = 0.0; + apparent->rErr = 0.0; + apparent->dErr = 0.0; + } + */ + if (apparent != NULL) { + psFree(apparent); + } // bend the actual vector away from the sun vector by deflection angle. // XXX: Not sure how to do this. Dave thinks the formula should be: @@ -414,17 +422,37 @@ theta = 0.0; double phi = 0.0; - // deflection = SEC_TO_RAD(deflection); + // deflection = SEC_TO_RAD(deflection) * 1e6; + // deflection *= M_PI * 1e-2; theta = atan(r0/PS_AU) * tan(deflection); - printf(" Theta = %.13g\n", theta); - printf(" deflection = %.13g\n", deflection); // phi = sqrt( deflection*deflection - theta*theta ); - phi = deflection * cos(asin(theta/deflection)) * 3e-2; + // phi = deflection * cos(asin(theta/deflection)); + + // phi = sqrt(theta*theta - deflection*deflection); + // phi = deflection * cos(asin(theta/deflection)) * 3e-2; // phi = cos(asin(theta/deflection)); // phi = asin(theta/deflection); - apparent->r = theta; - apparent->d = phi; + + // apparent->r = theta; + // apparent->d = phi; + /* + actualVector->x += actualVector->x*deflection; + actualVector->y += actualVector->y*deflection; + actualVector->z += actualVector->z*deflection; + apparent = psCubeToSphere(actualVector); + */ + theta = tan(sun->r - actual->r) * deflection; + phi = tan(sun->d - actual->d) * deflection; + + printf(" Theta = %.13g\n", theta); + printf(" phi = %.13g\n", phi); + + temp->r = theta; + temp->d = phi; + apparent = psSphereSetOffset(actual, temp, PS_SPHERICAL, PS_RADIAN); psFree(actualVector); psFree(sunVector); + psFree(temp); + return apparent; } @@ -807,5 +835,6 @@ } -psSphereRot* psSphereRot_TEOtoCEO(const psTime *time) +psSphereRot* psSphereRot_TEOtoCEO(const psTime *time, + psEarthPole *tidalCorr) { PS_ASSERT_PTR_NON_NULL(time,NULL); @@ -815,10 +844,17 @@ in = psTimeConvert(in, PS_TIME_UT1); } - // double T = (double)(in->sec) + (double)(in->nsec / 1e9); + if (tidalCorr != NULL && tidalCorr->s != 0.0) { + int nsec = in->nsec + (int)(tidalCorr->s * 1e9); + if (nsec < 0.0) { + in->sec += -1; + in->nsec = (int)(1e9) - nsec; + } else { + in->nsec = nsec; + } + } double T = psTimeToJD(in); T += -2451545.0; double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T); psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0); - // psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); psFree(in); @@ -919,4 +955,9 @@ out->x = SEC_TO_RAD(xOut); out->y = SEC_TO_RAD(yOut); + // psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time); + // out->x += polarTideCorr->x; + // out->y += polarTideCorr->y; + // psFree(polarTideCorr); + // out->s = SEC_TO_RAD(sOut); @@ -1034,5 +1075,5 @@ CORX = SEC_TO_RAD(CORX); CORY = SEC_TO_RAD(CORY); - CORZ = SEC_TO_RAD(CORZ); + // CORZ = SEC_TO_RAD(CORZ); out->x = CORX; Index: /trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.h (revision 6029) +++ /trunk/psLib/src/astro/psEarthOrientation.h (revision 6030) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-17 03:59:05 $ +* @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-18 00:41:29 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -62,11 +62,11 @@ psEarthPole *psEarthPoleAlloc(void); -/** Calculates the actual position of a star, given its apparent position and the +/** Calculates the apparent position of a star, given its actual position and the * velocity vector of the observer. * * The actual and apparent positions are represented as psSphere entries, as is the * direction of motion. The speed in that direction is given in units of the speed - * of light. If the value of actual is NULL, a new psSphere is allocated, otherwise - * the point to actual is used for the result. + * of light. If the value of apparent is NULL, a new psSphere is allocated, otherwise + * the point to apparent is used for the result. * * @return psSphere*: the actual position of a star. @@ -79,10 +79,10 @@ ); -/** Calculates the actual position of a star, given its apparent position and the +/** Calculates the apparent position of a star, given its actual position and the * position of the sun. * * The actual and apparent positions are represented as psSphere entries, as is - * position of the sun. If the value of actual is NULL, a new psSphere is allocated, - * otherwise the point to actual is used for the result. + * position of the sun. If the value of apparent is NULL, a new psSphere is allocated, + * otherwise the point to apparent is used for the result. * * @return psSphere*: the apparent position of a star. @@ -133,5 +133,6 @@ */ psSphereRot *psSphereRot_TEOtoCEO( - const psTime *time ///< specified time + const psTime *time, ///< specified time + psEarthPole *tidalCorr ///< UT1 polar tide correction or NULL ); Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 6029) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 6030) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-01-12 20:40:13 $ +* @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-18 00:41:29 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -85,8 +85,8 @@ { double r1, r2, d1, d2; - r1 = in1->r; - r2 = in2->r; - d1 = in1->d; - d2 = in2->d; + d1 = in1->r; + d2 = in2->r; + r1 = in1->d; + r2 = in2->d; double out = 0.0; double c1, c2, cd, s1, s2, sum, ac; @@ -101,5 +101,6 @@ ac = acos(sum); - printf("\n c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac); + out = ac; + // printf("\n c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac); // out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2)); // if (out != 0.0) printf("\n in greatCircle, out = %lf\n", out); @@ -109,17 +110,7 @@ psS32 testAberration(void) { - psSphere *apparent = NULL; - // psSphere *actual = psSphereAlloc(); - // psSphere *direction = psSphereAlloc(); psSphere *empty = NULL; - psCube *actualCube = psCubeAlloc(); - // actual->r = DEG_TO_RAD(122.9153182445501); - // actual->d = DEG_TO_RAD(48.562968978679194); - // actualCube->x = -0.3596195125758298; - // actualCube->y = 0.5555613903455866; - // actualCube->z = 0.7496834983724809; - //values from After gravity deflection// actualCube->x = -0.35961949760293604; @@ -127,5 +118,4 @@ actualCube->z = 0.7496835020836093; psSphere *actual = psCubeToSphere(actualCube); - psCube *cubeDir = psCubeAlloc(); cubeDir->x = 5148.713262821658; @@ -135,15 +125,8 @@ cubeDir->y += 248.46429758174693; cubeDir->z += 0.09694774143797581; - double length = sqrt(cubeDir->x*cubeDir->x+cubeDir->y*cubeDir->y+cubeDir->z*cubeDir->z); - cubeDir->x = cubeDir->x/length; - cubeDir->y = cubeDir->y/length; - cubeDir->z = cubeDir->z/length; psSphere *direction = psCubeToSphere(cubeDir); - double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z); - // Speed of light in vacuum (src:NIST) - double c = 299792458.0; /* m/s */ - // speed = speed / c; - speed = length / c; + double c = 299792458.0; // Speed of light in vacuum (src:NIST) /* m/s */ + speed /= c; empty = psAberration(empty, apparent, direction, speed); @@ -161,9 +144,6 @@ apparent = psAberration(apparent, actual, direction, speed); - - // printf("\nSphere Difference = r,d = %.13g, %.13g\n", - // (actual->r - apparent->r), (actual->d - apparent->d)); psCube *outCube = psSphereToCube(apparent); - printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", + printf("\n -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube->x, outCube->y, outCube->z); //expected cube values @@ -173,9 +153,8 @@ z = 0.7497078321908413; - printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g -- \n", x, y, z); + printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g -- \n", x, y, z); printf("Cube Difference = x,y,z = %.13g, %.13g, %.13g \n\n", (x - outCube->x), (y - outCube->y), (z - outCube->z) ); psFree(actual); - // psFree(actualCube); actualCube->x = x; actualCube->y = y; @@ -190,11 +169,10 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psAberration returned incorrect values.\n"); - printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n", - (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) ); - printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n", - (actualCube->x - outCube->x), (actualCube->y - outCube->y), - (actualCube->z - outCube->z) ); + // printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n", + // (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) ); + // printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n", + // (actualCube->x - outCube->x), (actualCube->y - outCube->y), + // (actualCube->z - outCube->z) ); psFree(actual); - // psFree(actualCube); actualCube->x = x; actualCube->y = y; @@ -209,5 +187,4 @@ psFree(outCube); psFree(actualCube); - psFree(cubeDir); psFree(apparent); @@ -220,10 +197,6 @@ psS32 testGravityDeflect(void) { - - // psSphere *actual = psSphereAlloc(); psSphere *apparent = NULL; - // psSphere *sun = psSphereAlloc(); psSphere *empty = NULL; - psCube *actualCube = psCubeAlloc(); actualCube->x = -0.3596195125758298; @@ -235,5 +208,12 @@ 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; + printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n", sunCube->x, sunCube->y, sunCube->z); psSphere *sun = psCubeToSphere(sunCube); + printf("sunSphere = r, d = %.13g, %.13g\n", sun->r, sun->d); + psCube *outCube = psCubeAlloc(); empty = psGravityDeflection(apparent, empty, sun); @@ -251,15 +231,16 @@ apparent = psGravityDeflection(NULL, actual, sun); - psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); - apparent->r *= -1.0; - apparent->d *= -1.0; - psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + // 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); 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(result2); + // psCube *outCube = psSphereToCube(result); + // printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", + // outCube->x, outCube->y, outCube->z); + psCube *outCube2 = psSphereToCube(apparent); printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g -- \n", outCube2->x, outCube2->y, outCube2->z); @@ -275,20 +256,23 @@ // psError(PS_ERR_BAD_PARAMETER_VALUE, false, // "psGravityDeflection returned incorrect values.\n"); - printf("expect-actual= x,y,z = %.11g, %.11g, %.11g \n", + printf("expect-actual= x,y,z = %.13g, %.13g, %.13g \n", (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) ); - printf("result-actual= x,y,z = %.11g, %.11g, %.11g \n", - (outCube->x - actualCube->x), (outCube->y - actualCube->y), (outCube->z - actualCube->z) ); - printf("expect-result= x,y,z = %.11g, %.11g, %.11g \n", - (x - outCube->x), (y - outCube->y), (z - outCube->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); + // psFree(result2); outCube2->x = x; outCube2->y = y; outCube2->z = z; result2 = psCubeToSphere(outCube2); - psFree(result); - result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN); - printf("The apparent output sphere = r,d = %.13g, %.13g\n", apparent->r, apparent->d); + // psFree(result); + psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN); + psFree(result2); + result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + 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); @@ -349,8 +333,8 @@ return 4; } - printf("\n Precession Model output = x,y,s = %.8g, %.8g, %.8g\n", + printf("\n PrecessionModel output = x,y,s = %.13g, %.13g, %.13g\n", pmodel->x, pmodel->y, pmodel->s); - printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - printf(" A difference of: %.13g, %.13g, %.13g\n\n", + printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); + printf(" A difference of: %.13g, %.13g, %.13g\n", (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) ); } @@ -402,10 +386,8 @@ return 7; } else { - printf("\nPrecession Correction output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n\n", + printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", pcorr->x, pcorr->y, pcorr->s); } psFree(pcorr); - - // printf("\n >>MJD value = %lf \n", psTimeToMJD(time2)); //Check values from IERS table B @@ -428,12 +410,14 @@ 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\n", xx, yy, ss); - printf(" A difference of: %.13g, %.13g, %.13g\n\n", + 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; @@ -445,11 +429,6 @@ pcorr->x += xCorr; pcorr->y += yCorr; - psEarthPole *precess = psEOC_PrecessionModel(time2); - precess->x += xCorr; - precess->y += yCorr; psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess); psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv); - // pcorr->x += 3.05224300720406e-7; - // pcorr->y += -1.39441339235822e-7; double q0, q1, q2; q0 = -1.1984522406756289e-5; @@ -461,9 +440,11 @@ printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n"); } - printf(" Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n", - pni->q0, pni->q1, pni->q2, pni->q3); - printf(" Expected sphere rotation = %.13g,%.13g,%.13g\n", q0,q1,q2); - printf(" MY sphere rotation = %.13g,%.13g,%.13g\n", + // printf(" Output from CEOtoGCRS only = %.13g,%.13g,%.13g,%.13g\n", + // pni->q0, pni->q1, pni->q2, pni->q3); + 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; @@ -491,17 +472,17 @@ y = 0.5555012823608123; z = 0.7496183628158023; - 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; + 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); - printf("<x, objC->y, objC->z); - printf(" Difference = %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z); + printf("<x, objC->y, objC->z); + printf(" Difference = %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z); double xx = greatCircle(result, expect); @@ -562,23 +543,21 @@ return 4; } - if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON - || fabs(polarMotion->s - s) > DBL_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "psEOC_GetPolarMotion returned incorrect values.\n"); - printf("\n <>PolarMotion output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", - polarMotion->x, polarMotion->y, polarMotion->s); - printf("\n <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", - x, y, s); - // return 5; - } + // if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + // || fabs(polarMotion->s - s) > DBL_EPSILON) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // "psEOC_GetPolarMotion returned incorrect values.\n"); + printf(" <>PolarMotion output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + // printf(" <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", + // x, y, s); + // return 5; + // } psFree(polarMotion); - //Return valid values for correct input time. Test IERS Bulletin B. polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); - x = -6.45381397904e-07; - y = 2.112819726698e-06; - s = 0.0; - + // x = -6.45381397904e-07; + // y = 2.112819726698e-06; + // s = 0.0; if (polarMotion == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, false, @@ -586,16 +565,14 @@ return 2; } - if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON - || fabs(polarMotion->s - s) > DBL_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "psEOC_GetPolarMotion returned incorrect values.\n"); - printf("\n <>PolarMotion output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", - polarMotion->x, polarMotion->y, polarMotion->s); - printf("\n <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", - x, y, s); - // return 3; - } - - + // if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + // || fabs(polarMotion->s - s) > DBL_EPSILON) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // "psEOC_GetPolarMotion returned incorrect values.\n"); + printf(" <>PolarMotion output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + // printf(" <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + // x, y, s); + // return 3; + // } psEarthPole *nutationCorr = psEOC_NutationCorr(in); if (nutationCorr == NULL) { @@ -603,33 +580,24 @@ return 6; } - polarMotion->x += nutationCorr->x; polarMotion->y += nutationCorr->y; polarMotion->s += nutationCorr->s; + psEarthPole *polarTide = psEOC_PolarTideCorr(in); + polarMotion->x += polarTide->x; + polarMotion->y += polarTide->y; + psFree(polarTide); x = -6.43607313124045e-7; y = 2.11351436973568e-6; s = -7.39617581324646e-12; - - if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON - || fabs(polarMotion->s - s) > DBL_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - " psEOC_GetPolarMotion returned incorrect values.\n"); - printf("\n PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n", - polarMotion->x, polarMotion->y, polarMotion->s); - printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - printf("\n Output Nutation = x,y,s = %.13g, %.13g, %.13g\n", - nutationCorr->x, nutationCorr->y, nutationCorr->s); - // printf(" A difference of: %.13g, %.13g, %.13g\n\n", - // (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) ); - // return 10; - } - - psEarthPole *polarTide = psEOC_PolarTideCorr(in); - polarMotion->x += polarTide->x; - polarMotion->y += polarTide->y; - polarMotion->s += polarTide->s; - psFree(polarTide); - printf("\n PolarMotion + PolarTideCorr out = x,y,s = %.13g, %.13g, %.13g\n", + // if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + // || fabs(polarMotion->s - s) > DBL_EPSILON) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // " psEOC_GetPolarMotion returned incorrect values.\n"); + printf("\n PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n", polarMotion->x, polarMotion->y, polarMotion->s); + printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); + printf(" Difference = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x - x, polarMotion->y - y, polarMotion->s - s); + // } if (!p_psEOCFinalize() ) { @@ -679,6 +647,4 @@ { psTime *in = psTimeAlloc(PS_TIME_UTC); - // in->sec = 1131579114; - // in->nsec = 498489000; in->sec = timesec; in->nsec = 0; @@ -728,7 +694,4 @@ psSphereRot *rot = NULL; psTime *empty = NULL; - // psTime *now = psTimeAlloc(PS_TIME_UT1); - // now->sec = 1128530000; - // now->nsec = 931154510; psTime *time = psTimeAlloc(PS_TIME_UT1); time->sec = timesec-1; @@ -736,8 +699,6 @@ time->leapsecond = false; - // psSphereRot *teoceo = psSphereRot_TEOtoCEO(now); - //return NULL for NULL input time - rot = psSphereRot_TEOtoCEO(empty); + rot = psSphereRot_TEOtoCEO(empty, NULL); if (rot != NULL) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, @@ -746,5 +707,6 @@ } - psSphereRot *teoceo = psSphereRot_TEOtoCEO(time); + psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time); + psSphereRot *teoceo = psSphereRot_TEOtoCEO(time, polarTideCorr); //Make sure values match for other psTime type empty = psTimeAlloc(PS_TIME_UTC); @@ -752,9 +714,14 @@ empty->nsec = 0; empty->leapsecond = false; - rot = psSphereRot_TEOtoCEO(empty); + + rot = psSphereRot_TEOtoCEO(empty, polarTideCorr); if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON || fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n"); + printf("\n Output Rotation1 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n", + teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 ); + printf("\n Output Rotation2 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n", + rot->q0, rot->q1, rot->q2, rot->q3 ); return 2; } @@ -783,5 +750,6 @@ return 3; } - // psFree(now); + + psFree(polarTideCorr); psFree(rot); psFree(empty); @@ -800,20 +768,7 @@ psEarthPole *empty = NULL; psSphereRot *rot = NULL; - - // in->x = M_PI / 4.0; - // in->y = M_PI / 6.0; - // in->s = M_PI / 8.0; - in->x = 2.857175590089105e-4; in->y = 2.3968739377734732e-5; in->s = -1.3970066457904322e-8; - // in->x += 0.06295703125; - // in->y += -0.0287618408203125; - // in->s += 0.0; - // in->x += 3.05224300720406e-7; - // in->y += -1.39441339235822e-7; - in->s += 0.0; - // in->x += -2.491942320903e-10; - // in->y += -1.648366515772e-11; double q0,q1,q2,q3; @@ -837,19 +792,15 @@ } - if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON || - fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) { + printf("\n Output sphere rotation = %.13g, %.13g, %.13g, %.13g\n", + rot->q0, rot->q1, rot->q2, rot->q3); + printf(" Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3); + printf(" difference: %.13g, %.13g, %.13g \n", + (rot->q0-q0), (rot->q1-q1), (rot->q2-q2) ); + if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON || + fabs(rot->q2-q2) > FLT_EPSILON || fabs(rot->q3+q3) > FLT_EPSILON) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psSphereRot_CEOtoGCRS failed to return expected values.\n"); - printf("\n Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n", - rot->q0, rot->q1, rot->q2, rot->q3); - printf(" Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3); - printf(" a difference: %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1), - (rot->q2-q2), (rot->q3-q3) ); - printf(" abs difference: %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1), - (rot->q2+q2), (rot->q3+q3) ); - //return 3; - } - - + return 3; + } psCube *tempCube = psCubeAlloc(); tempCube->x = -0.35963388069046304; @@ -858,7 +809,4 @@ obj = psCubeToSphere(tempCube); psFree(tempCube); - // obj = psSphereAlloc(); - // obj->r = objR; - // obj->d = objD; psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot); psSphere *result = psSphereRotApply(NULL, precessionNutation, obj); @@ -868,11 +816,8 @@ y = 0.5555012823608123; z = 0.7496183628158023; - printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n", cube->x, cube->y, cube->z); - printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z); - printf("A difference of: %.13g, %.13g, %.13g\n\n", + printf("\n Output cube = x,y,z = %.13g, %.13g, %.13g\n", cube->x, cube->y, cube->z); + printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z); + printf(" A difference of: %.13g, %.13g, %.13g\n\n", (x-cube->x), (y-cube->y), (z-cube->z)); - - //great sphere difference between result and expected - //Re = 6.38x10^6m psCube *expect = psCubeAlloc(); expect->x = x; @@ -880,26 +825,16 @@ expect->z = z; psSphere *expected = psCubeToSphere(expect); + double d = greatCircle(result, expected); + printf(" The great circle angular distance between expected & result = %.13g\n", d); + psFree(expect); - double d = acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r)); - printf("Great circle difference of: %.13g \n", d); - psFree(expected); psFree(obj); - - - - // psSphere *test = psSphereAlloc(); - // test->r = 0.0; - // test->d = 0.0; - // test = psSphereRotApply(test, rot, test); - // printf("\n Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d); - psFree(precessionNutation); psFree(result); psFree(cube); - - // psFree(test); psFree(rot); psFree(in); + return 0; } @@ -910,12 +845,7 @@ psEarthPole *empty = NULL; psSphereRot *rot = NULL; - - // in->x = M_PI / 4.0; - // in->y = M_PI / 6.0; - // in->s = M_PI / 8.0; in->x = -0.13275353774074533; in->y = 0.4359436319739848; in->s = -4.2376965863576153e-10; - in->x = SEC_TO_RAD(in->x); in->y = SEC_TO_RAD(in->y); @@ -928,5 +858,4 @@ return 1; } - rot = psSphereRot_ITRStoTEO(in); if (rot == NULL) { @@ -958,7 +887,4 @@ printf(" Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3); - // psSphere *test = psSphereAlloc(); - // test->r = 0.0; - // test->d = 0.0; psCube *temp = psCubeAlloc(); // temp->x = -0.3596195125758298; @@ -969,5 +895,4 @@ temp->z = 0.7496183628158023; - obj = psCubeToSphere(temp); psSphere *test = NULL; @@ -975,5 +900,4 @@ psSphereRot *newRot = psSphereRotConjugate(NULL, rot); test = psSphereRotApply(NULL, newRot, obj); - printf("\n Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d); temp = psSphereToCube(test); printf("\n Cube -test- has x,y,z = %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z); @@ -988,8 +912,7 @@ psSphere *sphere = psCubeToSphere(temp); double d = greatCircle(sphere, test); - printf("Great circle difference of: %.13g \n", d); + psFree(sphere); - psFree(newRot); psFree(obj); @@ -998,5 +921,6 @@ psFree(rot); psFree(in); - return 0; -} - + + return 0; +} +