Changeset 5319
- Timestamp:
- Oct 13, 2005, 2:07:37 PM (21 years ago)
- Location:
- trunk/psLib
- Files:
-
- 6 edited
-
src/astro/psSphereOps.c (modified) (5 diffs)
-
src/astro/psSphereOps.h (modified) (2 diffs)
-
test/astro/Makefile.am (modified) (2 diffs)
-
test/astro/tst_psSphereOps.c (modified) (6 diffs)
-
test/astro/verified/tst_psSphereOps.stderr (modified) (3 diffs)
-
test/astro/verified/tst_psSphereOps.stdout (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psSphereOps.c
r5306 r5319 6 6 * 7 7 * @author Robert DeSonia, MHPCC 8 * @author Dave Robbins, MHPCC 8 9 * 9 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-10-1 3 20:23:57 $10 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2005-10-14 00:07:37 $ 11 12 * 12 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 62 63 // calculate t*s*r. 63 64 psSphereRot* result = psSphereRotCombine(NULL,&t,&s); 64 psMemSetDeallocator(result, (psFreeFunc)sphereRotFree);65 65 psSphereRotCombine(result,result,&r); 66 66 … … 79 79 double q3) 80 80 { 81 psSphereRot* rot = psAlloc(sizeof(psSphereRot)); 81 psSphereRot* rot = (psSphereRot*)psAlloc(sizeof(psSphereRot)); 82 psMemSetDeallocator(rot, (psFreeFunc)sphereRotFree); 82 83 83 84 double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); … … 113 114 sin(coord->d), 114 115 0.0); 115 psSphereRot* coordQuatConjugate = psSphereRotQuat(116 coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3);117 coordQuat = psSphereRotInvert(coordQuat);116 // psSphereRot* coordQuatConjugate = psSphereRotQuat( 117 // coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3); 118 // coordQuat = psSphereRotInvert(coordQuat); 118 119 119 120 // calculate q=(rp)r' 120 coordQuat = psSphereRotCombine(coordQuat, transform, coordQuat);121 coordQuat = psSphereRotCombine(coordQuat, coordQuat, coordQuatConjugate);121 // coordQuat = psSphereRotCombine(coordQuat, transform, coordQuat); 122 // coordQuat = psSphereRotCombine(coordQuat, coordQuat, coordQuatConjugate); 122 123 // N.B., we can recycle coordQuat right away due to the implementation of 123 124 // psSphereRotCombine; it puts the input values in a local variable first 124 125 125 out->r = atan2(coordQuat->q1,coordQuat->q0); 126 out->d = asin(coordQuat->q2); 126 // out->r = atan2(coordQuat->q1,coordQuat->q0); 127 // out->d = asin(coordQuat->q2); 128 129 130 // psSphereRot *inv = psSphereRotInvert(transform); 131 psSphereRot *inv = (psSphereRot*)psAlloc(sizeof(psSphereRot)); 132 *inv = *transform; 133 psSphereRot *result = psSphereRotCombine(NULL, transform, coordQuat); 134 psSphereRotInvert(inv); 135 psSphereRotCombine(result, result, inv); 136 out->r = atan2(result->q1, result->q0); 137 out->d = asin(result->q2); 138 out->rErr = 0.0; 139 out->dErr = 0.0; 140 psFree(inv); 141 psFree(result); 127 142 128 143 psFree(coordQuat); 129 psFree(coordQuatConjugate);144 // psFree(coordQuatConjugate); 130 145 131 146 return out; … … 141 156 if (out == NULL) { 142 157 out = (psSphereRot* ) psAlloc(sizeof(psSphereRot)); 158 psMemSetDeallocator(out, (psFreeFunc)sphereRotFree); 143 159 } 144 160 -
trunk/psLib/src/astro/psSphereOps.h
r5306 r5319 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1. 4$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-10-1 3 20:23:57 $9 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-10-14 00:07:37 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 73 73 ); 74 74 75 psSphereRot* psSphereRotAlloc2(76 double alphaP, ///< north pole latitude77 double deltaP, ///< north pole longitude78 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).79 );80 81 psSphereRot* psSphereRotAlloc3(82 double alphaP, ///< north pole latitude83 double deltaP, ///< north pole longitude84 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).85 );86 87 psSphereRot* psSphereRotAlloc4(88 double alphaP, ///< north pole latitude89 double deltaP, ///< north pole longitude90 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).91 );92 93 psSphereRot* psSphereRotAlloc5(94 double alphaP, ///< north pole latitude95 double deltaP, ///< north pole longitude96 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).97 );98 99 psSphereRot* psSphereRotAlloc6(100 double alphaP, ///< north pole latitude101 double deltaP, ///< north pole longitude102 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).103 );104 105 psSphereRot* psSphereRotAlloc7(106 double alphaP, ///< north pole latitude107 double deltaP, ///< north pole longitude108 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).109 );110 111 psSphereRot* psSphereRotAlloc8(112 double alphaP, ///< north pole latitude113 double deltaP, ///< north pole longitude114 double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares).115 );116 117 118 75 /** Checks the type of a particular pointer. 119 76 * -
trunk/psLib/test/astro/Makefile.am
r5299 r5319 12 12 tst_psTime_04 \ 13 13 tst_psCoord \ 14 tst_psEarthOrientation 15 #tst_psSphereOps14 tst_psEarthOrientation \ 15 tst_psSphereOps 16 16 17 17 tst_psTime_01_SOURCES = tst_psTime_01.c … … 20 20 tst_psTime_04_SOURCES = tst_psTime_04.c 21 21 tst_psCoord_SOURCES = tst_psCoord.c 22 #tst_psSphereOps_SOURCES = tst_psSphereOps.c22 tst_psSphereOps_SOURCES = tst_psSphereOps.c 23 23 tst_psCoord01_SOURCES = tst_psCoord01.c 24 24 tst_psEarthOrientation_SOURCES = tst_psEarthOrientation.c -
trunk/psLib/test/astro/tst_psSphereOps.c
r5306 r5319 1 /** @file tst_ps Coord.c1 /** @file tst_psSphereOps.c 2 2 * 3 * @brief The code will ... 3 * @brief The code will ..... Work .... 4 4 * 5 5 * 6 * @author GLG, MHPCC6 * @author d-Rob, MHPCC 7 7 * 8 * @version $Revision: 1. 4$ $Name: not supported by cvs2svn $9 * @date $Date: 2005-10-1 3 20:23:57 $8 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2005-10-14 00:07:37 $ 10 10 * 11 11 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 15 15 #include "pslib_strict.h" 16 16 static psS32 testSphereRotAlloc(void); 17 static psS32 testSphereRotQuat(void); 17 18 static psS32 testSphereRotApply1(void); 18 static psS32 testSphereRotApply2(void);19 19 static psS32 testSphereRotApplyCelestial(void); 20 20 static psS32 testSphereRotPrecess(void); … … 22 22 testDescription tests[] = { 23 23 {testSphereRotAlloc, 819, "psSphereRotAlloc()", 0, false}, 24 {testSphereRot Apply1, 820, "psSphereRotApply()", 0, false},25 {testSphereRotApply 2, 820, "psSphereRotApply()", 0, false},26 {testSphereRotApplyCelestial, 82 0, "psSphereRotApply()", 0, false},27 // {testSphereRotApplyCelestial, 821, "psSphereRotICRSToEcliptic()", true},28 // {testSphereRotApplyCelestial, 822, "psSphereRotEclipticToICRS()", 0, true},29 // {testSphereRotApplyCelestial, 824, "psSphereRotICRSToGalactic()", 0, true},30 // {testSphereRotApplyCelestial, 823, "psSphereRotGalacticToICRS()", 0, true},31 {testSphereRotPrecess, 825, "psSphereRotPrecess()", 0, false},24 {testSphereRotQuat, 820, "psSphereRotQuat()", 0, false}, 25 {testSphereRotApply1, 821, "psSphereRotApply()", 0, false}, 26 {testSphereRotApplyCelestial, 822, "psSphereRotApplyCel()", 0, false}, 27 {testSphereRotPrecess, 823, "psSphereRotPrecess()", 0, false}, 28 // {testSphereRotApplyCelestial, 821, "psSphereRotICRSToEcliptic()", true}, 29 // {testSphereRotApplyCelestial, 822, "psSphereRotEclipticToICRS()", 0, true}, 30 // {testSphereRotApplyCelestial, 824, "psSphereRotICRSToGalactic()", 0, true}, 31 // {testSphereRotApplyCelestial, 823, "psSphereRotGalacticToICRS()", 0, true}, 32 32 {NULL} 33 33 }; … … 50 50 #define ALPHA_P M_PI/6 51 51 #define DELTA_P M_PI/4 52 #define PHI_P M_PI/ 252 #define PHI_P M_PI/3 53 53 54 54 psS32 testSphereRotAlloc( void ) … … 73 73 double q2 = sin(a2)*cos(DELTA_P/2); 74 74 75 printf("\n q0=%lf, q1=%lf, q2=%lf, q3=%lf,\n", q0, q1, q2, q3);76 printf("myq0=%lf, myq1=%lf, myq2=%lf, myq3=%lf\n", myST->q0, myST->q1, myST->q2, myST->q3);75 // printf("\n q0=%lf, q1=%lf, q2=%lf, q3=%lf,\n", q0, q1, q2, q3); 76 // printf("myq0=%lf, myq1=%lf, myq2=%lf, myq3=%lf\n", myST->q0, myST->q1, myST->q2, myST->q3); 77 77 if (FLT_EPSILON < fabs(q0 - myST->q0)) { 78 78 psError(PS_ERR_UNKNOWN,true,"myST->q0 is %lf, should be %lf\n", myST->q0, q0); … … 98 98 } 99 99 100 psS32 testSphereRotQuat(void) 101 { 102 double a0 = (ALPHA_P - PHI_P)/2.0; 103 double a1 = (ALPHA_P - PHI_P)/2.0; 104 double a2 = (ALPHA_P + PHI_P)/2.0; 105 double a3 = (ALPHA_P + PHI_P)/2.0; 106 107 double q3 = cos(a3)*cos(DELTA_P/2); 108 double q0 = sin(a0)*sin(DELTA_P/2); 109 double q1 = cos(a1)*sin(DELTA_P/2); 110 double q2 = sin(a2)*cos(DELTA_P/2); 111 112 psSphereRot *myST = psSphereRotQuat(q0*2.0, q1*2.0, q2*2.0, q3*2.0); 113 // Verify null not returned 114 if(myST == NULL) { 115 psError(PS_ERR_UNKNOWN,true,"Returned NULL with valid parameters"); 116 return 1; 117 } 118 119 // printf("\n q0=%lf, q1=%lf, q2=%lf, q3=%lf,\n", q0, q1, q2, q3); 120 // printf("myq0=%lf, myq1=%lf, myq2=%lf, myq3=%lf\n", myST->q0, myST->q1, myST->q2, myST->q3); 121 if (FLT_EPSILON < fabs(q0 - myST->q0)) { 122 psError(PS_ERR_UNKNOWN,true,"myST->q0 is %lf, should be %lf\n", myST->q0, q0); 123 return 2; 124 } 125 if (FLT_EPSILON < fabs(q1 - myST->q1)) { 126 psError(PS_ERR_UNKNOWN,true,"myST->q1 is %f, should be %f\n", myST->q1, q1); 127 return 3; 128 } 129 if (FLT_EPSILON < fabs(q2 - myST->q2)) { 130 psError(PS_ERR_UNKNOWN,true,"myST->q2 is %f, should be %f\n", myST->q2, q2); 131 return 4; 132 } 133 if (FLT_EPSILON < fabs(q3 - myST->q3)) { 134 psError(PS_ERR_UNKNOWN,true,"myST->q0 is %f, should be %f\n", myST->q3, q3); 135 return 5; 136 } 137 138 // Free data structure 139 psFree(myST); 140 141 return 0; 142 } 143 144 100 145 // We do a simple identity transformation on a few RA, DEC pairs. 101 146 psS32 testSphereRotApply1( void ) 102 147 { 103 /* 104 psSphere *in = psSphereAlloc(); 105 psSphere *out = psSphereAlloc(); 106 psSphere *temp = NULL; 107 psSphere *rc = NULL; 108 psSphereRot *myST = psSphereRotAlloc(0.0, 0.0, 0.0); 109 110 for (float r=0.0;r<180.0;r+=DEG_INC) { 111 for (float d=0.0;d<90.0;d+=DEG_INC) { 112 in->r = DEG_TO_RAD(r); 113 in->d = DEG_TO_RAD(d); 114 in->rErr = 0.0; 115 in->dErr = 0.0; 116 117 if(psSphereRotApply(out, myST, in) != out) { 118 psError(PS_ERR_UNKNOWN,true,"Did not return output pointer."); 119 return 1; 120 } 121 if (ERROR_TOL < fabs(out->r - in->r)) { 122 psError(PS_ERR_UNKNOWN,true,"out->r is %f, should be %f\n", out->r, in->r); 123 return 2; 124 } 125 if (ERROR_TOL < fabs(out->d - in->d)) { 126 psError(PS_ERR_UNKNOWN,true,"out->d is %f, should be %f\n", out->d, in->d); 127 return 3; 128 } 129 } 130 } 131 132 // Verify new sphere object is created if out parameter NULL 133 temp = psSphereRotApply(NULL, myST, in); 134 if ( temp == NULL) { 135 psError(PS_ERR_UNKNOWN,true,"Returned NULL when out parameter was null"); 136 return 4; 137 } 138 psFree(temp); 139 140 // Verify NULL returned if transform structure null 141 psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error"); 142 rc = psSphereRotApply(NULL, NULL, in); 143 if (rc != NULL) { 144 psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() did not return NULL."); 145 return 5; 146 } 147 148 // Verify NULL returned when input sphere is NULL 149 psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error"); 150 rc = psSphereRotApply(NULL, myST, NULL); 151 if (rc != NULL) { 152 psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() did not return NULL"); 153 return 6; 154 } 155 156 psFree(myST); 157 psFree(out); 158 psFree(in); 159 */ 148 149 psSphere *in = psSphereAlloc(); 150 psSphere *out = psSphereAlloc(); 151 psSphere *temp = NULL; 152 psSphere *rc = NULL; 153 // psSphereRot *myST = psSphereRotAlloc(0.0, 0.0, 0.0); 154 psSphereRot *myST = psSphereRotAlloc(ALPHA_P, DELTA_P, PHI_P); 155 156 for (float r=0.0;r<180.0;r+=DEG_INC) { 157 for (float d=0.0;d<90.0;d+=DEG_INC) { 158 in->r = DEG_TO_RAD(r); 159 in->d = DEG_TO_RAD(d); 160 in->rErr = 0.0; 161 in->dErr = 0.0; 162 163 if(psSphereRotApply(out, myST, in) != out) { 164 psError(PS_ERR_UNKNOWN,true,"Did not return output pointer."); 165 return 1; 166 } 167 psSphereRotInvert(myST); 168 psSphereRotApply(out, myST, out); 169 170 if (ERROR_TOL < fabs(out->r - in->r)) { 171 psError(PS_ERR_UNKNOWN,true,"out->r is %f, should be %f\n", out->r, in->r); 172 return 2; 173 } 174 if (ERROR_TOL < fabs(out->d - in->d)) { 175 psError(PS_ERR_UNKNOWN,true,"out->d is %f, should be %f\n", out->d, in->d); 176 return 3; 177 } 178 } 179 } 180 // Verify new sphere object is created if out parameter NULL 181 temp = psSphereRotApply(NULL, myST, in); 182 if ( temp == NULL) { 183 psError(PS_ERR_UNKNOWN,true,"Returned NULL when out parameter was null"); 184 return 4; 185 } 186 psFree(temp); 187 188 // Verify NULL returned if transform structure null 189 psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error"); 190 rc = psSphereRotApply(NULL, NULL, in); 191 if (rc != NULL) { 192 psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() did not return NULL."); 193 return 5; 194 } 195 196 // Verify NULL returned when input sphere is NULL 197 psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error"); 198 rc = psSphereRotApply(NULL, myST, NULL); 199 if (rc != NULL) { 200 psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() did not return NULL"); 201 return 6; 202 } 203 204 psFree(myST); 205 psFree(out); 206 psFree(in); 207 160 208 return 0; 161 209 } 162 210 163 /******************************************************************************164 testSphereRotApply2(): This test verifies that psSphereRotApply()165 works properly. We create two psSphereRots: a forward transform and a166 reverse transform (which is the mathematical inverse of the forward transform).167 We apply both transforms to several spherical coordinates and ensure that the168 original input coordinate is obtained after applying both transforms.169 170 XXX: We currently test the alpha and delta offsets independently. Attempts to171 test them both concurrently failed. Determine why this is. Are the following172 spherical transforms not mathematical inverses?173 psSphereRotAlloc(X, Y, 0.0)174 psSphereRotAlloc(-X, -Y, 0.0)175 *****************************************************************************/176 211 #define ERROR_PERCENT 0.01 177 psS32 testSphereRotApply2( void )178 {179 /*180 psS32 testStatus = 0;181 psSphere in;182 psSphere out;183 psSphere out2;184 psSphereRot *mySphereRotForward = NULL;185 psSphereRot *mySphereRotReverse = NULL;186 187 188 mySphereRotForward = psSphereRotAlloc(DEG_TO_RAD(22.0),189 0.0,190 0.0);191 mySphereRotReverse = psSphereRotAlloc(DEG_TO_RAD(-22.0),192 0.0,193 0.0);194 195 for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {196 for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {197 in.r = DEG_TO_RAD(r);198 in.d = DEG_TO_RAD(d);199 in.rErr = 0.0;200 in.dErr = 0.0;201 202 psSphereRotApply(&out, mySphereRotForward, &in);203 psSphereRotApply(&out2, mySphereRotReverse, &out);204 205 if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||206 (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {207 printf("ERROR: \n");208 printf("Input coords (R, D) are (%f, %f)\n", in.r, in.d);209 printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);210 testStatus = 4;211 }212 }213 }214 psFree(mySphereRotForward);215 psFree(mySphereRotReverse);216 217 mySphereRotForward = psSphereRotAlloc(0.0,218 DEG_TO_RAD(33.0),219 0.0);220 mySphereRotReverse = psSphereRotAlloc(0.0,221 DEG_TO_RAD(-33.0),222 0.0);223 for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {224 for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {225 in.r = DEG_TO_RAD(r);226 in.d = DEG_TO_RAD(d);227 in.rErr = 0.0;228 in.dErr = 0.0;229 230 psSphereRotApply(&out, mySphereRotForward, &in);231 psSphereRotApply(&out2, mySphereRotReverse, &out);232 233 if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||234 (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {235 printf("ERROR: \n");236 printf("Input coords (R, D) are (%f, %f)\n", in.r, in.d);237 printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);238 testStatus = 4;239 }240 }241 }242 psFree(mySphereRotForward);243 psFree(mySphereRotReverse);244 245 return(testStatus);246 */247 return 0;248 }249 212 250 213 psS32 testSphereRotApplyCelestial( void) -
trunk/psLib/test/astro/verified/tst_psSphereOps.stderr
r5235 r5319 10 10 /***************************** TESTPOINT ******************************************\ 11 11 * TestFile: tst_psSphereOps.c * 12 * TestPoint: psCoord{psSphereRot Apply()}*12 * TestPoint: psCoord{psSphereRotQuat()} * 13 13 * TestType: Positive * 14 14 \**********************************************************************************/ 15 15 16 16 17 ---> TESTPOINT PASSED (psCoord{psSphereRot Apply()} | tst_psSphereOps.c)17 ---> TESTPOINT PASSED (psCoord{psSphereRotQuat()} | tst_psSphereOps.c) 18 18 19 19 /***************************** TESTPOINT ******************************************\ … … 23 23 \**********************************************************************************/ 24 24 25 <DATE><TIME>|<HOST>|I|testSphereRotApply1 26 Following should generate an error 27 <DATE><TIME>|<HOST>|E|psSphereRotApply (FILE:LINENO) 28 Unallowable operation: transform is NULL. 29 <DATE><TIME>|<HOST>|I|testSphereRotApply1 30 Following should generate an error 31 <DATE><TIME>|<HOST>|E|psSphereRotApply (FILE:LINENO) 32 Unallowable operation: coord is NULL. 25 33 26 34 ---> TESTPOINT PASSED (psCoord{psSphereRotApply()} | tst_psSphereOps.c) … … 28 36 /***************************** TESTPOINT ******************************************\ 29 37 * TestFile: tst_psSphereOps.c * 30 * TestPoint: psCoord{psSphereRotApply ()}*38 * TestPoint: psCoord{psSphereRotApplyCel()} * 31 39 * TestType: Positive * 32 40 \**********************************************************************************/ 33 41 34 42 35 ---> TESTPOINT PASSED (psCoord{psSphereRotApply ()} | tst_psSphereOps.c)43 ---> TESTPOINT PASSED (psCoord{psSphereRotApplyCel()} | tst_psSphereOps.c) 36 44 37 45 /***************************** TESTPOINT ******************************************\ -
trunk/psLib/test/astro/verified/tst_psSphereOps.stdout
r5306 r5319 1 2 q0=-0.191342, q1=0.331414, q2=0.800103, q3=0.461940,3 myq0=-0.191342, myq1=0.331414, myq2=0.800103, myq3=0.461940
Note:
See TracChangeset
for help on using the changeset viewer.
