lib.frame_order.pseudo

1 ############################################################################### 2 # # 3 # Copyright (C) 2009-2013 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 # Dependency check module. 26 import dep_check 27 28 # Python module imports. 29 from math import cos, pi, sin, sqrt 30 from numpy import sinc 31 if dep_check.scipy_module: 32 from scipy.integrate import quad, tplquad 33 34 # relax module imports. 35 from lib.geometry.pec import pec 36 from lib.frame_order.matrix_ops import pcs_pivot_motion_full, pcs_pivot_motion_full_qrint, rotate_daeg 37 38

39 -def compile_1st_matrix_pseudo_ellipse(matrix, theta_x, theta_y, sigma_max):

40 """Generate the 1st degree Frame Order matrix for the pseudo-ellipse. 41 42 @param matrix: The Frame Order matrix, 1st degree to be populated. 43 @type matrix: numpy 3D, rank-2 array 44 @param theta_x: The cone opening angle along x. 45 @type theta_x: float 46 @param theta_y: The cone opening angle along y. 47 @type theta_y: float 48 @param sigma_max: The maximum torsion angle. 49 @type sigma_max: float 50 """ 51 52 # The surface area normalisation factor. 53 fact = 1.0 / (2.0 * sigma_max * pec(theta_x, theta_y)) 54 55 # Numerical integration of phi of each element. 56 matrix[0, 0] = fact * quad(part_int_daeg1_pseudo_ellipse_xx, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0] 57 matrix[1, 1] = fact * quad(part_int_daeg1_pseudo_ellipse_yy, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0] 58 matrix[2, 2] = fact * quad(part_int_daeg1_pseudo_ellipse_zz, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]

59 60

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

62 """Generate the 2nd degree Frame Order matrix for the pseudo-ellipse. 63 64 @param matrix: The Frame Order matrix, 2nd degree to be populated. 65 @type matrix: numpy 9D, rank-2 array 66 @param Rx2_eigen: The Kronecker product of the eigenframe rotation matrix with itself. 67 @type Rx2_eigen: numpy 9D, rank-2 array 68 @param theta_x: The cone opening angle along x. 69 @type theta_x: float 70 @param theta_y: The cone opening angle along y. 71 @type theta_y: float 72 @param sigma_max: The maximum torsion angle. 73 @type sigma_max: float 74 """ 75 76 # The surface area normalisation factor. 77 fact = 12.0 * pec(theta_x, theta_y) 78 79 # Invert. 80 if fact == 0.0: 81 fact = 1e100 82 else: 83 fact = 1.0 / fact 84 85 # Sigma_max part. 86 if sigma_max == 0.0: 87 fact2 = 1e100 88 else: 89 fact2 = fact / (2.0 * sigma_max) 90 91 # Diagonal. 92 matrix[0, 0] = fact * (4.0*pi*(sinc(2.0*sigma_max/pi) + 2.0) + quad(part_int_daeg2_pseudo_ellipse_00, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 93 matrix[1, 1] = fact * (4.0*pi*sinc(2.0*sigma_max/pi) + quad(part_int_daeg2_pseudo_ellipse_11, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 94 matrix[2, 2] = fact * 2.0*sinc(sigma_max/pi) * (5.0*pi - quad(part_int_daeg2_pseudo_ellipse_22, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 95 matrix[3, 3] = matrix[1, 1] 96 matrix[4, 4] = fact * (4.0*pi*(sinc(2.0*sigma_max/pi) + 2.0) + quad(part_int_daeg2_pseudo_ellipse_44, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 97 matrix[5, 5] = fact * 2.0*sinc(sigma_max/pi) * (5.0*pi - quad(part_int_daeg2_pseudo_ellipse_55, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 98 matrix[6, 6] = matrix[2, 2] 99 matrix[7, 7] = matrix[5, 5] 100 matrix[8, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_88, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 101 102 # Off diagonal set 1. 103 matrix[0, 4] = fact * (4.0*pi*(2.0 - sinc(2.0*sigma_max/pi)) + quad(part_int_daeg2_pseudo_ellipse_04, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 104 matrix[4, 0] = fact * (4.0*pi*(2.0 - sinc(2.0*sigma_max/pi)) + quad(part_int_daeg2_pseudo_ellipse_40, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 105 matrix[0, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_08, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 106 matrix[8, 0] = fact * (8.0*pi + quad(part_int_daeg2_pseudo_ellipse_80, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 107 matrix[4, 8] = 4.0 * fact * (2.0*pi - quad(part_int_daeg2_pseudo_ellipse_48, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 108 matrix[8, 4] = fact * (8.0*pi - quad(part_int_daeg2_pseudo_ellipse_84, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 109 110 # Off diagonal set 2. 111 matrix[1, 3] = matrix[3, 1] = fact * (4.0*pi*sinc(2.0*sigma_max/pi) + quad(part_int_daeg2_pseudo_ellipse_13, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 112 matrix[2, 6] = matrix[6, 2] = -fact * 4.0 * sinc(sigma_max/pi) * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_26, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 113 matrix[5, 7] = matrix[7, 5] = -fact * 4.0 * sinc(sigma_max/pi) * (2.0*pi + quad(part_int_daeg2_pseudo_ellipse_57, -pi, pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]) 114 115 # Rotate and return the frame order matrix. 116 return rotate_daeg(matrix, Rx2_eigen)

117 118

119 -def part_int_daeg1_pseudo_ellipse_xx(phi, x, y, smax):

120 """The theta-sigma partial integral of the 1st degree Frame Order matrix element xx for the pseudo-ellipse. 121 122 @param phi: The azimuthal tilt-torsion angle. 123 @type phi: float 124 @param x: The cone opening angle along x. 125 @type x: float 126 @param y: The cone opening angle along y. 127 @type y: float 128 @param smax: The maximum torsion angle. 129 @type smax: float 130 @return: The theta-sigma partial integral. 131 @rtype: float 132 """ 133 134 # Theta max. 135 tmax = tmax_pseudo_ellipse(phi, x, y) 136 137 # The theta-sigma integral. 138 return sin(smax) * (2 * (1 - cos(tmax)) * sin(phi)**2 + cos(phi)**2 * sin(tmax)**2)

139 140

141 -def part_int_daeg1_pseudo_ellipse_yy(phi, x, y, smax):

142 """The theta-sigma partial integral of the 1st degree Frame Order matrix element yy for the pseudo-ellipse. 143 144 @param phi: The azimuthal tilt-torsion angle. 145 @type phi: float 146 @param x: The cone opening angle along x. 147 @type x: float 148 @param y: The cone opening angle along y. 149 @type y: float 150 @param smax: The maximum torsion angle. 151 @type smax: 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 sin(smax) * (2.0 * cos(phi)**2 * (1.0 - cos(tmax)) + sin(phi)**2 * sin(tmax)**2)

161 162

163 -def part_int_daeg1_pseudo_ellipse_zz(phi, x, y, smax):

164 """The theta-sigma partial integral of the 1st degree Frame Order matrix element zz 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 @param smax: The maximum torsion angle. 173 @type smax: float 174 @return: The theta-sigma partial integral. 175 @rtype: float 176 """ 177 178 # Theta max. 179 tmax = tmax_pseudo_ellipse(phi, x, y) 180 181 # The theta-sigma integral. 182 return smax * sin(tmax)**2

183 184

185 -def part_int_daeg2_pseudo_ellipse_00(phi, x, y, smax):

186 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 11 for the pseudo-ellipse. 187 188 @param phi: The azimuthal tilt-torsion angle. 189 @type phi: float 190 @param x: The cone opening angle along x. 191 @type x: float 192 @param y: The cone opening angle along y. 193 @type y: float 194 @param smax: The maximum torsion angle. 195 @type smax: float 196 @return: The theta-sigma partial integral. 197 @rtype: float 198 """ 199 200 # Theta max. 201 tmax = tmax_pseudo_ellipse(phi, x, y) 202 203 # Repetitive trig. 204 cos_tmax = cos(tmax) 205 sin_tmax2 = sin(tmax)**2 206 cos_phi2 = cos(phi)**2 207 208 # Components. 209 a = sinc(2.0*smax/pi) * (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)) 210 b = 2.0*cos_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 211 212 # Return the theta-sigma integral. 213 return a + b

214 215

216 -def part_int_daeg2_pseudo_ellipse_04(phi, x, y, smax):

217 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 22 for the pseudo-ellipse. 218 219 @param phi: The azimuthal tilt-torsion angle. 220 @type phi: float 221 @param x: The cone opening angle along x. 222 @type x: float 223 @param y: The cone opening angle along y. 224 @type y: float 225 @param smax: The maximum torsion angle. 226 @type smax: float 227 @return: The theta-sigma partial integral. 228 @rtype: float 229 """ 230 231 # Theta max. 232 tmax = tmax_pseudo_ellipse(phi, x, y) 233 234 # Repetitive trig. 235 cos_tmax = cos(tmax) 236 sin_tmax2 = sin(tmax)**2 237 cos_phi2 = cos(phi)**2 238 sin_phi2 = sin(phi)**2 239 240 # Components. 241 a = sinc(2.0*smax/pi) * (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)) 242 b = 2.0*cos_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 243 244 # Return the theta-sigma integral. 245 return a + b

246 247

248 -def part_int_daeg2_pseudo_ellipse_08(phi, x, y, smax):

249 """The theta-sigma partial integral of the 2nd degree Frame Order matrix element 33 for the pseudo-ellipse. 250 251 @param phi: The azimuthal tilt-torsion angle. 252 @type phi: float 253 @param x: The cone opening angle along x. 254 @type x: float 255 @param y: The cone opening angle along y. 256 @type y: float 257 @param smax: The maximum torsion angle. 258 @type smax: float 259 @return: The theta-sigma partial integral. 260 @rtype: float 261 """ 262 263 # Theta max. 264 tmax = tmax_pseudo_ellipse(phi, x, y) 265 266 # The theta-sigma integral. 267 return cos(tmax) * cos(phi)**2 * (sin(tmax)**2 + 2.0)

268 269

270 -def part_int_daeg2_pseudo_ellipse_11(phi, x, y, smax):

271 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 272 273 @param phi: The azimuthal tilt-torsion angle. 274 @type phi: float 275 @param x: The cone opening angle along x. 276 @type x: float 277 @param y: The cone opening angle along y. 278 @type y: float 279 @param smax: The maximum torsion angle. 280 @type smax: float 281 @return: The theta-sigma partial integral. 282 @rtype: float 283 """ 284 285 # Theta max. 286 tmax = tmax_pseudo_ellipse(phi, x, y) 287 288 # Repetitive trig. 289 cos_tmax = cos(tmax) 290 cos_phi2 = cos(phi)**2 291 sin_phi2 = sin(phi)**2 292 293 # The integral. 294 a = sinc(2.0*smax/pi) * ((4.0*cos_phi2*((1.0 - cos_phi2)*cos_tmax + 3.0*(cos_phi2-1)) + 3.0)*sin(tmax)**2 - 16.0*cos_phi2*sin_phi2*cos_tmax) + 3.0*sin(tmax)**2 295 296 # The theta-sigma integral. 297 return a

298 299

300 -def part_int_daeg2_pseudo_ellipse_13(phi, x, y, smax):

301 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 302 303 @param phi: The azimuthal tilt-torsion angle. 304 @type phi: float 305 @param x: The cone opening angle along x. 306 @type x: float 307 @param y: The cone opening angle along y. 308 @type y: float 309 @param smax: The maximum torsion angle. 310 @type smax: float 311 @return: The theta-sigma partial integral. 312 @rtype: float 313 """ 314 315 # Theta max. 316 tmax = tmax_pseudo_ellipse(phi, x, y) 317 318 # Repetitive trig. 319 cos_tmax = cos(tmax) 320 sin_tmax2 = sin(tmax)**2 321 sinc_2smax = sinc(2.0*smax/pi) 322 cos_sin_phi2 = cos(phi)**2*sin(phi)**2 323 324 # The theta-sigma integral. 325 return sinc_2smax * (sin_tmax2 * (4*cos_sin_phi2*cos_tmax - 12*cos_sin_phi2 + 3) - 16*cos_sin_phi2*cos_tmax) - 3.0*sin_tmax2

326 327

328 -def part_int_daeg2_pseudo_ellipse_22(phi, x, y, smax):

329 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 330 331 @param phi: The azimuthal tilt-torsion angle. 332 @type phi: float 333 @param x: The cone opening angle along x. 334 @type x: float 335 @param y: The cone opening angle along y. 336 @type y: float 337 @param smax: The maximum torsion angle. 338 @type smax: float 339 @return: The theta-sigma partial integral. 340 @rtype: float 341 """ 342 343 # Theta max. 344 tmax = tmax_pseudo_ellipse(phi, x, y) 345 346 # Repetitive trig. 347 cos_tmax = cos(tmax) 348 cos_phi2 = cos(phi)**2 349 350 # Components. 351 a = 2.0*cos_phi2*cos_tmax**3 + 3.0*(1.0 - cos_phi2)*cos_tmax**2 352 353 # Return the theta-sigma integral. 354 return a

355 356

357 -def part_int_daeg2_pseudo_ellipse_26(phi, x, y, smax):

358 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 359 360 @param phi: The azimuthal tilt-torsion angle. 361 @type phi: float 362 @param x: The cone opening angle along x. 363 @type x: float 364 @param y: The cone opening angle along y. 365 @type y: float 366 @param smax: The maximum torsion angle. 367 @type smax: float 368 @return: The theta-sigma partial integral. 369 @rtype: float 370 """ 371 372 # Theta max. 373 tmax = tmax_pseudo_ellipse(phi, x, y) 374 375 # The theta-sigma integral. 376 return cos(phi)**2 * (cos(tmax)**3 - 3.0*cos(tmax))

377 378

379 -def part_int_daeg2_pseudo_ellipse_40(phi, x, y, smax):

380 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 381 382 @param phi: The azimuthal tilt-torsion angle. 383 @type phi: float 384 @param x: The cone opening angle along x. 385 @type x: float 386 @param y: The cone opening angle along y. 387 @type y: float 388 @param smax: The maximum torsion angle. 389 @type smax: float 390 @return: The theta-sigma partial integral. 391 @rtype: float 392 """ 393 394 # Theta max. 395 tmax = tmax_pseudo_ellipse(phi, x, y) 396 397 # Repetitive trig. 398 cos_tmax = cos(tmax) 399 sin_tmax2 = sin(tmax)**2 400 cos_phi2 = cos(phi)**2 401 sin_phi2 = sin(phi)**2 402 403 # Components. 404 a = sinc(2.0*smax/pi) * (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)) 405 b = 2.0*sin_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 406 407 # Return the theta-sigma integral. 408 return a + b

409 410

411 -def part_int_daeg2_pseudo_ellipse_44(phi, x, y, smax):

412 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 413 414 @param phi: The azimuthal tilt-torsion angle. 415 @type phi: float 416 @param x: The cone opening angle along x. 417 @type x: float 418 @param y: The cone opening angle along y. 419 @type y: float 420 @param smax: The maximum torsion angle. 421 @type smax: float 422 @return: The theta-sigma partial integral. 423 @rtype: float 424 """ 425 426 # Theta max. 427 tmax = tmax_pseudo_ellipse(phi, x, y) 428 429 # Repetitive trig. 430 cos_tmax = cos(tmax) 431 sin_tmax2 = sin(tmax)**2 432 cos_phi2 = cos(phi)**2 433 sin_phi2 = sin(phi)**2 434 435 # Components. 436 a = sinc(2.0*smax/pi) * (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)) 437 b = 2.0*sin_phi2*cos_tmax*(sin_tmax2 + 2.0) - 6.0*cos_tmax 438 439 # Return the theta-sigma integral. 440 return a + b

441 442

443 -def part_int_daeg2_pseudo_ellipse_48(phi, x, y, smax):

444 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 445 446 @param phi: The azimuthal tilt-torsion angle. 447 @type phi: float 448 @param x: The cone opening angle along x. 449 @type x: float 450 @param y: The cone opening angle along y. 451 @type y: float 452 @param smax: The maximum torsion angle. 453 @type smax: float 454 @return: The theta-sigma partial integral. 455 @rtype: float 456 """ 457 458 # Theta max. 459 tmax = tmax_pseudo_ellipse(phi, x, y) 460 461 # The theta-sigma integral. 462 return cos(tmax) * sin(phi)**2 * (sin(tmax)**2 + 2.0)

463 464

465 -def part_int_daeg2_pseudo_ellipse_55(phi, x, y, smax):

466 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 467 468 @param phi: The azimuthal tilt-torsion angle. 469 @type phi: float 470 @param x: The cone opening angle along x. 471 @type x: float 472 @param y: The cone opening angle along y. 473 @type y: float 474 @param smax: The maximum torsion angle. 475 @type smax: float 476 @return: The theta-sigma partial integral. 477 @rtype: float 478 """ 479 480 # Theta max. 481 tmax = tmax_pseudo_ellipse(phi, x, y) 482 483 # Repetitive trig. 484 cos_tmax = cos(tmax) 485 sin_phi2 = sin(phi)**2 486 487 # Return the theta-sigma integral. 488 return 2.0*sin_phi2*cos_tmax**3 + 3.0*(1.0 - sin_phi2)*cos_tmax**2

489 490

491 -def part_int_daeg2_pseudo_ellipse_57(phi, x, y, smax):

492 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 493 494 @param phi: The azimuthal tilt-torsion angle. 495 @type phi: float 496 @param x: The cone opening angle along x. 497 @type x: float 498 @param y: The cone opening angle along y. 499 @type y: float 500 @param smax: The maximum torsion angle. 501 @type smax: float 502 @return: The theta-sigma partial integral. 503 @rtype: float 504 """ 505 506 # Theta max. 507 tmax = tmax_pseudo_ellipse(phi, x, y) 508 509 # The theta-sigma integral. 510 return sin(phi)**2 * (cos(tmax)**3 - 3.0*cos(tmax))

511 512

513 -def part_int_daeg2_pseudo_ellipse_80(phi, x, y, smax):

514 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 515 516 @param phi: The azimuthal tilt-torsion angle. 517 @type phi: float 518 @param x: The cone opening angle along x. 519 @type x: float 520 @param y: The cone opening angle along y. 521 @type y: float 522 @param smax: The maximum torsion angle. 523 @type smax: float 524 @return: The theta-sigma partial integral. 525 @rtype: float 526 """ 527 528 # Theta max. 529 tmax = tmax_pseudo_ellipse(phi, x, y) 530 531 # Repetitive trig. 532 cos_tmax = cos(tmax) 533 sin_tmax2 = sin(tmax)**2 534 cos_phi2 = cos(phi)**2 535 536 # The theta-sigma integral. 537 return sinc(2.0*smax/pi) * (2.0*(1.0 - 2.0*cos_phi2)*cos_tmax*(sin_tmax2 + 2.0)) + 2.0*cos_tmax**3 - 6.0*cos_tmax

538 539

540 -def part_int_daeg2_pseudo_ellipse_84(phi, x, y, smax):

541 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 542 543 @param phi: The azimuthal tilt-torsion angle. 544 @type phi: float 545 @param x: The cone opening angle along x. 546 @type x: float 547 @param y: The cone opening angle along y. 548 @type y: float 549 @param smax: The maximum torsion angle. 550 @type smax: float 551 @return: The theta-sigma partial integral. 552 @rtype: float 553 """ 554 555 # Theta max. 556 tmax = tmax_pseudo_ellipse(phi, x, y) 557 558 # Repetitive trig. 559 cos_tmax = cos(tmax) 560 sin_tmax2 = sin(tmax)**2 561 cos_phi2 = cos(phi)**2 562 563 # The theta-sigma integral. 564 return sinc(2.0*smax/pi) * (2.0*(1.0 - 2.0*cos_phi2)*cos_tmax*(sin_tmax2 + 2.0)) - 2.0*cos_tmax**3 + 6.0*cos_tmax

565 566

567 -def part_int_daeg2_pseudo_ellipse_88(phi, x, y, smax):

568 """The theta-sigma partial integral of the 2nd degree Frame Order matrix for the pseudo-ellipse. 569 570 @param phi: The azimuthal tilt-torsion angle. 571 @type phi: float 572 @param x: The cone opening angle along x. 573 @type x: float 574 @param y: The cone opening angle along y. 575 @type y: float 576 @param smax: The maximum torsion angle. 577 @type smax: float 578 @return: The theta-sigma partial integral. 579 @rtype: float 580 """ 581 582 # Theta max. 583 tmax = tmax_pseudo_ellipse(phi, x, y) 584 585 # The theta-sigma integral. 586 return cos(tmax)**3

587 588

589 -def pcs_numeric_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):

590 """Determine the averaged PCS value via numerical integration. 591 592 @keyword theta_x: The x-axis half cone angle. 593 @type theta_x: float 594 @keyword theta_y: The y-axis half cone angle. 595 @type theta_y: float 596 @keyword sigma_max: The maximum torsion angle. 597 @type sigma_max: float 598 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 599 @type c: float 600 @keyword r_pivot_atom: The pivot point to atom vector. 601 @type r_pivot_atom: numpy rank-1, 3D array 602 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 603 @type r_ln_pivot: numpy rank-1, 3D array 604 @keyword A: The full alignment tensor of the non-moving domain. 605 @type A: numpy rank-2, 3D array 606 @keyword R_eigen: The eigenframe rotation matrix. 607 @type R_eigen: numpy rank-2, 3D array 608 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 609 @type RT_eigen: numpy rank-2, 3D array 610 @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. 611 @type Ri_prime: numpy rank-2, 3D array 612 @return: The averaged PCS value. 613 @rtype: float 614 """ 615 616 def pseudo_ellipse(theta, phi): 617 """The pseudo-ellipse wrapper formula.""" 618 619 return tmax_pseudo_ellipse(phi, theta_x, theta_y)

620 621 # Perform numerical integration. 622 result = tplquad(pcs_pivot_motion_full, -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)) 623 624 # The surface area normalisation factor. 625 SA = 2.0 * sigma_max * pec(theta_x, theta_y) 626 627 # Return the value. 628 return c * result[0] / SA 629 630

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

632 """The pseudo-ellipse tilt-torsion polar angle. 633 634 @param phi: The azimuthal tilt-torsion angle. 635 @type phi: float 636 @param theta_x: The cone opening angle along x. 637 @type theta_x: float 638 @param theta_y: The cone opening angle along y. 639 @type theta_y: float 640 @return: The theta max angle for the given phi angle. 641 @rtype: float 642 """ 643 644 # Zero points. 645 if theta_x == 0.0: 646 return 0.0 647 elif theta_y == 0.0: 648 return 0.0 649 650 # Return the maximum angle. 651 return theta_x * theta_y / sqrt((cos(phi)*theta_y)**2 + (sin(phi)*theta_x)**2)

652 653

654 -def pcs_numeric_int_pseudo_ellipse_qrint(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, error_flag=False):

655 """Determine the averaged PCS value via numerical integration. 656 657 @keyword points: The Sobol points in the torsion-tilt angle space. 658 @type points: numpy rank-2, 3D array 659 @keyword theta_x: The x-axis half cone angle. 660 @type theta_x: float 661 @keyword theta_y: The y-axis half cone angle. 662 @type theta_y: float 663 @keyword sigma_max: The maximum torsion angle. 664 @type sigma_max: float 665 @keyword c: The PCS constant (without the interatomic distance and in Angstrom units). 666 @type c: float 667 @keyword full_in_ref_frame: An array of flags specifying if the tensor in the reference frame is the full or reduced tensor. 668 @type full_in_ref_frame: numpy rank-1 array 669 @keyword r_pivot_atom: The pivot point to atom vector. 670 @type r_pivot_atom: numpy rank-1, 3D array 671 @keyword r_pivot_atom_rev: The reversed pivot point to atom vector. 672 @type r_pivot_atom_rev: numpy rank-2, 3D array 673 @keyword r_ln_pivot: The lanthanide position to pivot point vector. 674 @type r_ln_pivot: numpy rank-1, 3D array 675 @keyword A: The full alignment tensor of the non-moving domain. 676 @type A: numpy rank-2, 3D array 677 @keyword R_eigen: The eigenframe rotation matrix. 678 @type R_eigen: numpy rank-2, 3D array 679 @keyword RT_eigen: The transpose of the eigenframe rotation matrix (for faster calculations). 680 @type RT_eigen: numpy rank-2, 3D array 681 @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. 682 @type Ri_prime: numpy rank-2, 3D array 683 @keyword pcs_theta: The storage structure for the back-calculated PCS values. 684 @type pcs_theta: numpy rank-2 array 685 @keyword pcs_theta_err: The storage structure for the back-calculated PCS errors. 686 @type pcs_theta_err: numpy rank-2 array 687 @keyword missing_pcs: A structure used to indicate which PCS values are missing. 688 @type missing_pcs: numpy rank-2 array 689 @keyword error_flag: A flag which if True will cause the PCS errors to be estimated and stored in pcs_theta_err. 690 @type error_flag: bool 691 """ 692 693 694 # Clear the data structures. 695 for i in range(len(pcs_theta)): 696 for j in range(len(pcs_theta[i])): 697 pcs_theta[i, j] = 0.0 698 pcs_theta_err[i, j] = 0.0 699 700 # Loop over the samples. 701 num = 0 702 for i in range(len(points)): 703 # Unpack the point. 704 theta, phi, sigma = points[i] 705 706 # Calculate theta_max. 707 theta_max = tmax_pseudo_ellipse(phi, theta_x, theta_y) 708 709 # Outside of the distribution, so skip the point. 710 if theta > theta_max: 711 continue 712 if sigma > sigma_max or sigma < -sigma_max: 713 continue 714 715 # Calculate the PCSs for this state. 716 pcs_pivot_motion_full_qrint(theta_i=theta, phi_i=phi, sigma_i=sigma, 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, R_eigen=R_eigen, RT_eigen=RT_eigen, Ri_prime=Ri_prime, pcs_theta=pcs_theta, pcs_theta_err=pcs_theta_err, missing_pcs=missing_pcs) 717 718 # Increment the number of points. 719 num += 1 720 721 # Calculate the PCS and error. 722 for i in range(len(pcs_theta)): 723 for j in range(len(pcs_theta[i])): 724 # The average PCS. 725 pcs_theta[i, j] = c[i] * pcs_theta[i, j] / float(num) 726 727 # The error. 728 if error_flag: 729 pcs_theta_err[i, j] = abs(pcs_theta_err[i, j] / float(num) - pcs_theta[i, j]**2) / float(num) 730 pcs_theta_err[i, j] = c[i] * sqrt(pcs_theta_err[i, j]) 731 print("%8.3f +/- %-8.3f" % (pcs_theta[i, j]*1e6, pcs_theta_err[i, j]*1e6))

732

Source Code for Module lib.frame_order.pseudo_ellipse