Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 2982) +++ /trunk/psLib/src/astro/psCoord.c (revision 2983) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-10 19:47:10 $ +* @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-01-13 22:45:28 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,4 +24,5 @@ #include "psConstants.h" #include "psError.h" +#include "psLogMsg.h" #include "psAstronomyErrors.h" #include @@ -183,4 +184,23 @@ /****************************************************************************** +XXX: Private Function. + +piNormalize(): take an input angle in radians and convert it to the range 0:PI. + *****************************************************************************/ +psF32 piNormalize(psF32 angle) +{ + while (angle < FLT_EPSILON) { + angle+=M_PI; + } + + while (angle >= M_PI) { + angle-=M_PI; + } + return(angle); +} + +/****************************************************************************** +XXX: We convert Right Ascension angles to the range 0:PI. Is that acceptable? +XXX: Should we do something for Declination as well? *****************************************************************************/ psSphere* psSphereTransformApply(psSphere* out, @@ -210,5 +230,6 @@ psF64 phi = atan2(eq56, eq57) + transform->alphaP; - out->r = phi; + + out->r = piNormalize(phi); out->d = theta; @@ -281,4 +302,5 @@ XXX: Waiting for the definition of the PS_PROJ_PAR projection. XXX: Waiting for the definition of the PS_PROJ_GLS projection. +XXX: Must apply scaling at the end. *****************************************************************************/ psPlane* psProject(const psSphere* coord, @@ -288,13 +310,44 @@ PS_PTR_CHECK_NULL(projection, NULL); - // float R = 0.0; - // float alpha = 0.0; psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane)); - if (projection->type == PS_PROJ_TAN) { - tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r); - tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r); + + /* + // ******************************************** + // From the ADD + // ******************************************** + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); + // return(NULL); + psF32 sinTheta; + psF32 cosThetaCosPhi; + psF32 cosThetaSinPhi; + + sinTheta = (sin(coord->d) * sin(projection->D)) + + (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R))); + cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) - + (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R))); + cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R); + + tmp->x = -cosThetaSinPhi / sinTheta; + tmp->y = cosThetaCosPhi / sinTheta; + + */ + + // ******************************************** + // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html + // delta_0 and phi_1 are the projection centers, not the point being projected. + // (delta_0, phi_1) == (projection->R, projection->D) + // ******************************************** + psF32 cosC = (sin(projection->D) * sin(coord->d)) + + (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R)); + tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC; + tmp->y = ((cos(projection->D) * sin(coord->d)) - + (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) / + cosC; } else if (projection->type == PS_PROJ_SIN) { + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); + // return(NULL); + tmp->x = cos(coord->d) * sin(coord->d); tmp->y = -cos(coord->d) * cos(coord->d); @@ -332,6 +385,8 @@ /****************************************************************************** +XXX: Waiting for algorithms. XXX: Waiting for the definition of the PS_PROJ_PAR projection. XXX: Waiting for the definition of the PS_PROJ_GLS projection. +XXX: Must apply scaling at the beginning. *****************************************************************************/ psSphere* psDeproject(const psPlane* coord, @@ -340,5 +395,6 @@ PS_PTR_CHECK_NULL(coord, NULL); PS_PTR_CHECK_NULL(projection, NULL); - + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work."); + // return(NULL); float R = 0.0; float chu = 0.0; @@ -348,7 +404,39 @@ if (projection->type == PS_PROJ_TAN) { - tmp->r = atan2(-coord->y, coord->x); - R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = atan(1.0 / R); + /* + // ******************************************** + // From the ADD + // ******************************************** + psF32 phi; + psF32 theta; + psF32 sinDelta; + psF32 cosDeltaCosAlphaMinusAlphaP; + psF32 cosDeltaSinAlphaMinusAlphaP; + + phi = atan2(-coord->y, coord->x); + R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); + theta = atan(1.0 / R); + + sinDelta = (sin(theta) * sin(projection->D)) + + (cos(theta) * cos(projection->D) * cos(phi)); + cosDeltaCosAlphaMinusAlphaP = + (sin(theta) * cos(projection->D)) - + (cos(theta) * sin(projection->D) * cos(phi)); + cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi); + tmp->d = asin(sinDelta); + tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R; + */ + // ******************************************** + // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html + // delta_0 and phi_1 are the projection centers, not the point being projected. + // (delta_0, phi_1) == (projection->R, projection->D) + // ******************************************** + psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); + psF32 C = atan(row); + tmp->d = asin((cos(C) * sin(projection->D)) + ((coord->y * sin(C) * cos(projection->D)) / row)); + psF32 tmpAtan = atan((coord->x * sin(C)) / + ((row * cos(projection->D) * cos(C)) - (coord->y * sin(projection->D) * sin(C)))); + tmp->r = projection->R + tmpAtan; + } else if (projection->type == PS_PROJ_SIN) { @@ -415,4 +503,13 @@ PS_PTR_CHECK_NULL(position2, NULL); + if (position1->d >= DEG_TO_RAD(90.0)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees. Returning NULL.\n"); + return(NULL); + } + if (position2->d >= DEG_TO_RAD(90.0)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees. Returning NULL.\n"); + return(NULL); + } + psPlane* lin; psProjection proj; @@ -426,15 +523,8 @@ proj.type = PS_PROJ_TAN; - printf("------------ position 1 is (%f, %f)\n", position1->r, position1->d); - printf("------------ position 2 is (%f, %f)\n", position2->r, position2->d); lin = psProject(position1, &proj); - printf("------------ projected position 1 is (%f, %f)\n", lin->y, lin->x); lin = psProject(position2, &proj); - printf("------------ projected position 2 is (%f, %f)\n", lin->y, lin->x); tmp = psDeproject(lin, &proj); - printf("------------ deprojected position 2 is (%f, %f)\n", tmp->r, tmp->d); - psFree(lin); - // XXX: Do we need to convert units in tmp? return (tmp); Index: /trunk/psLib/src/astronomy/psAstrometry.c =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.c (revision 2982) +++ /trunk/psLib/src/astronomy/psAstrometry.c (revision 2983) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.54 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-01-07 23:52:15 $ + * @version $Revision: 1.55 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-01-13 22:45:28 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -566,5 +566,5 @@ } } - + printf("Hmmm, returning NULL\n"); return (NULL); } Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 2982) +++ /trunk/psLib/src/astronomy/psCoord.c (revision 2983) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-10 19:47:10 $ +* @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-01-13 22:45:28 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,4 +24,5 @@ #include "psConstants.h" #include "psError.h" +#include "psLogMsg.h" #include "psAstronomyErrors.h" #include @@ -183,4 +184,23 @@ /****************************************************************************** +XXX: Private Function. + +piNormalize(): take an input angle in radians and convert it to the range 0:PI. + *****************************************************************************/ +psF32 piNormalize(psF32 angle) +{ + while (angle < FLT_EPSILON) { + angle+=M_PI; + } + + while (angle >= M_PI) { + angle-=M_PI; + } + return(angle); +} + +/****************************************************************************** +XXX: We convert Right Ascension angles to the range 0:PI. Is that acceptable? +XXX: Should we do something for Declination as well? *****************************************************************************/ psSphere* psSphereTransformApply(psSphere* out, @@ -210,5 +230,6 @@ psF64 phi = atan2(eq56, eq57) + transform->alphaP; - out->r = phi; + + out->r = piNormalize(phi); out->d = theta; @@ -281,4 +302,5 @@ XXX: Waiting for the definition of the PS_PROJ_PAR projection. XXX: Waiting for the definition of the PS_PROJ_GLS projection. +XXX: Must apply scaling at the end. *****************************************************************************/ psPlane* psProject(const psSphere* coord, @@ -288,13 +310,44 @@ PS_PTR_CHECK_NULL(projection, NULL); - // float R = 0.0; - // float alpha = 0.0; psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane)); - if (projection->type == PS_PROJ_TAN) { - tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r); - tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r); + + /* + // ******************************************** + // From the ADD + // ******************************************** + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); + // return(NULL); + psF32 sinTheta; + psF32 cosThetaCosPhi; + psF32 cosThetaSinPhi; + + sinTheta = (sin(coord->d) * sin(projection->D)) + + (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R))); + cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) - + (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R))); + cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R); + + tmp->x = -cosThetaSinPhi / sinTheta; + tmp->y = cosThetaCosPhi / sinTheta; + + */ + + // ******************************************** + // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html + // delta_0 and phi_1 are the projection centers, not the point being projected. + // (delta_0, phi_1) == (projection->R, projection->D) + // ******************************************** + psF32 cosC = (sin(projection->D) * sin(coord->d)) + + (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R)); + tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC; + tmp->y = ((cos(projection->D) * sin(coord->d)) - + (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) / + cosC; } else if (projection->type == PS_PROJ_SIN) { + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); + // return(NULL); + tmp->x = cos(coord->d) * sin(coord->d); tmp->y = -cos(coord->d) * cos(coord->d); @@ -332,6 +385,8 @@ /****************************************************************************** +XXX: Waiting for algorithms. XXX: Waiting for the definition of the PS_PROJ_PAR projection. XXX: Waiting for the definition of the PS_PROJ_GLS projection. +XXX: Must apply scaling at the beginning. *****************************************************************************/ psSphere* psDeproject(const psPlane* coord, @@ -340,5 +395,6 @@ PS_PTR_CHECK_NULL(coord, NULL); PS_PTR_CHECK_NULL(projection, NULL); - + // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work."); + // return(NULL); float R = 0.0; float chu = 0.0; @@ -348,7 +404,39 @@ if (projection->type == PS_PROJ_TAN) { - tmp->r = atan2(-coord->y, coord->x); - R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = atan(1.0 / R); + /* + // ******************************************** + // From the ADD + // ******************************************** + psF32 phi; + psF32 theta; + psF32 sinDelta; + psF32 cosDeltaCosAlphaMinusAlphaP; + psF32 cosDeltaSinAlphaMinusAlphaP; + + phi = atan2(-coord->y, coord->x); + R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); + theta = atan(1.0 / R); + + sinDelta = (sin(theta) * sin(projection->D)) + + (cos(theta) * cos(projection->D) * cos(phi)); + cosDeltaCosAlphaMinusAlphaP = + (sin(theta) * cos(projection->D)) - + (cos(theta) * sin(projection->D) * cos(phi)); + cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi); + tmp->d = asin(sinDelta); + tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R; + */ + // ******************************************** + // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html + // delta_0 and phi_1 are the projection centers, not the point being projected. + // (delta_0, phi_1) == (projection->R, projection->D) + // ******************************************** + psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); + psF32 C = atan(row); + tmp->d = asin((cos(C) * sin(projection->D)) + ((coord->y * sin(C) * cos(projection->D)) / row)); + psF32 tmpAtan = atan((coord->x * sin(C)) / + ((row * cos(projection->D) * cos(C)) - (coord->y * sin(projection->D) * sin(C)))); + tmp->r = projection->R + tmpAtan; + } else if (projection->type == PS_PROJ_SIN) { @@ -415,4 +503,13 @@ PS_PTR_CHECK_NULL(position2, NULL); + if (position1->d >= DEG_TO_RAD(90.0)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees. Returning NULL.\n"); + return(NULL); + } + if (position2->d >= DEG_TO_RAD(90.0)) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees. Returning NULL.\n"); + return(NULL); + } + psPlane* lin; psProjection proj; @@ -426,15 +523,8 @@ proj.type = PS_PROJ_TAN; - printf("------------ position 1 is (%f, %f)\n", position1->r, position1->d); - printf("------------ position 2 is (%f, %f)\n", position2->r, position2->d); lin = psProject(position1, &proj); - printf("------------ projected position 1 is (%f, %f)\n", lin->y, lin->x); lin = psProject(position2, &proj); - printf("------------ projected position 2 is (%f, %f)\n", lin->y, lin->x); tmp = psDeproject(lin, &proj); - printf("------------ deprojected position 2 is (%f, %f)\n", tmp->r, tmp->d); - psFree(lin); - // XXX: Do we need to convert units in tmp? return (tmp); Index: /trunk/psLib/test/astronomy/tst_psAstrometry01.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 2982) +++ /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 2983) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-11-30 21:56:20 $ +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-01-13 22:45:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -21,5 +21,5 @@ {test1, 0000, "Tests psFixedPatternAlloc()", 0, false}, {test2, 0000, "Tests psGrommitAlloc()", 0, false}, - {test3, 0000, "psFunctionBar", 0, false}, + {test3, 0000, "MISC", 0, false}, {NULL} }; @@ -189,5 +189,5 @@ psLogSetLevel( PS_LOG_INFO ); - return ( ! runTestSuite( stderr, "psImage", tests, argc, argv ) ); + return ( ! runTestSuite( stderr, "psAstrometry", tests, argc, argv ) ); } @@ -440,5 +440,10 @@ if ((myCell == NULL) || (tmpChip == NULL)) { - printf("ERROR: NULL\n"); + if (myCell == NULL) { + printf("ERROR: myCell is NULL\n"); + } + if (tmpChip == NULL) { + printf("ERROR: tmpChip is NULL\n"); + } } else { psCoordFPAToChip(&chipCoord, &fpaCoord, tmpChip); @@ -517,5 +522,7 @@ } + printf("HERE 00\n"); psFree(myFPA); + printf("HERE 01\n"); return(0); Index: /trunk/psLib/test/astronomy/tst_psCoord.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psCoord.c (revision 2982) +++ /trunk/psLib/test/astronomy/tst_psCoord.c (revision 2983) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-01-08 00:12:44 $ +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-01-13 22:45:40 $ * * XXX: Must test psSpherePrecess. @@ -18,4 +18,5 @@ #include "psTest.h" #include "pslib.h" +/* static psS32 test1( void ); static psS32 test1b( void ); @@ -24,4 +25,5 @@ static psS32 test3( void ); static psS32 test4( void ); +static psS32 test4b( void ); static psS32 test5( void ); static psS32 test6( void ); @@ -32,19 +34,27 @@ static psS32 test40( void ); static psS32 test41( void ); +static psS32 test43( void ); +*/ +static psS32 test42( void ); testDescription tests[] = { - {test1, 0000, "psSphereTransformAlloc()", 0, false}, - {test1b, 0000, "psPlaneTransformAlloc()", 0, false}, - {test1c, 0000, "psPlaneDistortAlloc()", 0, false}, - {test2, 0000, "psPlaneTransformApply()", 0, false}, - {test3, 0000, "psPlaneDistortApply()", 0, false}, - {test4, 0000, "psPSphereTransformApply()", 0, false}, - {test5, 0000, "psSphereTransformICRSToEcliptic()", 0, false}, - {test6, 0000, "psSphereTransformEclipticToICRS()", 0, false}, - {test7, 0000, "psSphereTransformICRSToGalatic()", 0, false}, - {test8, 0000, "psSphereTransformGalaticToICRS()", 0, false}, - {test20, 0000, "psProject()", 0, false}, - {test21, 0000, "psDeProject()", 0, false}, - {test40, 0000, "psSphereGetOffset()", 0, false}, - {test41, 0000, "psSphereSetOffset()", 0, false}, + /* + {test1, 0000, "psSphereTransformAlloc()", 0, false}, + {test1b, 0000, "psPlaneTransformAlloc()", 0, false}, + {test1c, 0000, "psPlaneDistortAlloc()", 0, false}, + {test2, 0000, "psPlaneTransformApply()", 0, false}, + {test3, 0000, "psPlaneDistortApply()", 0, false}, + {test4, 0000, "psPSphereTransformApply()", 0, false}, + {test4b, 0000, "psPSphereTransformApply()", 0, false}, + {test5, 0000, "psSphereTransformICRSToEcliptic()", 0, false}, + {test6, 0000, "psSphereTransformEclipticToICRS()", 0, false}, + {test7, 0000, "psSphereTransformICRSToGalatic()", 0, false}, + {test8, 0000, "psSphereTransformGalaticToICRS()", 0, false}, + {test20, 0000, "psProject()", 0, false}, + {test21, 0000, "psDeProject()", 0, false}, + {test40, 0000, "psSphereGetOffset()", 0, false}, + {test41, 0000, "psSphereSetOffset()", 0, false}, + {test43, 0000, "psProject(), psDeproject", 0, false}, + */ + {test42, 0000, "psProject(), psDeproject", 0, false}, {NULL} }; @@ -436,6 +446,6 @@ psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0); - in.r = 45.0 * MY_DEG_TO_RAD; - in.d = 30.0 * MY_DEG_TO_RAD; + in.r = DEG_TO_RAD(45.0); + in.d = DEG_TO_RAD(30.0); in.rErr = 0.0; in.dErr = 0.0; @@ -443,6 +453,6 @@ for (float r=0.0;r<180.0;r+=DEG_INC) { for (float d=0.0;d<90.0;d+=DEG_INC) { - in.r = r * MY_DEG_TO_RAD; - in.d = d * MY_DEG_TO_RAD; + in.r = DEG_TO_RAD(r); + in.d = DEG_TO_RAD(d); in.rErr = 0.0; in.dErr = 0.0; @@ -483,4 +493,88 @@ printf("-------------------------------------------------------------------\n"); psFree(myST); + return(testStatus); +} + +/****************************************************************************** +test4b(): This test verifies that psSphereTransformApply() works properly. We +create two psSphereTransforms: a forward transform and a reverse transform +(which is the mathematical inverse of the forward transform). We apply both +transforms to several spherical coordinates and ensure that the original input +coordinate is obtained after applying both transforms. + +XXX: We currently test the alpha and delta offsets independently. Attempts to +test them both concurrently failed. Determine why this is. Are the following +spherical transforms not mathematical inverses? + psSphereTransformAlloc(X, Y, 0.0) + psSphereTransformAlloc(-X, -Y, 0.0) + *****************************************************************************/ +#define ERROR_PERCENT 0.01 +psS32 test4b( void ) +{ + psS32 testStatus = 0; + psSphere in; + psSphere out; + psSphere out2; + psSphereTransform *mySphereTransformForward = NULL; + psSphereTransform *mySphereTransformReverse = NULL; + + + mySphereTransformForward = psSphereTransformAlloc(DEG_TO_RAD(22.0), + 0.0, + 0.0); + mySphereTransformReverse = psSphereTransformAlloc(DEG_TO_RAD(-22.0), + 0.0, + 0.0); + + for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) { + for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) { + in.r = DEG_TO_RAD(r); + in.d = DEG_TO_RAD(d); + in.rErr = 0.0; + in.dErr = 0.0; + + psSphereTransformApply(&out, mySphereTransformForward, &in); + psSphereTransformApply(&out2, mySphereTransformReverse, &out); + + if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) || + (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) { + printf("ERROR: \n"); + printf("Input coords (R, D) are (%f, %f)\n", in.r, in.d); + printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d); + testStatus = 4; + } + } + } + psFree(mySphereTransformForward); + psFree(mySphereTransformReverse); + + mySphereTransformForward = psSphereTransformAlloc(0.0, + DEG_TO_RAD(33.0), + 0.0); + mySphereTransformReverse = psSphereTransformAlloc(0.0, + DEG_TO_RAD(-33.0), + 0.0); + for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) { + for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) { + in.r = DEG_TO_RAD(r); + in.d = DEG_TO_RAD(d); + in.rErr = 0.0; + in.dErr = 0.0; + + psSphereTransformApply(&out, mySphereTransformForward, &in); + psSphereTransformApply(&out2, mySphereTransformReverse, &out); + + if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) || + (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) { + printf("ERROR: \n"); + printf("Input coords (R, D) are (%f, %f)\n", in.r, in.d); + printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d); + testStatus = 4; + } + } + } + psFree(mySphereTransformForward); + psFree(mySphereTransformReverse); + return(testStatus); } @@ -677,6 +771,6 @@ psPlane *out; psProjection myProjection; - myProjection.R = 20.0 * MY_DEG_TO_RAD; - myProjection.D = 10.0 * MY_DEG_TO_RAD; + myProjection.R = DEG_TO_RAD(20.0); + myProjection.D = DEG_TO_RAD(10.0); myProjection.Xs = 1.0; myProjection.Ys = 1.0; @@ -684,6 +778,6 @@ for (float r=0.0;r<180.0;r+=DEG_INC) { for (float d=0.0;d<90.0;d+=DEG_INC) { - in.r = r * MY_DEG_TO_RAD; - in.d = d * MY_DEG_TO_RAD; + in.r = DEG_TO_RAD(r); + in.d = DEG_TO_RAD(d); in.rErr = 0.0; in.dErr = 0.0; @@ -789,6 +883,6 @@ psSphere *out; psProjection myProjection; - myProjection.R = 20.0 * MY_DEG_TO_RAD; - myProjection.D = 10.0 * MY_DEG_TO_RAD; + myProjection.R = DEG_TO_RAD(20.0); + myProjection.D = DEG_TO_RAD(10.0); myProjection.Xs = 1.0; myProjection.Ys = 1.0; @@ -904,6 +998,6 @@ psSphere *offset = NULL; - position1.r = 90.0 * MY_DEG_TO_RAD; - position1.d = 45.0 * MY_DEG_TO_RAD; + position1.r = DEG_TO_RAD(90.0); + position1.d = DEG_TO_RAD(45.0); position1.rErr = 0.0; position1.dErr = 0.0; @@ -911,6 +1005,6 @@ for (float r=0.0;r<180.0;r+=DEG_INC) { for (float d=0.0;d<90.0;d+=DEG_INC) { - position2.r = r * MY_DEG_TO_RAD; - position2.d = d * MY_DEG_TO_RAD; + position2.r = DEG_TO_RAD(r); + position2.d = DEG_TO_RAD(d); position2.rErr = 0.0; position2.dErr = 0.0; @@ -924,6 +1018,6 @@ offset = psSphereGetOffset(&position1, &position2, PS_SPHERICAL, PS_DEGREE); - offset->r = offset->r * MY_DEG_TO_RAD; - offset->d = offset->d * MY_DEG_TO_RAD; + offset->r = DEG_TO_RAD(offset->r); + offset->d = DEG_TO_RAD(offset->d); PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2); PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3); @@ -946,14 +1040,15 @@ psFree(offset); - /* XXX: This code does not work correctly. - offset = psSphereGetOffset(&position1, &position2, - PS_LINEAR, 0); - printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n", - position1.r, position1.d, position2.r, position2.d, - offset->r, offset->d); + //HEY /* XXX: This code does not work correctly. + printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n", + position1.r, position1.d, position2.r, position2.d, + offset->r, offset->d); + /* + offset = psSphereGetOffset(&position1, &position2, PS_LINEAR, 0); PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1); PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1); - psFree(offset); */ + psFree(offset); + // */ } } @@ -1007,6 +1102,6 @@ psSphere tmpOffset; - position1.r = 90.0 * MY_DEG_TO_RAD; - position1.d = 45.0 * MY_DEG_TO_RAD; + position1.r = DEG_TO_RAD(90.0); + position1.d = DEG_TO_RAD(45.0); position1.rErr = 0.0; position1.dErr = 0.0; @@ -1014,6 +1109,6 @@ for (float r=0.0;r<180.0;r+=DEG_INC) { for (float d=0.0;d<90.0;d+=DEG_INC) { - offset.r = r * MY_DEG_TO_RAD; - offset.d = d * MY_DEG_TO_RAD; + offset.r = DEG_TO_RAD(offset.r); + offset.d = DEG_TO_RAD(offset.d); offset.rErr = 0.0; offset.dErr = 0.0; @@ -1025,6 +1120,6 @@ psFree(position2); - tmpOffset.r = offset.r * MY_RAD_TO_DEG; - tmpOffset.d = offset.d * MY_RAD_TO_DEG; + tmpOffset.r = DEG_TO_RAD(offset.r); + tmpOffset.d = DEG_TO_RAD(offset.d); tmpOffset.rErr = 0.0; tmpOffset.dErr = 0.0; @@ -1108,3 +1203,106 @@ } +/****************************************************************************** +test43(): Not really a test. We simply project many points and print the +results. + *****************************************************************************/ +psS32 test43( void ) +{ + psS32 testStatus = 0; + psSphere projCenter; + psSphere startCoords; + psPlane *planeCoords; + + projCenter.r = DEG_TO_RAD(33.0); + projCenter.d = DEG_TO_RAD(22.0); + projCenter.r = 0.0; + projCenter.d = 0.0; + // Create the psProjection plane: + psProjection myProj; + myProj.R = projCenter.r; + myProj.D = projCenter.d; + myProj.Xs = 1.0; + myProj.Ys = 1.0; + myProj.type = PS_PROJ_TAN; + + printf(" projCenter is (%f, %f)\n", projCenter.r, projCenter.d); + for (float r=-90.0;r<90.0;r+=DEG_INC/5.0) { + for (float d=-90.0;d<90.0;d+=DEG_INC/5.0) { + printf("----------------------------------------------------------\n"); + printf("Angles (R, D) are (%.2f, %.2f)\n", r, d); + startCoords.r = DEG_TO_RAD(r); + startCoords.d = DEG_TO_RAD(d); + + planeCoords = psProject(&startCoords, &myProj); + printf(" startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d); + printf(" plane Coords (X, Y) is (%f, %f)\n", planeCoords->x, planeCoords->y); + psFree(planeCoords); + } + } + printf("----------------------------------------------------------\n"); + return(testStatus); +} + +#define NUM_DEGREES 5 +/****************************************************************************** +test42(): This test attempts to verify the psProject() and psDeproject() +functions. We create an arbitrary psProjection. We then project several +spherical coordinates onto the plane and then immediately deproject them; we +ensure that the resulting spherical coordinates are identical to the start +coordinates. + *****************************************************************************/ +psS32 test42( void ) +{ + psS32 testStatus = 0; + psSphere projCenter; + psSphere startCoords; + psSphere *endCoords; + psPlane *planeCoords; + + #define PROG_CENTER_R_DEG 20.0 + #define PROG_CENTER_D_DEG 20.0 + + projCenter.r = DEG_TO_RAD(PROG_CENTER_R_DEG); + projCenter.d = DEG_TO_RAD(PROG_CENTER_D_DEG); + // Create the psProjection plane: + psProjection myProj; + myProj.R = projCenter.r; + myProj.D = projCenter.d; + myProj.Xs = 1.0; + myProj.Ys = 1.0; + myProj.type = PS_PROJ_TAN; + + #define RADIUS_DEG 30.0 + #define NUM_ANGLES 20 + #define MY_DEG_INC ((2.0 * RADIUS_DEG) / NUM_ANGLES) + + for (float r= (PROG_CENTER_R_DEG - RADIUS_DEG); + r<=(PROG_CENTER_R_DEG + RADIUS_DEG); + r+=MY_DEG_INC) { + printf("------------------------ r is %.1f ------------------------\n", r); + for (float d= (PROG_CENTER_D_DEG - RADIUS_DEG); + d<=(PROG_CENTER_D_DEG + RADIUS_DEG); + d+=MY_DEG_INC) { + startCoords.r = DEG_TO_RAD(r); + startCoords.d = DEG_TO_RAD(d); + printf("----- startCoords (R, D) is (%f, %f) is deg (%.1f, %.1f) -----\n", startCoords.r, startCoords.d, r, d); + + planeCoords = psProject(&startCoords, &myProj); + endCoords = psDeproject(planeCoords, &myProj); + + if ((fabs((endCoords->r - startCoords.r) / startCoords.r) > ERROR_PERCENT) || + (fabs((endCoords->d - startCoords.d) / startCoords.d) > ERROR_PERCENT)) { + printf("ERROR: \n"); + printf(" startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d); + // printf(" plane Coords is (%f, %f)\n", planeCoords->x, planeCoords->y); + printf(" endCoords (R, D) is (%f, %f)\n", endCoords->r, endCoords->d); + } + + psFree(planeCoords); + psFree(endCoords); + } + } + return(testStatus); +} + // This code will ...