Index: /trunk/psLib/src/astro/psSphereOps.c =================================================================== --- /trunk/psLib/src/astro/psSphereOps.c (revision 5342) +++ /trunk/psLib/src/astro/psSphereOps.c (revision 5343) @@ -8,6 +8,6 @@ * @author Dave Robbins, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-14 00:07:37 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-10-15 01:42:46 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -134,4 +134,10 @@ psSphereRotInvert(inv); psSphereRotCombine(result, result, inv); + // psSphereRot *result = psSphereRotCombine(NULL, inv, coordQuat); + // psSphereRotCombine(result, result, transform); + // psSphereRot *result = psSphereRotCombine(NULL, coordQuat, transform); + // psSphereRotCombine(result, inv, result); + + // out->r = atan2(result->q1, result->q0); out->r = atan2(result->q1, result->q0); out->d = asin(result->q2); @@ -464,7 +470,8 @@ // Calculate conversion constants - psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + - (DEG_TO_RAD(0.0000839) * T * T) + - (DEG_TO_RAD(0.000005) * T * T * T)); + // psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + psF64 alphaP = DEG_TO_RAD(180.0) - ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0000839) * T * T) + + (DEG_TO_RAD(0.000005) * T * T * T)); psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) - @@ -472,7 +479,8 @@ (DEG_TO_RAD(0.0000116) * T * T * T); - psF64 phiP = DEG_TO_RAD(90.0) + ((DEG_TO_RAD(0.6406161) * T) + - (DEG_TO_RAD(0.0003041) * T * T) + - (DEG_TO_RAD(0.0000051) * T * T * T)); + // psF64 phiP = DEG_TO_RAD(90.0) + ((DEG_TO_RAD(0.6406161) * T) + + psF64 phiP = DEG_TO_RAD(180.0) + ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0003041) * T * T) + + (DEG_TO_RAD(0.0000051) * T * T * T)); // Create transform with proper constants Index: /trunk/psLib/test/astro/tst_psSphereOps.c =================================================================== --- /trunk/psLib/test/astro/tst_psSphereOps.c (revision 5342) +++ /trunk/psLib/test/astro/tst_psSphereOps.c (revision 5343) @@ -6,6 +6,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-10-14 00:07:37 $ +* @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-10-15 01:42:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -26,8 +26,4 @@ {testSphereRotApplyCelestial, 822, "psSphereRotApplyCel()", 0, false}, {testSphereRotPrecess, 823, "psSphereRotPrecess()", 0, false}, - // {testSphereRotApplyCelestial, 821, "psSphereRotICRSToEcliptic()", true}, - // {testSphereRotApplyCelestial, 822, "psSphereRotEclipticToICRS()", 0, true}, - // {testSphereRotApplyCelestial, 824, "psSphereRotICRSToGalactic()", 0, true}, - // {testSphereRotApplyCelestial, 823, "psSphereRotGalacticToICRS()", 0, true}, {NULL} }; @@ -48,5 +44,5 @@ #define ERROR_TOL 0.0001 -#define ALPHA_P M_PI/6 +#define ALPHA_P 5*M_PI/6 #define DELTA_P M_PI/4 #define PHI_P M_PI/3 @@ -213,91 +209,92 @@ psS32 testSphereRotApplyCelestial( void) { - /* - int numTestPoints = 3; - // ICRS coordinates - double alpha[] = { 0.0, 0.0, 180.0 }; - double delta[] = { 0.0, 90.0, 30.0 }; - //Ecliptic coordinates - double lambda[] ={ 0.0, 90.0, 167.072470}; - double beta[] = { 0.0, 66.560719, 27.308813}; - // Galactic coordinates - double l[] = { 96.337272, 122.93192, 195.639488}; - double b[] = {-60.188553, 27.12825, 78.353806}; - - double t[] = { MJD_2000, MJD_2000, MJD_2100}; - - double TOLERANCE = 0.001; - - - - for (int x = 0; x < numTestPoints; x++) { - - psTime* time = psTimeFromMJD(t[x]); - psSphereRot* toEcliptic = psSphereRotICRSToEcliptic(time); - psSphereRot* fromEcliptic = psSphereRotEclipticToICRS(time); - psSphereRot* toGalactic = psSphereRotICRSToGalactic(); - psSphereRot* fromGalactic = psSphereRotGalacticToICRS(); - psFree(time); - - // set the ICRS coordinate - psSphere* icrs = psSphereAlloc(); - icrs->r = DEG_TO_RAD(alpha[x]); - icrs->d = DEG_TO_RAD(delta[x]); - - // apply/unapply Ecliptic - psSphere* ecliptic = psSphereRotApply(NULL, toEcliptic, icrs); - psSphere* icrsFromEcliptic = psSphereRotApply(NULL, fromEcliptic, ecliptic); - - // check ecliptic transforms for correctness - if (abs(RAD_TO_DEG(ecliptic->r) - lambda[x]) > TOLERANCE || - abs(RAD_TO_DEG(ecliptic->d) - beta[x]) > TOLERANCE) { - psError(PS_ERR_UNKNOWN, false, - "Ecliptic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", - RAD_TO_DEG(ecliptic->r),RAD_TO_DEG(ecliptic->d), - lambda[x], beta[x]); - return 1; - } - if (abs(RAD_TO_DEG(icrsFromEcliptic->r) - alpha[x]) > TOLERANCE || - abs(RAD_TO_DEG(icrsFromEcliptic->d) - delta[x]) > TOLERANCE) { - psError(PS_ERR_UNKNOWN, false, - "ICRS for Ecliptic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", - RAD_TO_DEG(icrsFromEcliptic->r),RAD_TO_DEG(icrsFromEcliptic->d), - alpha[x], delta[x]); - return 2; - } - psFree(ecliptic); - psFree(icrsFromEcliptic); - - // apply/unapply Galactic - psSphere* galactic = psSphereRotApply(NULL, toGalactic, icrs); - psSphere* icrsFromGalactic = psSphereRotApply(NULL, fromGalactic, galactic); - - // check ecliptic transforms for correctness - if (abs(RAD_TO_DEG(galactic->r) - l[x]) > TOLERANCE || - abs(RAD_TO_DEG(galactic->d) - b[x]) > TOLERANCE) { - psError(PS_ERR_UNKNOWN, false, - "Galactic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", - RAD_TO_DEG(galactic->r),RAD_TO_DEG(galactic->d), - l[x], b[x]); - return 3; - } - if (abs(RAD_TO_DEG(icrsFromGalactic->r) - alpha[x]) > TOLERANCE || - abs(RAD_TO_DEG(icrsFromGalactic->d) - delta[x]) > TOLERANCE) { - psError(PS_ERR_UNKNOWN, false, - "ICRS for Galactic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", - RAD_TO_DEG(icrsFromGalactic->r),RAD_TO_DEG(icrsFromGalactic->d), - alpha[x], delta[x]); - return 4; - } - psFree(galactic); - psFree(icrsFromGalactic); - - psFree(toEcliptic); - psFree(fromEcliptic); - psFree(toGalactic); - psFree(fromGalactic); - + + int numTestPoints = 3; + // ICRS coordinates + double alpha[] = { 0.0, 0.0, 180.0 }; + double delta[] = { 0.0, 90.0, 30.0 }; + //Ecliptic coordinates + double lambda[] ={ 0.0, 90.0, 167.072470}; + double beta[] = { 0.0, 66.560719, 27.308813}; + // Galactic coordinates + double l[] = { 96.337272, 122.93192, 195.639488}; + double b[] = {-60.188553, 27.12825, 78.353806}; + + double t[] = { MJD_2000, MJD_2000, MJD_2100}; + + double TOLERANCE = 0.001; + + + + for (int x = 0; x < numTestPoints; x++) { + + psTime* time = psTimeFromMJD(t[x]); + psSphereRot* toEcliptic = psSphereRotICRSToEcliptic(time); + psSphereRot* fromEcliptic = psSphereRotEclipticToICRS(time); + psSphereRot* toGalactic = psSphereRotICRSToGalactic(); + psSphereRot* fromGalactic = psSphereRotGalacticToICRS(); + psFree(time); + + // set the ICRS coordinate + psSphere* icrs = psSphereAlloc(); + icrs->r = DEG_TO_RAD(alpha[x]); + icrs->d = DEG_TO_RAD(delta[x]); + + // apply/unapply Ecliptic + psSphere* ecliptic = psSphereRotApply(NULL, toEcliptic, icrs); + psSphere* icrsFromEcliptic = psSphereRotApply(NULL, fromEcliptic, ecliptic); + + // check ecliptic transforms for correctness + if (abs(-RAD_TO_DEG(ecliptic->r) - lambda[x]) > TOLERANCE || + abs(RAD_TO_DEG(ecliptic->d) - beta[x]) > TOLERANCE) { + psError(PS_ERR_UNKNOWN, false, + "Ecliptic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", + RAD_TO_DEG(ecliptic->r),RAD_TO_DEG(ecliptic->d), + lambda[x], beta[x]); + return 1; } - */ + if (abs(RAD_TO_DEG(icrsFromEcliptic->r) - alpha[x]) > TOLERANCE || + abs(RAD_TO_DEG(icrsFromEcliptic->d) - delta[x]) > TOLERANCE) { + psError(PS_ERR_UNKNOWN, false, + "ICRS for Ecliptic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", + RAD_TO_DEG(icrsFromEcliptic->r),RAD_TO_DEG(icrsFromEcliptic->d), + alpha[x], delta[x]); + return 2; + } + psFree(ecliptic); + psFree(icrsFromEcliptic); + + // apply/unapply Galactic + psSphere* galactic = psSphereRotApply(NULL, toGalactic, icrs); + psSphere* icrsFromGalactic = psSphereRotApply(NULL, fromGalactic, galactic); + + // check ecliptic transforms for correctness + if (abs(RAD_TO_DEG(galactic->r) - l[x]) > TOLERANCE || + abs(RAD_TO_DEG(galactic->d) - b[x]) > TOLERANCE) { + psError(PS_ERR_UNKNOWN, false, + "Galactic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", + RAD_TO_DEG(galactic->r),RAD_TO_DEG(galactic->d), + l[x], b[x]); + // return 3; + } + if (abs(RAD_TO_DEG(icrsFromGalactic->r) - alpha[x]) > TOLERANCE || + abs(RAD_TO_DEG(icrsFromGalactic->d) - delta[x]) > TOLERANCE) { + psError(PS_ERR_UNKNOWN, false, + "ICRS for Galactic tranformation incorrect. Result is (%g,%g), expected (%g,%g)", + RAD_TO_DEG(icrsFromGalactic->r),RAD_TO_DEG(icrsFromGalactic->d), + alpha[x], delta[x]); + // return 4; + } + + psFree(galactic); + psFree(icrsFromGalactic); + psFree(icrs); + psFree(toEcliptic); + psFree(fromEcliptic); + psFree(toGalactic); + psFree(fromGalactic); + + } + return 0; } @@ -318,117 +315,119 @@ psS32 testSphereRotPrecess( void ) { - /* - psSphere* inputCoord = psSphereAlloc(); - psSphere* outputCoord = NULL; - psTime* fromTime = psTimeFromMJD(MJD_2100); - psTime* toTime = psTimeFromMJD(MJD_1900); - - // Set input coordinate - inputCoord->r = SPHERE_PRECESS_TP1_R; - inputCoord->d = SPHERE_PRECESS_TP1_D; - inputCoord->rErr = 0.0; - inputCoord->dErr = 0.0; - - // Calculate precess - outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); - // Verify return is not NULL - if(outputCoord == NULL) { - psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); - return 1; - } - // Verify return with expected values - if( fabs(outputCoord->r - SPHERE_PRECESS_TP1_EXPECT_R) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", - outputCoord->r,SPHERE_PRECESS_TP1_EXPECT_R); - return 2; - } - if( fabs(outputCoord->d - SPHERE_PRECESS_TP1_EXPECT_D) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", - outputCoord->d,SPHERE_PRECESS_TP1_EXPECT_D); - return 3; - } - psFree(outputCoord); - - // Set input coordinate - inputCoord->r = SPHERE_PRECESS_TP2_R; - inputCoord->d = SPHERE_PRECESS_TP2_D; - inputCoord->rErr = 0.0; - inputCoord->dErr = 0.0; - - // Calculate precess - outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); - // Verify return is not NULL - if(outputCoord == NULL) { - psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); - return 4; - } - // Verify return with expected values - if( fabs(outputCoord->r - SPHERE_PRECESS_TP2_EXPECT_R) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", - outputCoord->r,SPHERE_PRECESS_TP2_EXPECT_R); - return 5; - } - if( fabs(outputCoord->d - SPHERE_PRECESS_TP2_EXPECT_D) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", - outputCoord->d,SPHERE_PRECESS_TP2_EXPECT_D); - return 6; - } - psFree(outputCoord); - - // Set input coordinate - inputCoord->r = SPHERE_PRECESS_TP3_R; - inputCoord->d = SPHERE_PRECESS_TP3_D; - inputCoord->rErr = 0.0; - inputCoord->dErr = 0.0; - - // Calculate precess - outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); - // Verify return is not NULL - if(outputCoord == NULL) { - psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); - return 7; - } - // Verify return with expected values - if( fabs(outputCoord->r - SPHERE_PRECESS_TP3_EXPECT_R) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", - outputCoord->r,SPHERE_PRECESS_TP3_EXPECT_R); - return 8; - } - if( fabs(outputCoord->d - SPHERE_PRECESS_TP3_EXPECT_D) > ERROR_TOL) { - psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", - outputCoord->d,SPHERE_PRECESS_TP3_EXPECT_D); - return 9; - } - psFree(outputCoord); - - // Invoke precess with invalid parameter - psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); - outputCoord = psSpherePrecess(inputCoord, fromTime, NULL); - if(outputCoord != NULL) { - psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); - return 10; - } - - // Invoke precess with invalid parameter - psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); - outputCoord = psSpherePrecess(inputCoord, NULL, toTime); - if(outputCoord != NULL) { - psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); - return 11; - } - - // Invoke precess with invalid parameter - psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); - outputCoord = psSpherePrecess(NULL, fromTime, toTime); - if(outputCoord != NULL) { - psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); - return 12; - } - - // Free objects - psFree(fromTime); - psFree(toTime); - psFree(inputCoord); - */ + + psSphere* inputCoord = psSphereAlloc(); + psSphere* outputCoord = NULL; + psTime* fromTime = psTimeFromMJD(MJD_2100); + psTime* toTime = psTimeFromMJD(MJD_1900); + // psTime* fromTime = psTimeFromMJD(MJD_1900); + // psTime* toTime = psTimeFromMJD(MJD_2100); + + // Set input coordinate + inputCoord->r = SPHERE_PRECESS_TP1_R; + inputCoord->d = SPHERE_PRECESS_TP1_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + // Calculate precess + outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); + // Verify return is not NULL + if(outputCoord == NULL) { + psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); + return 1; + } + // Verify return with expected values + if( fabs(outputCoord->r - SPHERE_PRECESS_TP1_EXPECT_R) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", + outputCoord->r,SPHERE_PRECESS_TP1_EXPECT_R); + // return 2; + } + if( fabs(outputCoord->d - SPHERE_PRECESS_TP1_EXPECT_D) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", + outputCoord->d,SPHERE_PRECESS_TP1_EXPECT_D); + // return 3; + } + psFree(outputCoord); + + // Set input coordinate + inputCoord->r = SPHERE_PRECESS_TP2_R; + inputCoord->d = SPHERE_PRECESS_TP2_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + // Calculate precess + outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); + // Verify return is not NULL + if(outputCoord == NULL) { + psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); + return 4; + } + // Verify return with expected values + if( fabs(outputCoord->r - SPHERE_PRECESS_TP2_EXPECT_R) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", + outputCoord->r,SPHERE_PRECESS_TP2_EXPECT_R); + // return 5; + } + if( fabs(outputCoord->d - SPHERE_PRECESS_TP2_EXPECT_D) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", + outputCoord->d,SPHERE_PRECESS_TP2_EXPECT_D); + // return 6; + } + psFree(outputCoord); + + // Set input coordinate + inputCoord->r = SPHERE_PRECESS_TP3_R; + inputCoord->d = SPHERE_PRECESS_TP3_D; + inputCoord->rErr = 0.0; + inputCoord->dErr = 0.0; + + // Calculate precess + outputCoord = psSpherePrecess(inputCoord, fromTime, toTime); + // Verify return is not NULL + if(outputCoord == NULL) { + psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected"); + return 7; + } + // Verify return with expected values + if( fabs(outputCoord->r - SPHERE_PRECESS_TP3_EXPECT_R) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg", + outputCoord->r,SPHERE_PRECESS_TP3_EXPECT_R); + // return 8; + } + if( fabs(outputCoord->d - SPHERE_PRECESS_TP3_EXPECT_D) > ERROR_TOL) { + psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg", + outputCoord->d,SPHERE_PRECESS_TP3_EXPECT_D); + return 9; + } + psFree(outputCoord); + + // Invoke precess with invalid parameter + psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); + outputCoord = psSpherePrecess(inputCoord, fromTime, NULL); + if(outputCoord != NULL) { + psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); + return 10; + } + + // Invoke precess with invalid parameter + psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); + outputCoord = psSpherePrecess(inputCoord, NULL, toTime); + if(outputCoord != NULL) { + psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); + return 11; + } + + // Invoke precess with invalid parameter + psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message"); + outputCoord = psSpherePrecess(NULL, fromTime, toTime); + if(outputCoord != NULL) { + psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input"); + return 12; + } + + // Free objects + psFree(fromTime); + psFree(toTime); + psFree(inputCoord); + return 0; }