1 // test file for HSO3.hpp and HSO4.hpp
2
3 // (C) Copyright Hubert Holin 2001.
4 // Distributed under the Boost Software License, Version 1.0. (See
5 // accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7
8
9 #include <iostream>
10
11 #include <boost/math/quaternion.hpp>
12
13 #include "HSO3.hpp"
14 #include "HSO4.hpp"
15
16
17 const int number_of_intervals = 5;
18
19 const float pi = ::std::atan(1.0f)*4;
20
21
22
23 void test_SO3();
24
25 void test_SO4();
26
27
main()28 int main()
29
30 {
31 test_SO3();
32
33 test_SO4();
34
35 ::std::cout << "That's all folks!" << ::std::endl;
36 }
37
38
39 //
40 // Test of quaternion and R^3 rotation relationship
41 //
42
test_SO3_spherical()43 void test_SO3_spherical()
44 {
45 ::std::cout << "Testing spherical:" << ::std::endl;
46 ::std::cout << ::std::endl;
47
48 const float rho = 1.0f;
49
50 float theta;
51 float phi1;
52 float phi2;
53
54 for (int idxphi2 = 0; idxphi2 <= number_of_intervals; idxphi2++)
55 {
56 phi2 = (-pi/2)+(idxphi2*pi)/number_of_intervals;
57
58 for (int idxphi1 = 0; idxphi1 <= number_of_intervals; idxphi1++)
59 {
60 phi1 = (-pi/2)+(idxphi1*pi)/number_of_intervals;
61
62 for (int idxtheta = 0; idxtheta <= number_of_intervals; idxtheta++)
63 {
64 theta = -pi+(idxtheta*(2*pi))/number_of_intervals;
65
66 //::std::cout << "theta = " << theta << " ; ";
67 //::std::cout << "phi1 = " << phi1 << " ; ";
68 //::std::cout << "phi2 = " << phi2;
69 //::std::cout << ::std::endl;
70
71 ::boost::math::quaternion<float> q = ::boost::math::spherical(rho, theta, phi1, phi2);
72
73 //::std::cout << "q = " << q << ::std::endl;
74
75 R3_matrix<float> rot = quaternion_to_R3_rotation(q);
76
77 //::std::cout << "rot = ";
78 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
79 //::std::cout << "\t";
80 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
81 //::std::cout << "\t";
82 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
83
84 ::boost::math::quaternion<float> p = R3_rotation_to_quaternion(rot, &q);
85
86 //::std::cout << "p = " << p << ::std::endl;
87
88 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
89
90 //::std::cout << ::std::endl;
91 }
92 }
93 }
94
95 ::std::cout << ::std::endl;
96 }
97
98
test_SO3_semipolar()99 void test_SO3_semipolar()
100 {
101 ::std::cout << "Testing semipolar:" << ::std::endl;
102 ::std::cout << ::std::endl;
103
104 const float rho = 1.0f;
105
106 float alpha;
107 float theta1;
108 float theta2;
109
110 for (int idxalpha = 0; idxalpha <= number_of_intervals; idxalpha++)
111 {
112 alpha = (idxalpha*(pi/2))/number_of_intervals;
113
114 for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
115 {
116 theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
117
118 for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
119 {
120 theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
121
122 //::std::cout << "alpha = " << alpha << " ; ";
123 //::std::cout << "theta1 = " << theta1 << " ; ";
124 //::std::cout << "theta2 = " << theta2;
125 //::std::cout << ::std::endl;
126
127 ::boost::math::quaternion<float> q = ::boost::math::semipolar(rho, alpha, theta1, theta2);
128
129 //::std::cout << "q = " << q << ::std::endl;
130
131 R3_matrix<float> rot = quaternion_to_R3_rotation(q);
132
133 //::std::cout << "rot = ";
134 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
135 //::std::cout << "\t";
136 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
137 //::std::cout << "\t";
138 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
139
140 ::boost::math::quaternion<float> p = R3_rotation_to_quaternion(rot, &q);
141
142 //::std::cout << "p = " << p << ::std::endl;
143
144 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
145
146 //::std::cout << ::std::endl;
147 }
148 }
149 }
150
151 ::std::cout << ::std::endl;
152 }
153
154
test_SO3_multipolar()155 void test_SO3_multipolar()
156 {
157 ::std::cout << "Testing multipolar:" << ::std::endl;
158 ::std::cout << ::std::endl;
159
160 float rho1;
161 float rho2;
162
163 float theta1;
164 float theta2;
165
166 for (int idxrho = 0; idxrho <= number_of_intervals; idxrho++)
167 {
168 rho1 = (idxrho*1.0f)/number_of_intervals;
169 rho2 = ::std::sqrt(1.0f-rho1*rho1);
170
171 for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
172 {
173 theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
174
175 for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
176 {
177 theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
178
179 //::std::cout << "rho1 = " << rho1 << " ; ";
180 //::std::cout << "theta1 = " << theta1 << " ; ";
181 //::std::cout << "theta2 = " << theta2;
182 //::std::cout << ::std::endl;
183
184 ::boost::math::quaternion<float> q = ::boost::math::multipolar(rho1, theta1, rho2, theta2);
185
186 //::std::cout << "q = " << q << ::std::endl;
187
188 R3_matrix<float> rot = quaternion_to_R3_rotation(q);
189
190 //::std::cout << "rot = ";
191 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
192 //::std::cout << "\t";
193 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
194 //::std::cout << "\t";
195 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
196
197 ::boost::math::quaternion<float> p = R3_rotation_to_quaternion(rot, &q);
198
199 //::std::cout << "p = " << p << ::std::endl;
200
201 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
202
203 //::std::cout << ::std::endl;
204 }
205 }
206 }
207
208 ::std::cout << ::std::endl;
209 }
210
211
test_SO3_cylindrospherical()212 void test_SO3_cylindrospherical()
213 {
214 ::std::cout << "Testing cylindrospherical:" << ::std::endl;
215 ::std::cout << ::std::endl;
216
217 float t;
218
219 float radius;
220 float longitude;
221 float latitude;
222
223 for (int idxt = 0; idxt <= number_of_intervals; idxt++)
224 {
225 t = -1.0f+(idxt*2.0f)/number_of_intervals;
226 radius = ::std::sqrt(1.0f-t*t);
227
228 for (int idxlatitude = 0; idxlatitude <= number_of_intervals; idxlatitude++)
229 {
230 latitude = (-pi/2)+(idxlatitude*pi)/number_of_intervals;
231
232 for (int idxlongitude = 0; idxlongitude <= number_of_intervals; idxlongitude++)
233 {
234 longitude = -pi+(idxlongitude*(2*pi))/number_of_intervals;
235
236 //::std::cout << "t = " << t << " ; ";
237 //::std::cout << "longitude = " << longitude;
238 //::std::cout << "latitude = " << latitude;
239 //::std::cout << ::std::endl;
240
241 ::boost::math::quaternion<float> q = ::boost::math::cylindrospherical(t, radius, longitude, latitude);
242
243 //::std::cout << "q = " << q << ::std::endl;
244
245 R3_matrix<float> rot = quaternion_to_R3_rotation(q);
246
247 //::std::cout << "rot = ";
248 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
249 //::std::cout << "\t";
250 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
251 //::std::cout << "\t";
252 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
253
254 ::boost::math::quaternion<float> p = R3_rotation_to_quaternion(rot, &q);
255
256 //::std::cout << "p = " << p << ::std::endl;
257
258 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
259
260 //::std::cout << ::std::endl;
261 }
262 }
263 }
264
265 ::std::cout << ::std::endl;
266 }
267
268
test_SO3_cylindrical()269 void test_SO3_cylindrical()
270 {
271 ::std::cout << "Testing cylindrical:" << ::std::endl;
272 ::std::cout << ::std::endl;
273
274 float r;
275 float angle;
276
277 float h1;
278 float h2;
279
280 for (int idxh2 = 0; idxh2 <= number_of_intervals; idxh2++)
281 {
282 h2 = -1.0f+(idxh2*2.0f)/number_of_intervals;
283
284 for (int idxh1 = 0; idxh1 <= number_of_intervals; idxh1++)
285 {
286 h1 = ::std::sqrt(1.0f-h2*h2)*(-1.0f+(idxh2*2.0f)/number_of_intervals);
287 r = ::std::sqrt(1.0f-h1*h1-h2*h2);
288
289 for (int idxangle = 0; idxangle <= number_of_intervals; idxangle++)
290 {
291 angle = -pi+(idxangle*(2*pi))/number_of_intervals;
292
293 //::std::cout << "angle = " << angle << " ; ";
294 //::std::cout << "h1 = " << h1;
295 //::std::cout << "h2 = " << h2;
296 //::std::cout << ::std::endl;
297
298 ::boost::math::quaternion<float> q = ::boost::math::cylindrical(r, angle, h1, h2);
299
300 //::std::cout << "q = " << q << ::std::endl;
301
302 R3_matrix<float> rot = quaternion_to_R3_rotation(q);
303
304 //::std::cout << "rot = ";
305 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
306 //::std::cout << "\t";
307 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
308 //::std::cout << "\t";
309 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
310
311 ::boost::math::quaternion<float> p = R3_rotation_to_quaternion(rot, &q);
312
313 //::std::cout << "p = " << p << ::std::endl;
314
315 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
316
317 //::std::cout << ::std::endl;
318 }
319 }
320 }
321
322 ::std::cout << ::std::endl;
323 }
324
325
test_SO3()326 void test_SO3()
327 {
328 ::std::cout << "Testing SO3:" << ::std::endl;
329 ::std::cout << ::std::endl;
330
331 test_SO3_spherical();
332
333 test_SO3_semipolar();
334
335 test_SO3_multipolar();
336
337 test_SO3_cylindrospherical();
338
339 test_SO3_cylindrical();
340 }
341
342
343 //
344 // Test of quaternion and R^4 rotation relationship
345 //
346
test_SO4_spherical()347 void test_SO4_spherical()
348 {
349 ::std::cout << "Testing spherical:" << ::std::endl;
350 ::std::cout << ::std::endl;
351
352 const float rho1 = 1.0f;
353 const float rho2 = 1.0f;
354
355 float theta1;
356 float phi11;
357 float phi21;
358
359 float theta2;
360 float phi12;
361 float phi22;
362
363 for (int idxphi21 = 0; idxphi21 <= number_of_intervals; idxphi21++)
364 {
365 phi21 = (-pi/2)+(idxphi21*pi)/number_of_intervals;
366
367 for (int idxphi22 = 0; idxphi22 <= number_of_intervals; idxphi22++)
368 {
369 phi22 = (-pi/2)+(idxphi22*pi)/number_of_intervals;
370
371 for (int idxphi11 = 0; idxphi11 <= number_of_intervals; idxphi11++)
372 {
373 phi11 = (-pi/2)+(idxphi11*pi)/number_of_intervals;
374
375 for (int idxphi12 = 0; idxphi12 <= number_of_intervals; idxphi12++)
376 {
377 phi12 = (-pi/2)+(idxphi12*pi)/number_of_intervals;
378
379 for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
380 {
381 theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
382
383 for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
384 {
385 theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
386
387 //::std::cout << "theta1 = " << theta1 << " ; ";
388 //::std::cout << "phi11 = " << phi11 << " ; ";
389 //::std::cout << "phi21 = " << phi21;
390 //::std::cout << "theta2 = " << theta2 << " ; ";
391 //::std::cout << "phi12 = " << phi12 << " ; ";
392 //::std::cout << "phi22 = " << phi22;
393 //::std::cout << ::std::endl;
394
395 ::boost::math::quaternion<float> p1 = ::boost::math::spherical(rho1, theta1, phi11, phi21);
396
397 //::std::cout << "p1 = " << p1 << ::std::endl;
398
399 ::boost::math::quaternion<float> q1 = ::boost::math::spherical(rho2, theta2, phi12, phi22);
400
401 //::std::cout << "q1 = " << q1 << ::std::endl;
402
403 ::std::pair< ::boost::math::quaternion<float> , ::boost::math::quaternion<float> > pq1 =
404 ::std::make_pair(p1,q1);
405
406 R4_matrix<float> rot = quaternions_to_R4_rotation(pq1);
407
408 //::std::cout << "rot = ";
409 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << "\t" << rot.a14 << ::std::endl;
410 //::std::cout << "\t";
411 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << "\t" << rot.a24 << ::std::endl;
412 //::std::cout << "\t";
413 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << "\t" << rot.a34 << ::std::endl;
414 //::std::cout << "\t";
415 //::std::cout << "\t" << rot.a41 << "\t" << rot.a42 << "\t" << rot.a43 << "\t" << rot.a44 << ::std::endl;
416
417 ::std::pair< ::boost::math::quaternion<float> , ::boost::math::quaternion<float> > pq2 =
418 R4_rotation_to_quaternions(rot, &pq1);
419
420 //::std::cout << "p1 = " << pq.first << ::std::endl;
421 //::std::cout << "p2 = " << pq.second << ::std::endl;
422
423 //::std::cout << "round trip discrepancy: " << ::std::sqrt(::boost::math::norm(pq1.first-pq2.first)+::boost::math::norm(pq1.second-pq2.second)) << ::std::endl;
424
425 //::std::cout << ::std::endl;
426 }
427 }
428 }
429 }
430 }
431 }
432
433 ::std::cout << ::std::endl;
434 }
435
436
test_SO4()437 void test_SO4()
438 {
439 ::std::cout << "Testing SO4:" << ::std::endl;
440 ::std::cout << ::std::endl;
441
442 test_SO4_spherical();
443 }
444
445
446