1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 from math import pi
24 from numpy import array, float64, zeros
25 from unittest import TestCase
26
27
28 import dep_check
29 from lib.frame_order.format import print_frame_order_2nd_degree
30 from lib.frame_order.free_rotor import compile_2nd_matrix_free_rotor
31 from lib.frame_order.iso_cone import compile_2nd_matrix_iso_cone
32 from lib.frame_order.iso_cone_free_rotor import compile_2nd_matrix_iso_cone_free_rotor
33 from lib.frame_order.iso_cone_torsionless import compile_2nd_matrix_iso_cone_torsionless
34 from lib.frame_order.pseudo_ellipse import compile_2nd_matrix_pseudo_ellipse
35 from lib.frame_order.pseudo_ellipse_free_rotor import compile_2nd_matrix_pseudo_ellipse_free_rotor
36 from lib.frame_order.pseudo_ellipse_torsionless import compile_2nd_matrix_pseudo_ellipse_torsionless
37 from lib.frame_order.rotor import compile_2nd_matrix_rotor
38 from lib.frame_order.matrix_ops import reduce_alignment_tensor
39 from lib.geometry.coord_transform import cartesian_to_spherical, spherical_to_cartesian
40 from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R
41 from lib.linear_algebra.kronecker_product import kron_prod, transpose_23
42 from status import Status; status = Status()
43
44
46 """Unit tests for the lib.frame_order_matrix_ops relax module."""
47
48 - def __init__(self, methodName='runTest'):
49 """Skip the tests if scipy is not installed.
50
51 @keyword methodName: The name of the test.
52 @type methodName: str
53 """
54
55
56 if not dep_check.scipy_module:
57
58 status.skipped_tests.append([methodName, 'Scipy', 'unit'])
59
60
61 super(Test_matrix_ops, self).__init__(methodName)
62
63
65 """Initialise a few data structures for the tests."""
66
67
68 self.f2_temp = zeros((9, 9), float64)
69 self.R_temp = zeros((3, 3), float64)
70 self.z_axis = array([0, 0, 1], float64)
71 self.cone_axis = zeros(3, float64)
72
73
74 self.setup_identity()
75
76
77 self.setup_identity_free_rotor()
78
79
80 self.out_of_frame_alpha = 1.10714871779
81 self.out_of_frame_beta = 0.841068670568
82 self.out_of_frame_gamma = 5.81953769818
83
84
86 """Calculate the Kronecker product of the eigenframe rotation for the z-axis based frame."""
87
88
89 spherical_to_cartesian([1.0, axis_theta, axis_phi], self.cone_axis)
90
91
92 two_vect_to_R(self.z_axis, self.cone_axis, self.R_temp)
93
94
95 return kron_prod(self.R_temp, self.R_temp)
96
97
99 """Calculate the Kronecker product of the eigenframe rotation for the full eigenframe."""
100
101
102 euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, self.R_temp)
103
104
105 return kron_prod(self.R_temp, self.R_temp)
106
107
109 """Set up a few identity matrices."""
110
111
112 self.I_order = zeros((9, 9), float64)
113 for i in range(9):
114 self.I_order[i, i] = 1.0
115
116
117 data = [[0, 0, 1.0/3.0],
118 [4, 4, 1.0/3.0],
119 [8, 8, 1.0/3.0],
120 [0, 4, 1.0/3.0],
121 [4, 0, 1.0/3.0],
122 [0, 8, 1.0/3.0],
123 [8, 0, 1.0/3.0],
124 [4, 8, 1.0/3.0],
125 [8, 4, 1.0/3.0]]
126 self.I_disorder = zeros((9, 9), float64)
127 for i, j, val in data:
128 self.I_disorder[i, j] = val
129
130
131 data = [[0, 0, 1.0/3.0],
132 [4, 4, 1.0/3.0],
133 [8, 8, 1.0/3.0],
134 [0, 4, 1.0/3.0],
135 [4, 0, 1.0/3.0],
136 [0, 8, 1.0/3.0],
137 [8, 0, 1.0/3.0],
138 [4, 8, 1.0/3.0],
139 [8, 4, 1.0/3.0],
140 [1, 3, -0.25],
141 [3, 1, -0.25],
142 [1, 1, 0.25],
143 [3, 3, 0.25]]
144 self.f2_half_cone = zeros((9, 9), float64)
145 for i, j, val in data:
146 self.f2_half_cone[i, j] = val
147
148
149 data = [[0, 0, 1.0/3.0],
150 [4, 4, 1.0/3.0],
151 [8, 8, 1.0/3.0],
152 [0, 4, 1.0/3.0],
153 [4, 0, 1.0/3.0],
154 [0, 8, 1.0/3.0],
155 [8, 0, 1.0/3.0],
156 [4, 8, 1.0/3.0],
157 [8, 4, 1.0/3.0],
158 [5, 7, -0.25],
159 [7, 5, -0.25],
160 [7, 7, 0.25],
161 [5, 5, 0.25]]
162 self.f2_half_cone_90_y = zeros((9, 9), float64)
163 for i, j, val in data:
164 self.f2_half_cone_90_y[i, j] = val
165
166
168 """Set up a few identity matrices."""
169
170
171 data = [[0, 0, 0.5],
172 [1, 1, 0.5],
173 [3, 3, 0.5],
174 [4, 4, 0.5],
175 [0, 4, 0.5],
176 [4, 0, 0.5],
177 [1, 3, -0.5],
178 [3, 1, -0.5],
179 [8, 8, 1.0]]
180 self.I_order_free_rotor = zeros((9, 9), float64)
181 for i, j, val in data:
182 self.I_order_free_rotor[i, j] = val
183
184
185 data = [[0, 0, 0.25],
186 [1, 1, 0.125],
187 [3, 3, 0.125],
188 [4, 4, 0.25],
189 [0, 4, 0.25],
190 [4, 0, 0.25],
191 [1, 3, -0.125],
192 [3, 1, -0.125],
193 [0, 8, 0.5],
194 [8, 0, 0.5],
195 [4, 8, 0.5],
196 [8, 4, 0.5]]
197 self.I_disorder_free_rotor = zeros((9, 9), float64)
198 for i, j, val in data:
199 self.I_disorder_free_rotor[i, j] = val
200
201
203 """Check the operation of the compile_2nd_matrix_free_rotor() function."""
204
205
206 real = array(
207 [[ 0.5001, 0.0001, 0, 0.0001, 0.4999, 0, 0, 0, 0],
208 [ -0.0001, 0.5001, 0, -0.4999, 0.0001, 0, 0, 0, 0],
209 [ 0, 0, 0.0006, 0, 0, -0.0005, 0, 0, 0],
210 [ -0.0001, -0.4999, 0, 0.5001, 0.0001, 0, 0, 0, 0],
211 [ 0.4999, -0.0001, 0, -0.0001, 0.5001, 0, 0, 0, 0],
212 [ 0, 0, 0.0005, 0, 0, 0.0006, 0, 0, 0],
213 [ 0, 0, 0, 0, 0, 0, 0.0006, -0.0005, 0],
214 [ 0, 0, 0, 0, 0, 0, 0.0005, 0.0006, 0],
215 [ 0, 0, 0, 0, 0, 0, 0, 0, 1.0000]])
216
217
218 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
219
220
221 f2 = compile_2nd_matrix_free_rotor(self.f2_temp, Rx2_eigen)
222
223
224 print_frame_order_2nd_degree(real, "real")
225 print_frame_order_2nd_degree(f2, "calculated")
226 print_frame_order_2nd_degree(real-f2, "difference")
227
228
229 for i in range(9):
230 for j in range(9):
231 print("Element %s, %s; diff %s." % (i, j, f2[i, j] - real[i, j]))
232 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
233
234
236 """Check the operation of the compile_2nd_matrix_free_rotor() function."""
237
238
239 real = array(
240 [[ 0.3367, -0.0100, -0.0307, -0.0100, 0.3521, -0.0152, -0.0307, -0.0152, 0.3112],
241 [ -0.0104, 0.3520, -0.0152, -0.2908, -0.0559, 0.2602, 0.1989, -0.1685, 0.0664],
242 [ -0.0306, -0.0155, 0.3112, 0.1991, -0.1683, 0.0666, 0.2399, 0.2092, 0.1989],
243 [ -0.0104, -0.2908, 0.1989, 0.3520, -0.0559, -0.1685, -0.0152, 0.2602, 0.0664],
244 [ 0.3520, -0.0563, -0.1684, -0.0563, 0.4362, -0.0841, -0.1684, -0.0841, 0.2118],
245 [ -0.0153, 0.2602, 0.0661, -0.1684, -0.0844, 0.2117, 0.2093, -0.0740, 0.0997],
246 [ -0.0306, 0.1991, 0.2399, -0.0155, -0.1683, 0.2092, 0.3112, 0.0666, 0.1989],
247 [ -0.0153, -0.1684, 0.2093, 0.2602, -0.0844, -0.0740, 0.0661, 0.2117, 0.0997],
248 [ 0.3113, 0.0663, 0.1991, 0.0663, 0.2117, 0.0993, 0.1991, 0.0993, 0.4770]])
249
250
251 r, theta, phi = cartesian_to_spherical([2, 1, 3])
252
253
254 Rx2_eigen = self.calc_Rx2_eigen_axis(theta, phi)
255
256
257 f2 = compile_2nd_matrix_free_rotor(self.f2_temp, Rx2_eigen)
258
259
260 print_frame_order_2nd_degree(real, "real")
261 print_frame_order_2nd_degree(f2, "calculated")
262 print_frame_order_2nd_degree(real-f2, "difference")
263
264
265 for i in range(9):
266 for j in range(9):
267 print("Element %s, %s; diff %s." % (i, j, f2[i, j] - real[i, j]))
268 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
269
270
272 """Check if compile_2nd_matrix_iso_cone() can return the identity matrix for disorder."""
273
274
275 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
276
277
278 f2 = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen, pi, pi)
279
280
281 print_frame_order_2nd_degree(self.I_disorder, "Identity for disorder")
282 print_frame_order_2nd_degree(f2, "Compiled frame order")
283
284
285 for i in range(9):
286 for j in range(9):
287 print("Element %s, %s." % (i, j))
288 self.assertAlmostEqual(f2[i, j], self.I_disorder[i, j])
289
290
292 """Check if compile_2nd_matrix_iso_cone() can return the matrix for a half cone."""
293
294
295 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
296
297
298 f2 = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen, pi/2.0, pi)
299
300
301 print_frame_order_2nd_degree(self.f2_half_cone, "The half cone frame order matrix")
302 print_frame_order_2nd_degree(f2, "Compiled frame order")
303
304
305 for i in range(9):
306 for j in range(9):
307 print("Element %s, %s." % (i, j))
308 self.assertAlmostEqual(f2[i, j], self.f2_half_cone[i, j])
309
310
312 """Check if compile_2nd_matrix_iso_cone() can return the matrix for a half cone rotated 90 degrees about y."""
313
314
315 Rx2_eigen = self.calc_Rx2_eigen_axis(pi/2.0, 0.0)
316
317
318 f2 = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen, pi/2.0, pi)
319
320
321 print_frame_order_2nd_degree(self.f2_half_cone_90_y, "The half cone frame order matrix")
322 print_frame_order_2nd_degree(f2, "Compiled frame order")
323
324
325 for i in range(9):
326 for j in range(9):
327 print("Element %s, %s." % (i, j))
328 self.assertAlmostEqual(f2[i, j], self.f2_half_cone_90_y[i, j])
329
330
332 """Check if compile_2nd_matrix_iso_cone() can return the identity matrix for order."""
333
334
335 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
336
337
338 f2 = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen, 1e-5, 1e-10)
339
340
341 print_frame_order_2nd_degree(self.I_order, "Identity for order")
342 print_frame_order_2nd_degree(f2, "Compiled frame order")
343
344
345 for i in range(9):
346 for j in range(9):
347 print("Element %s, %s." % (i, j))
348 self.assertAlmostEqual(f2[i, j], self.I_order[i, j])
349
350
352 """2nd check if compile_2nd_matrix_iso_cone() can return the identity matrix for order."""
353
354
355 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
356
357
358 f2 = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen, 0.0, 0.0)
359
360
361 print_frame_order_2nd_degree(self.I_order, "Identity for order")
362 print_frame_order_2nd_degree(f2, "Compiled frame order")
363
364
365 for i in range(9):
366 for j in range(9):
367 print("Element %s, %s." % (i, j))
368 self.assertAlmostEqual(f2[i, j], self.I_order[i, j])
369
370
372 """Check if compile_2nd_matrix_iso_cone() can approximate compile_2nd_matrix_iso_cone_free_rotor()."""
373
374
375 Rx2_eigen_a = self.calc_Rx2_eigen_axis(0.0, 0.0)
376 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 1.0)
377
378
379 f2a = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen_a, pi/4.6, pi)
380 f2b = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen_b, pi/4.6)
381
382
383 print_frame_order_2nd_degree(f2a, "Isotropic cone frame order")
384 print_frame_order_2nd_degree(f2b, "Free rotor isotropic cone frame order")
385 print_frame_order_2nd_degree(f2b-f2a, "difference")
386
387
388 for i in range(9):
389 for j in range(9):
390 print("Element %s, %s." % (i, j))
391 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
392
393
395 """Check if compile_2nd_matrix_iso_cone() can approximate compile_2nd_matrix_iso_cone_torsionless()."""
396
397
398 Rx2_eigen_a = self.calc_Rx2_eigen_axis(0.0, 0.0)
399 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 0.0)
400
401
402 f2a = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen_a, pi/4.6, 0)
403 f2b = compile_2nd_matrix_iso_cone_torsionless(self.f2_temp, Rx2_eigen_b, pi/4.6)
404
405
406 print_frame_order_2nd_degree(f2a, "Isotropic cone frame order")
407 print_frame_order_2nd_degree(f2b, "Torsionless isotropic cone frame order")
408 print_frame_order_2nd_degree(f2b-f2a, "difference")
409
410
411 for i in range(9):
412 for j in range(9):
413 print("Element %s, %s." % (i, j))
414 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
415
416
418 """Check if compile_2nd_matrix_iso_cone_free_rotor() can return the identity matrix for disorder."""
419
420
421 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 1.0)
422
423
424 f2 = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen, pi)
425
426
427 f2_ref = self.f2_half_cone
428
429
430 print_frame_order_2nd_degree(self.I_disorder, "Identity for disorder")
431 print_frame_order_2nd_degree(f2, "Compiled frame order")
432
433
434 for i in range(9):
435 for j in range(9):
436 print("Element %s, %s." % (i, j))
437 self.assertAlmostEqual(f2[i, j], self.I_disorder[i, j])
438
439
441 """Check if compile_2nd_matrix_iso_cone_free_rotor() can return the matrix for a half cone."""
442
443
444 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 1.0)
445
446
447 f2 = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen, pi/2.0)
448
449
450 print_frame_order_2nd_degree(self.f2_half_cone, "The half cone frame order matrix")
451 print_frame_order_2nd_degree(f2, "Compiled frame order")
452
453
454 for i in range(9):
455 for j in range(9):
456 print("Element %s, %s." % (i, j))
457 self.assertAlmostEqual(f2[i, j], self.f2_half_cone[i, j])
458
459
461 """Check if compile_2nd_matrix_iso_cone_free_rotor() can return the matrix for a half cone rotated 90 degrees about y."""
462
463
464 self.z_axis = array([1, 0, 0], float64)
465
466
467 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 1.0)
468
469
470 f2 = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen, pi/2.0)
471
472
473 print_frame_order_2nd_degree(self.f2_half_cone_90_y, "The half cone frame order matrix")
474 print_frame_order_2nd_degree(f2, "Compiled frame order")
475
476
477 for i in range(9):
478 for j in range(9):
479 print("Element %s, %s." % (i, j))
480 self.assertAlmostEqual(f2[i, j], self.f2_half_cone_90_y[i, j])
481
482
484 """Check if compile_2nd_matrix_iso_cone_free_rotor() can return the identity matrix for order."""
485
486
487 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 1.0)
488
489
490 f2 = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen, 0.0)
491
492
493 print_frame_order_2nd_degree(self.I_order_free_rotor, "Free rotor identity for order")
494 print_frame_order_2nd_degree(f2, "Compiled frame order")
495
496
497 for i in range(9):
498 for j in range(9):
499 print("Element %s, %s." % (i, j))
500 self.assertAlmostEqual(f2[i, j], self.I_order_free_rotor[i, j])
501
502
504 """Check the operation of the compile_2nd_matrix_pseudo_ellipse() function."""
505
506
507 real = array(
508 [[ 0.7901, 0, 0, 0, 0.7118, 0, 0, 0, 0.6851],
509 [ 0, 0.0816, 0, -0.0606, 0, 0, 0, 0, 0],
510 [ 0, 0, 0.1282, 0, 0, 0, -0.1224, 0, 0],
511 [ 0, -0.0606, 0, 0.0708, 0, 0, 0, 0, 0],
512 [ 0.7118, 0, 0, 0, 0.6756, 0, 0, 0, 0.6429],
513 [ 0, 0, 0, 0, 0, 0.2536, 0, -0.2421, 0],
514 [ 0, 0, -0.1224, 0, 0, 0, 0.1391, 0, 0],
515 [ 0, 0, 0, 0, 0, -0.2421, 0, 0.2427, 0],
516 [ 0.6851, 0, 0, 0, 0.6429, 0, 0, 0, 0.6182]], float64)
517 transpose_23(real)
518
519
520 x = pi/4.0
521 y = 3.0*pi/8.0
522 z = pi/6.0
523
524
525 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
526
527
528 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, x, y, z)
529
530
531 print_frame_order_2nd_degree(real, "real")
532 print_frame_order_2nd_degree(f2, "calculated")
533 print_frame_order_2nd_degree(real-f2, "difference")
534
535
536 for i in range(9):
537 for j in range(9):
538 print("Element %s, %s; diff %s." % (i, j, f2[i, j] - real[i, j]))
539 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-4)
540
541
543 """Check the operation of the compile_2nd_matrix_pseudo_ellipse() function."""
544
545
546 real = array(
547 [[ 0.7379, 0, 0, 0, 0.1338, 0, 0, 0, 0.1284],
548 [ 0, 0.6637, 0, -0.1085, 0, 0, 0, 0, 0],
549 [ 0, 0, 0.6603, 0, 0, 0, -0.1181, 0, 0],
550 [ 0, -0.1085, 0, 0.6637, 0, 0, 0, 0, 0],
551 [ 0.1154, 0, 0, 0, 0.6309, 0, 0, 0, 0.2536],
552 [ 0, 0, 0, 0, 0, 0.6196, 0, -0.2336, 0],
553 [ 0, 0, -0.1181, 0, 0, 0, 0.6603, 0, 0],
554 [ 0, 0, 0, 0, 0, -0.2336, 0, 0.6196, 0],
555 [ 0.1467, 0, 0, 0, 0.2353, 0, 0, 0, 0.6180]], float64)
556
557
558 x = pi/4.0
559 y = 3.0*pi/8.0
560 z = 40.0 / 360.0 * 2.0 * pi
561
562
563 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
564
565
566 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, x, y, z)
567
568
569 print_frame_order_2nd_degree(real, "real")
570 print_frame_order_2nd_degree(f2, "calculated")
571 print_frame_order_2nd_degree(real-f2, "difference")
572
573
574 for i in range(9):
575 for j in range(9):
576 print("Element %s, %s; diff %s." % (i, j, f2[i, j] - real[i, j]))
577 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
578
579
581 """Check the operation of the compile_2nd_matrix_pseudo_ellipse() function."""
582
583
584 real = array(
585 [[ 0.6314, 0.0232, -0.0344, 0.0232, 0.1558, -0.0222, -0.0344, -0.0222, 0.2128],
586 [ 0.0220, 0.6366, 0.0069, -0.1352, 0.0243, -0.0722, 0.0206, -0.0277, -0.0464],
587 [ -0.0332, 0.0097, 0.6137, 0.0222, 0.0668, 0.0173, -0.1967, 0.0489, -0.0336],
588 [ 0.0220, -0.1352, 0.0206, 0.6366, 0.0243, -0.0277, 0.0069, -0.0722, -0.0464],
589 [ 0.1554, 0.0233, 0.0669, 0.0233, 0.6775, 0.0113, 0.0669, 0.0113, 0.1671],
590 [ -0.0222, -0.0738, 0.0188, -0.0286, 0.0113, 0.6310, 0.0507, -0.1502, 0.0109],
591 [ -0.0332, 0.0222, -0.1967, 0.0097, 0.0668, 0.0489, 0.6137, 0.0173, -0.0336],
592 [ -0.0222, -0.0286, 0.0507, -0.0738, 0.0113, -0.1502, 0.0188, 0.6310, 0.0109],
593 [ 0.2132, -0.0465, -0.0324, -0.0465, 0.1667, 0.0110, -0.0324, 0.0110, 0.6201]], float64)
594
595
596 x = 60.0 / 360.0 * 2.0 * pi
597 y = 3.0 * pi / 8.0
598 z = pi / 6.0
599
600
601 Rx2_eigen = self.calc_Rx2_eigen_full(self.out_of_frame_alpha, self.out_of_frame_beta, self.out_of_frame_gamma)
602
603
604 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, x, y, z)
605
606
607 print_frame_order_2nd_degree(real, "real")
608 print_frame_order_2nd_degree(f2, "calculated")
609 print_frame_order_2nd_degree(real-f2, "difference")
610
611
612 for i in range(9):
613 for j in range(9):
614 print("Element %s, %s; diff %s." % (i, j, f2[i, j] - real[i, j]))
615 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
616
617
619 """Check if compile_2nd_matrix_pseudo_ellipse() can return the identity matrix for disorder."""
620
621
622 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
623
624
625 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, pi, pi, pi)
626
627
628 print_frame_order_2nd_degree(self.I_disorder, "Identity for disorder")
629 print_frame_order_2nd_degree(f2, "Compiled frame order")
630
631
632 for i in range(9):
633 for j in range(9):
634 print("Element %s, %s." % (i, j))
635 self.assertAlmostEqual(f2[i, j], self.I_disorder[i, j])
636
637
639 """Check if compile_2nd_matrix_pseudo_ellipse() can return the matrix for a half cone."""
640
641
642 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
643
644
645 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, pi/2.0, pi/2.0, pi)
646
647
648 print_frame_order_2nd_degree(self.f2_half_cone, "The half cone frame order matrix")
649 print_frame_order_2nd_degree(f2, "Compiled frame order")
650
651
652 for i in range(9):
653 for j in range(9):
654 print("Element %s, %s." % (i, j))
655 self.assertAlmostEqual(f2[i, j], self.f2_half_cone[i, j])
656
657
659 """Check if compile_2nd_matrix_pseudo_ellipse() can return the matrix for a half cone rotated 90 degrees about y."""
660
661
662 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, pi/2.0, 0.0)
663
664
665 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, pi/2.0, pi/2.0, pi)
666
667
668 print_frame_order_2nd_degree(self.f2_half_cone_90_y, "The half cone frame order matrix")
669 print_frame_order_2nd_degree(f2, "Compiled frame order")
670
671
672 for i in range(9):
673 for j in range(9):
674 print("Element %s, %s." % (i, j))
675 self.assertAlmostEqual(f2[i, j], self.f2_half_cone_90_y[i, j])
676
677
679 """Check if compile_2nd_matrix_pseudo_ellipse() can return the identity matrix for order."""
680
681
682 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
683
684
685 f2 = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen, 1e-2, 1e-2, 1e-10)
686
687
688 print_frame_order_2nd_degree(self.I_order, "Identity for order")
689 print_frame_order_2nd_degree(f2, "Compiled frame order")
690 print_frame_order_2nd_degree(f2-self.I_order, "difference")
691
692
693 for i in range(9):
694 for j in range(9):
695 print("Element %s, %s." % (i, j))
696 self.assertAlmostEqual(f2[i, j], self.I_order[i, j], 4)
697
698
700 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate compile_2nd_matrix_pseudo_ellipse_free_rotor()."""
701
702
703 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
704 Rx2_eigen_b = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
705
706
707 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/1.6, pi/5.8, pi)
708 f2b = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen_b, pi/1.6, pi/5.8)
709
710
711 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
712 print_frame_order_2nd_degree(f2b, "Free rotor pseudo-ellipse frame order")
713 print_frame_order_2nd_degree(f2b-f2a, "difference")
714
715
716 for i in range(9):
717 for j in range(9):
718 print("Element %s, %s." % (i, j))
719 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
720
721
723 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate a pi/2 rotated compile_2nd_matrix_pseudo_ellipse_free_rotor()."""
724
725
726 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
727 Rx2_eigen_b = self.calc_Rx2_eigen_full(pi/2, 0.0, 0.0)
728
729
730 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/1.6, pi/5.8, pi)
731 f2b = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen_b, pi/5.8, pi/1.6)
732
733
734 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
735 print_frame_order_2nd_degree(f2b, "pi/2 rotated free rotor pseudo-ellipse frame order")
736 print_frame_order_2nd_degree(f2b-f2a, "difference")
737
738
739 for i in range(9):
740 for j in range(9):
741 print("Element %s, %s." % (i, j))
742 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
743
744
746 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate compile_2nd_matrix_pseudo_ellipse_torsionless()."""
747
748
749 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
750 Rx2_eigen_b = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
751
752
753 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/2.1, pi/4.6, 0)
754 f2b = compile_2nd_matrix_pseudo_ellipse_torsionless(self.f2_temp, Rx2_eigen_b, pi/2.1, pi/4.6)
755
756
757 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
758 print_frame_order_2nd_degree(f2b, "Torsionless pseudo-ellipse frame order")
759 print_frame_order_2nd_degree(f2b-f2a, "difference")
760
761
762 for i in range(9):
763 for j in range(9):
764 print("Element %s, %s; diff %s." % (i, j, f2b[i, j] - f2a[i, j]))
765 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
766
767
769 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate compile_2nd_matrix_iso_cone()."""
770
771
772 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
773 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 0.0)
774
775
776 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/4.6, pi/4.6, 0.2)
777 f2b = compile_2nd_matrix_iso_cone(self.f2_temp, Rx2_eigen_b, pi/4.6, 0.2)
778
779
780 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
781 print_frame_order_2nd_degree(f2b, "Isotropic cone frame order")
782 print_frame_order_2nd_degree(f2b-f2a, "difference")
783
784
785 for i in range(9):
786 for j in range(9):
787 print("Element %s, %s." % (i, j))
788 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
789
790
792 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate compile_2nd_matrix_iso_cone_free_rotor()."""
793
794
795 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
796 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 1.0)
797
798
799 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/4.6, pi/4.6, pi)
800 f2b = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen_b, pi/4.6)
801
802
803 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
804 print_frame_order_2nd_degree(f2b, "Free rotor isotropic cone frame order")
805 print_frame_order_2nd_degree(f2b-f2a, "difference")
806
807
808 for i in range(9):
809 for j in range(9):
810 print("Element %s, %s." % (i, j))
811 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
812
813
815 """Check if compile_2nd_matrix_pseudo_ellipse() can approximate compile_2nd_matrix_iso_cone_torsionless()."""
816
817
818 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
819 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 0.0)
820
821
822 f2a = compile_2nd_matrix_pseudo_ellipse(self.f2_temp, Rx2_eigen_a, pi/8.6, pi/8.6, 0)
823 f2b = compile_2nd_matrix_iso_cone_torsionless(self.f2_temp, Rx2_eigen_b, pi/8.6)
824
825
826 print_frame_order_2nd_degree(f2a, "Pseudo-ellipse frame order")
827 print_frame_order_2nd_degree(f2b, "Torsionless isotropic cone frame order")
828 print_frame_order_2nd_degree(f2b-f2a, "difference")
829
830
831 for i in range(9):
832 for j in range(9):
833 print("Element %s, %s." % (i, j))
834 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
835
836
838 """Check if compile_2nd_matrix_pseudo_ellipse_free_rotor() can return the identity matrix for disorder."""
839
840
841 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
842
843
844 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, pi, pi)
845
846
847 print_frame_order_2nd_degree(self.I_disorder, "Identity for disorder")
848 print_frame_order_2nd_degree(f2, "Compiled frame order")
849
850
851 for i in range(9):
852 for j in range(9):
853 print("Element %s, %s." % (i, j))
854 self.assertAlmostEqual(f2[i, j], self.I_disorder[i, j])
855
856
858 """Check if compile_2nd_matrix_pseudo_ellipse() can return the matrix for a half cone."""
859
860
861 Rx2_eigen = self.calc_Rx2_eigen_full(pi, 0.0, pi)
862
863
864 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, pi/2.0, pi/2.0)
865
866
867 print_frame_order_2nd_degree(self.f2_half_cone, "The half cone frame order matrix")
868 print_frame_order_2nd_degree(f2, "Compiled frame order")
869
870
871 for i in range(9):
872 for j in range(9):
873 print("Element %s, %s." % (i, j))
874 self.assertAlmostEqual(f2[i, j], self.f2_half_cone[i, j])
875
876
878 """Check if compile_2nd_matrix_pseudo_ellipse() can return the matrix for a half cone rotated 90 degrees about y."""
879
880
881 Rx2_eigen = self.calc_Rx2_eigen_full(pi, pi/2.0, pi)
882
883
884 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, pi/2.0, pi/2.0)
885
886
887 print_frame_order_2nd_degree(self.f2_half_cone_90_y, "The half cone frame order matrix")
888 print_frame_order_2nd_degree(f2, "Compiled frame order")
889
890
891 for i in range(9):
892 for j in range(9):
893 print("Element %s, %s." % (i, j))
894 self.assertAlmostEqual(f2[i, j], self.f2_half_cone_90_y[i, j])
895
896
898 """Check if compile_2nd_matrix_pseudo_ellipse_free_rotor() can return the identity matrix for order."""
899
900
901 Rx2_eigen = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
902
903
904 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, 1e-10, 1e-10)
905
906
907 print_frame_order_2nd_degree(self.I_order_free_rotor, "Free rotor identity for order")
908 print_frame_order_2nd_degree(f2, "Compiled frame order")
909
910
911 for i in range(9):
912 for j in range(9):
913 print("Element %s, %s." % (i, j))
914 self.assertAlmostEqual(f2[i, j], self.I_order_free_rotor[i, j])
915
916
918 """Check the operation of the compile_2nd_matrix_pseudo_ellipse_free_rotor() function."""
919
920
921 real = array(
922 [[ 0.3428, -0.0193, 0.0389, -0.0193, 0.3137, -0.0194, 0.0389, -0.0194, 0.3435],
923 [ -0.0225, 0.2313, 0.0034, -0.1413, 0.0449, 0.2309, -0.1830, -0.1412, -0.0224],
924 [ 0.0417, 0.0091, 0.2142, -0.1767, -0.0838, 0.0092, 0.1211, -0.1770, 0.0421],
925 [ -0.0225, -0.1413, -0.1830, 0.2313, 0.0449, -0.1412, 0.0034, 0.2309, -0.0224],
926 [ 0.3124, 0.0418, -0.0840, 0.0418, 0.3758, 0.0418, -0.0840, 0.0418, 0.3118],
927 [ -0.0193, 0.2251, 0.0151, -0.1476, 0.0389, 0.2251, -0.1706, -0.1468, -0.0196],
928 [ 0.0417, -0.1767, 0.1211, 0.0091, -0.0838, -0.1770, 0.2142, 0.0092, 0.0421],
929 [ -0.0193, -0.1476, -0.1706, 0.2251, 0.0389, -0.1468, 0.0151, 0.2251, -0.0196],
930 [ 0.3448, -0.0225, 0.0450, -0.0225, 0.3104, -0.0224, 0.0450, -0.0224, 0.3447]], float64)
931
932
933 x = pi/4.0
934 y = 50.0 / 360.0 * 2.0 * pi
935
936
937 Rx2_eigen = self.calc_Rx2_eigen_full(self.out_of_frame_alpha, self.out_of_frame_beta, self.out_of_frame_gamma)
938
939
940 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, x, y)
941
942
943 print_frame_order_2nd_degree(real, "real")
944 print_frame_order_2nd_degree(f2, "calculated")
945 print_frame_order_2nd_degree(real-f2, "difference")
946
947
948 for i in range(9):
949 for j in range(9):
950 print("Element %s, %s." % (i, j))
951 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
952
953
955 """Check the operation of the compile_2nd_matrix_pseudo_ellipse_free_rotor() function."""
956
957
958 real = array(
959 [[ 0.3251, 0.0163, -0.0324, 0.0163, 0.3493, 0.0164, -0.0324, 0.0164, 0.3256],
960 [ -0.0248, 0.1481, -0.0480, -0.0500, 0.0492, 0.1475, -0.1472, -0.0500, -0.0244],
961 [ 0.0079, 0.0328, 0.0572, -0.0661, -0.0163, 0.0331, 0.0074, -0.0662, 0.0084],
962 [ -0.0248, -0.0500, -0.1472, 0.1481, 0.0492, -0.0500, -0.0480, 0.1475, -0.0244],
963 [ 0.3289, 0.0081, -0.0167, 0.0081, 0.3426, 0.0080, -0.0167, 0.0080, 0.3285],
964 [ 0.0163, 0.0669, 0.1139, -0.1322, -0.0324, 0.0662, 0.0157, -0.1307, 0.0161],
965 [ 0.0079, -0.0661, 0.0074, 0.0328, -0.0163, -0.0662, 0.0572, 0.0331, 0.0084],
966 [ 0.0163, -0.1322, 0.0157, 0.0669, -0.0324, -0.1307, 0.1139, 0.0662, 0.0161],
967 [ 0.3459, -0.0245, 0.0491, -0.0245, 0.3081, -0.0244, 0.0491, -0.0244, 0.3460]], float64)
968
969
970 x = pi / 4.0
971 y = 150.0 / 360.0 * 2.0 * pi
972
973
974 Rx2_eigen = self.calc_Rx2_eigen_full(self.out_of_frame_alpha, self.out_of_frame_beta, self.out_of_frame_gamma)
975
976
977 f2 = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen, x, y)
978
979
980 print_frame_order_2nd_degree(real, "real")
981 print_frame_order_2nd_degree(f2, "calculated")
982 print_frame_order_2nd_degree(real-f2, "difference")
983
984
985 for i in range(9):
986 for j in range(9):
987 print("Element %s, %s." % (i, j))
988 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-2)
989
990
992 """Check if compile_2nd_matrix_pseudo_ellipse_free_rotor() can approximate compile_2nd_matrix_iso_cone_free_rotor()."""
993
994
995 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
996 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 1.0)
997
998
999 f2a = compile_2nd_matrix_pseudo_ellipse_free_rotor(self.f2_temp, Rx2_eigen_a, pi/4.6, pi/4.6)
1000 f2b = compile_2nd_matrix_iso_cone_free_rotor(self.f2_temp, Rx2_eigen_b, pi/4.6)
1001
1002
1003 print_frame_order_2nd_degree(f2a, "Free rotor pseudo-ellipse frame order")
1004 print_frame_order_2nd_degree(f2b, "Free rotor isotropic cone frame order")
1005 print_frame_order_2nd_degree(f2b-f2a, "difference")
1006
1007
1008 for i in range(9):
1009 for j in range(9):
1010 print("Element %s, %s." % (i, j))
1011 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
1012
1013
1015 """Check if compile_2nd_matrix_pseudo_ellipse_torsionless() can approximate compile_2nd_matrix_iso_cone_torsionless()."""
1016
1017
1018 Rx2_eigen_a = self.calc_Rx2_eigen_full(0.0, 0.0, 0.0)
1019 Rx2_eigen_b = self.calc_Rx2_eigen_axis(0.0, 0.0)
1020
1021
1022 f2a = compile_2nd_matrix_pseudo_ellipse_torsionless(self.f2_temp, Rx2_eigen_a, pi/4.6, pi/4.6)
1023 f2b = compile_2nd_matrix_iso_cone_torsionless(self.f2_temp, Rx2_eigen_b, pi/4.6)
1024
1025
1026 print_frame_order_2nd_degree(f2a, "Torsionless pseudo-ellipse frame order")
1027 print_frame_order_2nd_degree(f2b, "Torsionless isotropic cone frame order")
1028 print_frame_order_2nd_degree(f2b-f2a, "difference")
1029
1030
1031 for i in range(9):
1032 for j in range(9):
1033 print("Element %s, %s." % (i, j))
1034 self.assertAlmostEqual(f2a[i, j], f2b[i, j])
1035
1036
1038 """Check the operation of the compile_2nd_matrix_rotor() function."""
1039
1040
1041 real = array(
1042 [[ 7.06775425e-01, 1.36710179e-04, 0.00000000e+00, 1.36710179e-04, 2.93224575e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1043 [ -1.36710179e-04, 7.06775425e-01, 0.00000000e+00, -2.93224575e-01, 1.36710179e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1044 [ 0.00000000e+00, 0.00000000e+00, 8.27014112e-01, 0.00000000e+00, 0.00000000e+00, 2.19539417e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1045 [ -1.36710179e-04, -2.93224575e-01, 0.00000000e+00, 7.06775425e-01, 1.36710179e-04, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1046 [ 2.93224575e-01, -1.36710179e-04, 0.00000000e+00, -1.36710179e-04, 7.06775425e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1047 [ 0.00000000e+00, 0.00000000e+00, -2.19539417e-04, 0.00000000e+00, 0.00000000e+00, 8.27014112e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
1048 [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.27014112e-01, 2.19539417e-04, 0.00000000e+00],
1049 [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -2.19539417e-04, 8.27014112e-01, 0.00000000e+00],
1050 [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])
1051
1052
1053 sigma_max = 60.0 / 180.0 * pi
1054
1055
1056 Rx2_eigen = self.calc_Rx2_eigen_axis(0.0, 0.0)
1057
1058
1059 f2 = compile_2nd_matrix_rotor(self.f2_temp, Rx2_eigen, sigma_max)
1060
1061
1062 print_frame_order_2nd_degree(real, "real")
1063 print_frame_order_2nd_degree(f2, "calculated")
1064 print_frame_order_2nd_degree(real-f2, "difference")
1065
1066
1067 for i in range(9):
1068 for j in range(9):
1069 print("Element %s, %s." % (i, j))
1070 self.assertTrue(abs(f2[i, j] - real[i, j]) < 1e-3)
1071
1072
1074 """Test the alignment tensor reduction for the order identity matrix."""
1075
1076
1077 A = array([1, 2, 3, 4, 5], float64)
1078 red = zeros(5, float64)
1079
1080
1081 reduce_alignment_tensor(self.I_order, A, red)
1082
1083
1084 for i in range(5):
1085 self.assertEqual(A[i], red[i])
1086
1087
1089 """Test the alignment tensor reduction for the order identity matrix."""
1090
1091
1092 A = array([1, 2, 3, 4, 5], float64)
1093 red = zeros(5, float64)
1094
1095
1096 reduce_alignment_tensor(self.I_disorder, A, red)
1097
1098
1099 for i in range(5):
1100 self.assertEqual(red[i], 0.0)
1101
1102
1104 """Test the alignment tensor reduction for the order identity matrix."""
1105
1106
1107 A = array([1, 2, 3, 4, 5], float64)
1108 red = zeros(5, float64)
1109
1110
1111 reduce_alignment_tensor(self.f2_half_cone, A, red)
1112
1113
1114 for i in range(5):
1115 self.assertEqual(red[i], 0.0)
1116