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 sinc 28 try: 29 from scipy.integrate import quad, tplquad 30 except ImportError: 31 pass 32 33 # relax module imports. 34 from lib.geometry.pec import pec 35 from lib.frame_order.matrix_ops import pcs_pivot_motion_full, pcs_pivot_motion_full_qrint, rotate_daeg 36 37

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

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

58 59

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

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

116 117

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

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

138 139

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

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

160 161

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

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

182 183

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

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

213 214

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

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

245 246

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

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

267 268

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

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

297 298

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

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

325 326

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

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

354 355

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

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

376 377

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

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

408 409

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

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

440 441

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

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

462 463

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

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

488 489

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

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

510 511

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

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

537 538

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

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

564 565

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

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

586 587

588 -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):

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

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

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

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

651 652

653 -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):

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

731

Source Code for Module lib.frame_order.pseudo_ellipse