lib.frame_order.pseudo

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for the pseudo-ellipse frame order model.""" 24 25 # Python module imports. 26 from math import cos, pi, sin, sqrt 27 from numpy import divide, dot, eye, float64, multiply, sinc, swapaxes, tensordot 28 from numpy import cos as np_cos 29 from numpy import sin as np_sin 30 from numpy import sqrt as np_sqrt 31 try: 32 from scipy.integrate import quad, tplquad 33 except ImportError: 34 pass 35 36 # relax module imports. 37 from lib.geometry.pec import pec 38 from lib.frame_order.matrix_ops import pcs_pivot_motion_full_qr_int, pcs_pivot_motion_full_quad_int, rotate_daeg 39 40

41 -def compile_1st_matrix_pseudo_ellipse(matrix, R_eigen, theta_x, theta_y, sigma_max):

42 """Generate the 1st degree Frame Order matrix for the pseudo-ellipse. 43 44 @param matrix: The Frame Order matrix, 1st degree to be populated. 45 @type matrix: numpy 3D, rank-2 array 46 @param R_eigen: The eigenframe rotation matrix. 47 @type R_eigen: numpy 3D, rank-2 array 48 @param theta_x: The cone opening angle along x. 49 @type theta_x: float 50 @param theta_y: The cone opening angle along y. 51 @type theta_y: float 52 @param sigma_max: The maximum torsion angle. 53 @type sigma_max: float 54 """ 55 56 # The surface area normalisation factor. 57 fact = 2.0 * pec(theta_x, theta_y) 58 59 # Invert. 60 if fact == 0.0: 61 fact = 1e100 62 else: 63 fact = 1.0 / fact 64 65 # Sinc value. 66 sinc_smax = sinc(sigma_max/pi) 67 68 # Numerical integration of phi of each element. 69 matrix[0, 0] = fact * sinc_smax * (2.0*pi + quad(part_int_daeg1_pseudo_ellipse_00, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 70 matrix[1, 1] = fact * sinc_smax * (2.0*pi + quad(part_int_daeg1_pseudo_ellipse_11, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 71 matrix[2, 2] = fact * quad(part_int_daeg1_pseudo_ellipse_22, -pi, pi, args=(theta_x, theta_y), full_output=1)[0] 72 73 # Rotate and return the frame order matrix. 74 return rotate_daeg(matrix, R_eigen)

75 76

77 -def compile_2nd_matrix_pseudo_ellipse(matrix, Rx2_eigen, theta_x, theta_y, sigma_max):

78 """Generate the 2nd degree Frame Order matrix for the pseudo-ellipse. 79 80 @param matrix: The Frame Order matrix, 2nd degree to be populated. 81 @type matrix: numpy 9D, rank-2 array 82 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 83 @type Rx2_eigen: numpy 9D, rank-2 array 84 @param theta_x: The cone opening angle along x. 85 @type theta_x: float 86 @param theta_y: The cone opening angle along y. 87 @type theta_y: float 88 @param sigma_max: The maximum torsion angle. 89 @type sigma_max: float 90 """ 91 92 # The rigid case. 93 if theta_x == 0.0 and sigma_max == 0.0: 94 # Set up the matrix as the identity. 95 matrix[:] = 0.0 96 for i in range(len(matrix)): 97 matrix[i, i] = 1.0 98 99 # Rotate and return the frame order matrix. 100 return rotate_daeg(matrix, Rx2_eigen) 101 102 # The surface area normalisation factor. 103 fact = 12.0 * pec(theta_x, theta_y) 104 105 # Invert. 106 if fact == 0.0: 107 fact = 1e100 108 else: 109 fact = 1.0 / fact 110 111 # Repetitive calculations. 112 sinc_smax = sinc(sigma_max/pi) 113 sinc_2smax = sinc(2.0*sigma_max/pi) 114 115 # Diagonal. 116 matrix[0, 0] = fact * (4.0*pi*(sinc_2smax + 2.0) + quad(part_int_daeg2_pseudo_ellipse_00, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 117 matrix[1, 1] = fact * (4.0*pi*sinc_2smax + quad(part_int_daeg2_pseudo_ellipse_11, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 118 matrix[2, 2] = fact * 2.0*sinc_smax * (5.0*pi - quad(part_int_daeg2_pseudo_ellipse_22, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 119 matrix[3, 3] = matrix[1, 1] 120 matrix[4, 4] = fact * (4.0*pi*(sinc_2smax + 2.0) + quad(part_int_daeg2_pseudo_ellipse_44, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 121 matrix[5, 5] = fact * 2.0*sinc_smax * (5.0*pi - quad(part_int_daeg2_pseudo_ellipse_55, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 122 matrix[6, 6] = matrix[2, 2] 123 matrix[7, 7] = matrix[5, 5] 124 matrix[8, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_88, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 125 126 # Off diagonal set 1. 127 matrix[0, 4] = fact * (4.0*pi*(2.0 - sinc_2smax) + quad(part_int_daeg2_pseudo_ellipse_04, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 128 matrix[4, 0] = fact * (4.0*pi*(2.0 - sinc_2smax) + quad(part_int_daeg2_pseudo_ellipse_40, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 129 matrix[0, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_08, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 130 matrix[8, 0] = fact * (8.0*pi + quad(part_int_daeg2_pseudo_ellipse_80, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 131 matrix[4, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_48, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 132 matrix[8, 4] = fact * (8.0*pi - quad(part_int_daeg2_pseudo_ellipse_84, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 133 134 # Off diagonal set 2. 135 matrix[1, 3] = matrix[3, 1] = fact * (4.0*pi*sinc_2smax + quad(part_int_daeg2_pseudo_ellipse_13, -pi, pi, args=(theta_x, theta_y, sinc_2smax), full_output=1)[0]) 136 matrix[2, 6] = matrix[6, 2] = -fact * 4.0 * sinc_smax * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_26, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 137 matrix[5, 7] = matrix[7, 5] = -fact * 4.0 * sinc_smax * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_57, -pi, pi, args=(theta_x, theta_y), full_output=1)[0]) 138 139 # Rotate and return the frame order matrix. 140 return rotate_daeg(matrix, Rx2_eigen)

141 142

143 -def part_int_daeg1_pseudo_ellipse_00(phi, x, y):

144 """The theta-sigma partial integral of the 1st degree Frame Order matrix element xx for the pseudo-ellipse. 145 146 @param phi: The azimuthal tilt-torsion angle. 147 @type phi: float 148 @param x: The cone opening angle along x. 149 @type x: float 150 @param y: The cone opening angle along y. 151 @type y: float 152 @return: The theta-sigma partial integral. 153 @rtype: float 154 """ 155 156 # Theta max. 157 tmax = tmax_pseudo_ellipse(phi, x, y) 158 159 # The theta-sigma integral. 160 return cos(phi)**2 * sin(tmax)**2 - 2.0 * sin(phi)**2 * cos(tmax)

161 162

163 -def part_int_daeg1_pseudo_ellipse_11(phi, x, y):

164 """The theta-sigma partial integral of the 1st degree Frame Order matrix element yy for the pseudo-ellipse. 165 166 @param phi: The azimuthal tilt-torsion angle. 167 @type phi: float 168 @param x: The cone opening angle along x. 169 @type x: float 170 @param y: The cone opening angle along y. 171 @type y: float 172 @return: The theta-sigma partial integral. 173 @rtype: float 174 """ 175 176 # Theta max. 177 tmax = tmax_pseudo_ellipse(phi, x, y) 178 179 # The theta-sigma integral. 180 return sin(phi)**2 * sin(tmax)**2 - 2.0 * cos(phi)**2 * cos(tmax)

181 182

183 -def part_int_daeg1_pseudo_ellipse_22(phi, x, y):

184 """The theta-sigma partial integral of the 1st degree Frame Order matrix element zz for the pseudo-ellipse. 185 186 @param phi: The azimuthal tilt-torsion angle. 187 @type phi: float 188 @param x: The cone opening angle along x. 189 @type x: float 190 @param y: The cone opening angle along y. 191 @type y: float 192 @return: The theta-sigma partial integral. 193 @rtype: float 194 """ 195 196 # Theta max. 197 tmax = tmax_pseudo_ellipse(phi, x, y) 198 199 # The theta-sigma integral. 200 return sin(tmax)**2

201 202

203 -def part_int_daeg2_pseudo_ellipse_00(phi, x, y, sinc_2smax):

204 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 11 for the pseudo-ellipse. 205 206 @param phi: The azimuthal tilt-torsion angle. 207 @type phi: float 208 @param x: The cone opening angle along x. 209 @type x: float 210 @param y: The cone opening angle along y. 211 @type y: float 212 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 213 @type sinc_2smax: float 214 @return: The theta-sigma partial integral. 215 @rtype: float 216 """ 217 218 # Theta max. 219 tmax = tmax_pseudo_ellipse(phi, x, y) 220 221 # Repetitive trig. 222 cos_tmax = cos(tmax) 223 sin_tmax2 = sin(tmax)**2 224 cos_phi2 = cos(phi)**2 225 226 # Components. 227 a = sinc_2smax * (2.0 * sin_tmax2 * cos_phi2 * ((2.0*cos_phi2 - 1.0)*cos_tmax - 6.0*(cos_phi2 - 1.0)) - 2.0*cos_tmax * (2.0*cos_phi2*(4.0*cos_phi2 - 5.0) + 3.0)) 228 b = 2.0*cos_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 229 230 # Return the theta-sigma integral. 231 return a + b

232 233

234 -def part_int_daeg2_pseudo_ellipse_04(phi, x, y, sinc_2smax):

235 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 22 for the pseudo-ellipse. 236 237 @param phi: The azimuthal tilt-torsion angle. 238 @type phi: float 239 @param x: The cone opening angle along x. 240 @type x: float 241 @param y: The cone opening angle along y. 242 @type y: float 243 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 244 @type sinc_2smax: float 245 @return: The theta-sigma partial integral. 246 @rtype: float 247 """ 248 249 # Theta max. 250 tmax = tmax_pseudo_ellipse(phi, x, y) 251 252 # Repetitive trig. 253 cos_tmax = cos(tmax) 254 sin_tmax2 = sin(tmax)**2 255 cos_phi2 = cos(phi)**2 256 sin_phi2 = sin(phi)**2 257 258 # Components. 259 a = sinc_2smax * (2.0 * sin_tmax2 * cos_phi2 * ((2.0*sin_phi2 - 1.0)*cos_tmax - 6.0*sin_phi2) + 2.0*cos_tmax * (2.0*cos_phi2*(4.0*cos_phi2 - 5.0) + 3.0)) 260 b = 2.0*cos_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 261 262 # Return the theta-sigma integral. 263 return a + b

264 265

266 -def part_int_daeg2_pseudo_ellipse_08(phi, x, y):

267 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 33 for the pseudo-ellipse. 268 269 @param phi: The azimuthal tilt-torsion angle. 270 @type phi: float 271 @param x: The cone opening angle along x. 272 @type x: float 273 @param y: The cone opening angle along y. 274 @type y: float 275 @return: The theta-sigma partial integral. 276 @rtype: float 277 """ 278 279 # Theta max. 280 tmax = tmax_pseudo_ellipse(phi, x, y) 281 282 # The theta-sigma integral. 283 return cos(tmax) * cos(phi)**2 * (sin(tmax)**2 + 2.0)

284 285

286 -def part_int_daeg2_pseudo_ellipse_11(phi, x, y, sinc_2smax):

287 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 288 289 @param phi: The azimuthal tilt-torsion angle. 290 @type phi: float 291 @param x: The cone opening angle along x. 292 @type x: float 293 @param y: The cone opening angle along y. 294 @type y: float 295 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 296 @type sinc_2smax: float 297 @return: The theta-sigma partial integral. 298 @rtype: float 299 """ 300 301 # Theta max. 302 tmax = tmax_pseudo_ellipse(phi, x, y) 303 304 # Repetitive trig. 305 cos_tmax = cos(tmax) 306 sin_tmax2 = sin(tmax)**2 307 cos_phi2 = cos(phi)**2 308 sin_phi2 = sin(phi)**2 309 cos_sin_phi2 = cos_phi2 * sin_phi2 310 311 # The integral. 312 a = sinc_2smax * ((4.0*cos_sin_phi2*(cos_tmax - 3.0) + 3.0)*sin_tmax2 - 16.0*cos_sin_phi2*cos_tmax) + 3.0*sin_tmax2 313 314 # The theta-sigma integral. 315 return a

316 317

318 -def part_int_daeg2_pseudo_ellipse_13(phi, x, y, sinc_2smax):

319 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 320 321 @param phi: The azimuthal tilt-torsion angle. 322 @type phi: float 323 @param x: The cone opening angle along x. 324 @type x: float 325 @param y: The cone opening angle along y. 326 @type y: float 327 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 328 @type sinc_2smax: float 329 @return: The theta-sigma partial integral. 330 @rtype: float 331 """ 332 333 # Theta max. 334 tmax = tmax_pseudo_ellipse(phi, x, y) 335 336 # Repetitive trig. 337 cos_tmax = cos(tmax) 338 sin_tmax2 = sin(tmax)**2 339 cos_sin_phi2 = cos(phi)**2 * sin(phi)**2 340 341 # The theta-sigma integral. 342 return sinc_2smax * (sin_tmax2 * (4.0*cos_sin_phi2*cos_tmax - 12.0*cos_sin_phi2 + 3) - 16.0*cos_sin_phi2*cos_tmax) - 3.0*sin_tmax2

343 344

345 -def part_int_daeg2_pseudo_ellipse_22(phi, x, y):

346 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 347 348 @param phi: The azimuthal tilt-torsion angle. 349 @type phi: float 350 @param x: The cone opening angle along x. 351 @type x: float 352 @param y: The cone opening angle along y. 353 @type y: float 354 @return: The theta-sigma partial integral. 355 @rtype: float 356 """ 357 358 # Theta max. 359 tmax = tmax_pseudo_ellipse(phi, x, y) 360 361 # Repetitive trig. 362 cos_tmax = cos(tmax) 363 cos_phi2 = cos(phi)**2 364 sin_phi2 = sin(phi)**2 365 366 # Components. 367 a = 2.0*cos_phi2*cos_tmax**3 + 3.0*sin_phi2*cos_tmax**2 368 369 # Return the theta-sigma integral. 370 return a

371 372

373 -def part_int_daeg2_pseudo_ellipse_26(phi, x, y):

374 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 375 376 @param phi: The azimuthal tilt-torsion angle. 377 @type phi: float 378 @param x: The cone opening angle along x. 379 @type x: float 380 @param y: The cone opening angle along y. 381 @type y: float 382 @return: The theta-sigma partial integral. 383 @rtype: float 384 """ 385 386 # Theta max. 387 tmax = tmax_pseudo_ellipse(phi, x, y) 388 389 # Repetitive trig. 390 cos_tmax = cos(tmax) 391 392 # The theta-sigma integral. 393 return cos(phi)**2 * (cos_tmax**3 - 3.0*cos_tmax)

394 395

396 -def part_int_daeg2_pseudo_ellipse_40(phi, x, y, sinc_2smax):

397 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 398 399 @param phi: The azimuthal tilt-torsion angle. 400 @type phi: float 401 @param x: The cone opening angle along x. 402 @type x: float 403 @param y: The cone opening angle along y. 404 @type y: float 405 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 406 @type sinc_2smax: float 407 @return: The theta-sigma partial integral. 408 @rtype: float 409 """ 410 411 # Theta max. 412 tmax = tmax_pseudo_ellipse(phi, x, y) 413 414 # Repetitive trig. 415 cos_tmax = cos(tmax) 416 sin_tmax2 = sin(tmax)**2 417 cos_phi2 = cos(phi)**2 418 sin_phi2 = sin(phi)**2 419 420 # Components. 421 a = sinc_2smax * (2.0 * sin_tmax2 * sin_phi2 * ((2.0*cos_phi2 - 1.0)*cos_tmax - 6.0*cos_phi2) + 2.0*cos_tmax * (2.0*sin_phi2*(4.0*sin_phi2 - 5.0) + 3.0)) 422 b = 2.0*sin_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 423 424 # Return the theta-sigma integral. 425 return a + b

426 427

428 -def part_int_daeg2_pseudo_ellipse_44(phi, x, y, sinc_2smax):

429 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 430 431 @param phi: The azimuthal tilt-torsion angle. 432 @type phi: float 433 @param x: The cone opening angle along x. 434 @type x: float 435 @param y: The cone opening angle along y. 436 @type y: float 437 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 438 @type sinc_2smax: float 439 @return: The theta-sigma partial integral. 440 @rtype: float 441 """ 442 443 # Theta max. 444 tmax = tmax_pseudo_ellipse(phi, x, y) 445 446 # Repetitive trig. 447 cos_tmax = cos(tmax) 448 sin_tmax2 = sin(tmax)**2 449 cos_phi2 = cos(phi)**2 450 sin_phi2 = sin(phi)**2 451 452 # Components. 453 a = sinc_2smax * (2.0 * sin_tmax2 * sin_phi2 * ((2.0*sin_phi2 - 1.0)*cos_tmax + 6.0*cos_phi2) - 2.0*cos_tmax * (2.0*sin_phi2*(4.0*sin_phi2 - 5.0) + 3.0)) 454 b = 2.0*sin_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 455 456 # Return the theta-sigma integral. 457 return a + b

458 459

460 -def part_int_daeg2_pseudo_ellipse_48(phi, x, y):

461 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 462 463 @param phi: The azimuthal tilt-torsion angle. 464 @type phi: float 465 @param x: The cone opening angle along x. 466 @type x: float 467 @param y: The cone opening angle along y. 468 @type y: float 469 @return: The theta-sigma partial integral. 470 @rtype: float 471 """ 472 473 # Theta max. 474 tmax = tmax_pseudo_ellipse(phi, x, y) 475 476 # The theta-sigma integral. 477 return cos(tmax) * sin(phi)**2 * (sin(tmax)**2 + 2.0)

478 479

480 -def part_int_daeg2_pseudo_ellipse_55(phi, x, y):

481 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 482 483 @param phi: The azimuthal tilt-torsion angle. 484 @type phi: float 485 @param x: The cone opening angle along x. 486 @type x: float 487 @param y: The cone opening angle along y. 488 @type y: float 489 @return: The theta-sigma partial integral. 490 @rtype: float 491 """ 492 493 # Theta max. 494 tmax = tmax_pseudo_ellipse(phi, x, y) 495 496 # Repetitive trig. 497 cos_tmax = cos(tmax) 498 sin_phi2 = sin(phi)**2 499 cos_phi2 = cos(phi)**2 500 501 # Return the theta-sigma integral. 502 return 2.0*sin_phi2*cos_tmax**3 + 3.0*cos_phi2*cos_tmax**2

503 504

505 -def part_int_daeg2_pseudo_ellipse_57(phi, x, y):

506 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 507 508 @param phi: The azimuthal tilt-torsion angle. 509 @type phi: float 510 @param x: The cone opening angle along x. 511 @type x: float 512 @param y: The cone opening angle along y. 513 @type y: float 514 @return: The theta-sigma partial integral. 515 @rtype: float 516 """ 517 518 # Theta max. 519 tmax = tmax_pseudo_ellipse(phi, x, y) 520 521 # Repetitive trig. 522 cos_tmax = cos(tmax) 523 524 # The theta-sigma integral. 525 return sin(phi)**2 * (cos_tmax**3 - 3.0*cos_tmax)

526 527

528 -def part_int_daeg2_pseudo_ellipse_80(phi, x, y, sinc_2smax):

529 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 530 531 @param phi: The azimuthal tilt-torsion angle. 532 @type phi: float 533 @param x: The cone opening angle along x. 534 @type x: float 535 @param y: The cone opening angle along y. 536 @type y: float 537 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 538 @type sinc_2smax: float 539 @return: The theta-sigma partial integral. 540 @rtype: float 541 """ 542 543 # Theta max. 544 tmax = tmax_pseudo_ellipse(phi, x, y) 545 546 # Repetitive trig. 547 cos_tmax = cos(tmax) 548 sin_tmax2 = sin(tmax)**2 549 cos_phi2 = cos(phi)**2 550 551 # The theta-sigma integral. 552 return sinc_2smax * (2.0*(1.0 - 2.0*cos_phi2)*cos_tmax*(sin_tmax2 + 2.0)) + 2.0*cos_tmax**3 - 6.0*cos_tmax

553 554

555 -def part_int_daeg2_pseudo_ellipse_84(phi, x, y, sinc_2smax):

556 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 557 558 @param phi: The azimuthal tilt-torsion angle. 559 @type phi: float 560 @param x: The cone opening angle along x. 561 @type x: float 562 @param y: The cone opening angle along y. 563 @type y: float 564 @param sinc_2smax: The pre-calculated value of sinc(2*sigma_max). 565 @type sinc_2smax: float 566 @return: The theta-sigma partial integral. 567 @rtype: float 568 """ 569 570 # Theta max. 571 tmax = tmax_pseudo_ellipse(phi, x, y) 572 573 # Repetitive trig. 574 cos_tmax = cos(tmax) 575 sin_tmax2 = sin(tmax)**2 576 cos_phi2 = cos(phi)**2 577 578 # The theta-sigma integral. 579 return sinc_2smax * (2.0*(1.0 - 2.0*cos_phi2)*cos_tmax*(sin_tmax2 + 2.0)) - 2.0*cos_tmax**3 + 6.0*cos_tmax

580 581

582 -def part_int_daeg2_pseudo_ellipse_88(phi, x, y):

583 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 584 585 @param phi: The azimuthal tilt-torsion angle. 586 @type phi: float 587 @param x: The cone opening angle along x. 588 @type x: float 589 @param y: The cone opening angle along y. 590 @type y: float 591 @return: The theta-sigma partial integral. 592 @rtype: float 593 """ 594 595 # Theta max. 596 tmax = tmax_pseudo_ellipse(phi, x, y) 597 598 # The theta-sigma integral. 599 return cos(tmax)**3

600 601

602 -def pcs_numeric_qr_int_pseudo_ellipse(points=None, max_points=None, theta_x=None, theta_y=None, sigma_max=None, c=None, full_in_ref_frame=None, r_pivot_atom=None, r_pivot_atom_rev=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None, pcs_theta=None, pcs_theta_err=None, missing_pcs=None):

603 """Determine the averaged PCS value via numerical integration. 604 605 @keyword points: The Sobol points in the torsion-tilt angle space. 606 @type points: numpy rank-2, 3D array 607 @keyword max_points: The maximum number of Sobol' points to use. Once this number is reached, the loop over the Sobol' torsion-tilt angles is terminated. 608 @type max_points: int 609 @keyword theta_x: The x-axis half cone angle. 610 @type theta_x: float 611 @keyword theta_y: The y-axis half cone angle. 612 @type theta_y: float 613 @keyword sigma_max: The maximum torsion angle. 614 @type sigma_max: float 615 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 616 @type c: float 617 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 618 @type full_in_ref_frame: numpy rank-1 array 619 @keyword r_pivot_atom: The pivot point to atom vector. 620 @type r_pivot_atom: numpy rank-1, 3D array 621 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 622 @type r_pivot_atom_rev: numpy rank-2, 3D array 623 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 624 @type r_ln_pivot: numpy rank-1, 3D array 625 @keyword A: The full alignment tensor of the non-moving domain. 626 @type A: numpy rank-2, 3D array 627 @keyword R_eigen: The eigenframe rotation matrix. 628 @type R_eigen: numpy rank-2, 3D array 629 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 630 @type RT_eigen: numpy rank-2, 3D array 631 @keyword Ri_prime: The array of pre-calculated rotation matrices for the in-frame pseudo-elliptic cone motion, used to calculate the PCS for each state i in the numerical integration. 632 @type Ri_prime: numpy rank-3, array of 3D arrays 633 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 634 @type pcs_theta: numpy rank-2 array 635 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 636 @type pcs_theta_err: numpy rank-2 array 637 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 638 @type missing_pcs: numpy rank-2 array 639 """ 640 641 # Clear the data structures. 642 pcs_theta[:] = 0.0 643 pcs_theta_err[:] = 0.0 644 645 # Fast frame shift. 646 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 647 Ri = swapaxes(Ri, 0, 1) 648 649 # Unpack the points. 650 theta, phi, sigma = points 651 652 # Calculate theta_max. 653 theta_max = tmax_pseudo_ellipse_array(phi, theta_x, theta_y) 654 655 # Loop over the samples. 656 num = 0 657 for i in range(len(points[0])): 658 # The maximum number of points has been reached (well, surpassed by one so exit the loop before it is used). 659 if num == max_points: 660 break 661 662 # Check the torsion angle first, for speed. 663 if abs(sigma[i]) > sigma_max: 664 continue 665 666 # As theta_x <= theta_y, check if theta is outside of the isotropic cone defined by theta_y to minimise calculations for speed. 667 if theta[i] > theta_y: 668 continue 669 670 # Outside of the distribution, so skip the point. 671 if theta[i] > theta_max[i]: 672 continue 673 674 # Calculate the PCSs for this state. 675 pcs_pivot_motion_full_qr_int(full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A, Ri=Ri[i], pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) 676 677 # Increment the number of points. 678 num += 1 679 680 # Default to the rigid state if no points lie in the distribution. 681 if num == 0: 682 # Fast identity frame shift. 683 Ri_prime = eye(3, dtype=float64) 684 Ri = dot(R_eigen, tensordot(Ri_prime, RT_eigen, axes=1)) 685 Ri = swapaxes(Ri, 0, 1) 686 687 # Calculate the PCSs for this state. 688 pcs_pivot_motion_full_qr_int(full_in_ref_frame=full_in_ref_frame, r_pivot_atom=r_pivot_atom, r_pivot_atom_rev=r_pivot_atom_rev, r_ln_pivot=r_ln_pivot, A=A, Ri=Ri, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) 689 690 # Multiply the constant. 691 multiply(c, pcs_theta, pcs_theta) 692 693 # Average the PCS. 694 else: 695 multiply(c, pcs_theta, pcs_theta) 696 divide(pcs_theta, float(num), pcs_theta)

697 698

699 -def pcs_numeric_quad_int_pseudo_ellipse(theta_x=None, theta_y=None, sigma_max=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None):

700 """Determine the averaged PCS value via numerical integration. 701 702 @keyword theta_x: The x-axis half cone angle. 703 @type theta_x: float 704 @keyword theta_y: The y-axis half cone angle. 705 @type theta_y: float 706 @keyword sigma_max: The maximum torsion angle. 707 @type sigma_max: float 708 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 709 @type c: float 710 @keyword r_pivot_atom: The pivot point to atom vector. 711 @type r_pivot_atom: numpy rank-1, 3D array 712 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 713 @type r_ln_pivot: numpy rank-1, 3D array 714 @keyword A: The full alignment tensor of the non-moving domain. 715 @type A: numpy rank-2, 3D array 716 @keyword R_eigen: The eigenframe rotation matrix. 717 @type R_eigen: numpy rank-2, 3D array 718 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 719 @type RT_eigen: numpy rank-2, 3D array 720 @keyword Ri_prime: The empty rotation matrix for the in-frame isotropic cone motion, used to calculate the PCS for each state i in the numerical integration. 721 @type Ri_prime: numpy rank-2, 3D array 722 @return: The averaged PCS value. 723 @rtype: float 724 """ 725 726 def pseudo_ellipse(theta, phi): 727 """The pseudo-ellipse wrapper formula.""" 728 729 return tmax_pseudo_ellipse(phi, theta_x, theta_y)

730 731 # Perform numerical integration. 732 result = tplquad(pcs_pivot_motion_full_quad_int, -sigma_max, sigma_max, lambda phi: -pi, lambda phi: pi, lambda theta, phi: 0.0, pseudo_ellipse, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime)) 733 734 # The surface area normalisation factor. 735 SA = 2.0 * sigma_max * pec(theta_x, theta_y) 736 737 # Return the value. 738 return c * result[0] / SA 739 740

741 -def tmax_pseudo_ellipse(phi, theta_x, theta_y):

742 """The pseudo-ellipse tilt-torsion polar angle. 743 744 @param phi: The azimuthal tilt-torsion angle. 745 @type phi: float 746 @param theta_x: The cone opening angle along x. 747 @type theta_x: float 748 @param theta_y: The cone opening angle along y. 749 @type theta_y: float 750 @return: The theta max angle for the given phi angle. 751 @rtype: float 752 """ 753 754 # Zero points. 755 if theta_x == 0.0: 756 return 0.0 757 elif theta_y == 0.0: 758 return 0.0 759 760 # Return the maximum angle. 761 return theta_x * theta_y / sqrt((cos(phi)*theta_y)**2 + (sin(phi)*theta_x)**2)

762 763

764 -def tmax_pseudo_ellipse_array(phi, theta_x, theta_y):

765 """The pseudo-ellipse tilt-torsion polar angle for numpy arrays. 766 767 @param phi: The azimuthal tilt-torsion angle. 768 @type phi: numpy rank-1 float64 array 769 @param theta_x: The cone opening angle along x. 770 @type theta_x: float 771 @param theta_y: The cone opening angle along y. 772 @type theta_y: float 773 @return: The array theta max angles for the given phi angle array. 774 @rtype: numpy rank-1 float64 array 775 """ 776 777 # Zero points. 778 if theta_x == 0.0: 779 return 0.0 * phi 780 elif theta_y == 0.0: 781 return 0.0 * phi 782 783 # Return the maximum angle. 784 return theta_x * theta_y / np_sqrt((np_cos(phi)*theta_y)**2 + (np_sin(phi)*theta_x)**2)

785

Source Code for Module lib.frame_order.pseudo_ellipse