generic_fns.diffusion

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2012 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax. # 6 # # 7 # relax 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 2 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # relax 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 relax; if not, write to the Free Software # 19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for the support of diffusion tensors.""" 25 26 # Python module imports. 27 from copy import deepcopy 28 from math import cos, pi, sin 29 from numpy import cross, float64, int32, ones, transpose, zeros 30 from numpy.linalg import norm, svd 31 from operator import itemgetter 32 from re import search 33 import string 34 35 # relax module imports. 36 from angles import wrap_angles 37 from data.diff_tensor import DiffTensorData 38 from generic_fns import pipes 39 from generic_fns.angles import fold_spherical_angles 40 from generic_fns.mol_res_spin import get_molecule_names, spin_loop 41 from maths_fns.coord_transform import cartesian_to_spherical 42 from maths_fns.rotation_matrix import R_to_euler_zyz 43 from physical_constants import element_from_isotope, number_from_isotope 44 from relax_errors import RelaxError, RelaxNoTensorError, RelaxStrError, RelaxTensorError, RelaxUnknownParamCombError, RelaxUnknownParamError 45 46

47 -def bmrb_read(star):

48 """Read the relaxation data from the NMR-STAR dictionary object. 49 50 @param star: The NMR-STAR dictionary object. 51 @type star: NMR_STAR instance 52 """ 53 54 # Get the diffusion tensor data. 55 found = 0 56 for data in star.tensor.loop(): 57 # No data. 58 if data == None: 59 continue 60 61 # Not a diffusion tensor. 62 if data['tensor_type'] != 'diffusion': 63 continue 64 65 # Found. 66 found = found + 1 67 68 # No diffusion tensor data. 69 if not found: 70 return 71 72 # Check. 73 if found != 1: 74 raise RelaxError("More than one diffusion tensor found.") 75 76 # Rebuild the tensor. 77 tensor = zeros((3, 3), float64) 78 tensor[0, 0] = data['tensor_11'][0] 79 tensor[0, 1] = data['tensor_12'][0] 80 tensor[0, 2] = data['tensor_13'][0] 81 tensor[1, 0] = data['tensor_21'][0] 82 tensor[1, 1] = data['tensor_22'][0] 83 tensor[1, 2] = data['tensor_23'][0] 84 tensor[2, 0] = data['tensor_31'][0] 85 tensor[2, 1] = data['tensor_32'][0] 86 tensor[2, 2] = data['tensor_33'][0] 87 88 # Eigenvalues. 89 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 90 91 # X-Y eigenvalue comparison. 92 xy_match = False 93 epsilon = 1e-1 94 if abs(Di[0] - Di[1]) < epsilon: 95 xy_match = True 96 97 # Y-Z eigenvalue comparison. 98 yz_match = False 99 if abs(Di[1] - Di[2]) < epsilon: 100 yz_match = True 101 102 # Determine the tensor type. 103 if xy_match and yz_match: 104 shape = ['sphere'] 105 elif xy_match: 106 shape = ['spheroid', 'prolate spheroid'] 107 type = 'prolate' 108 Dpar = Di[2] 109 Dper = Di[0] 110 elif yz_match: 111 shape = ['spheroid', 'oblate spheroid'] 112 type = 'oblate' 113 Dpar = Di[0] 114 Dper = Di[2] 115 else: 116 shape = ['ellipsoid'] 117 118 # Check the shape. 119 if data['geometric_shape'] not in shape: 120 raise RelaxError("The tensor with eigenvalues %s does not match the %s geometric shape." % (Di, shape[0])) 121 122 # Add the diff_tensor object to the data pipe. 123 cdp.diff_tensor = DiffTensorData() 124 125 # Set the fixed flag. 126 cdp.diff_tensor.fixed = True 127 128 # Sphere. 129 if data['geometric_shape'] == 'sphere': 130 sphere(params=Di[0], d_scale=1.0, param_types=1) 131 132 # Spheroid. 133 elif data['geometric_shape'] in ['spheroid', 'oblate spheroid', 'prolate spheroid']: 134 # The spherical angles. 135 theta = data['axial_sym_axis_polar_angle'][0] 136 phi = data['axial_sym_axis_azimuthal_angle'][0] 137 138 # Set up the tensor. 139 spheroid(params=(Dpar, Dper, theta, phi), d_scale=1.0, param_types=3, spheroid_type=type) 140 141 # Ellipsoid. 142 elif data['geometric_shape'] == 'ellipsoid': 143 ellipsoid(params=(Di[0], Di[1], Di[2], alpha, beta, gamma), d_scale=1.0, param_types=3)

144 145

146 -def bmrb_write(star):

147 """Generate the diffusion tensor saveframes for the NMR-STAR dictionary object. 148 149 @param star: The NMR-STAR dictionary object. 150 @type star: NMR_STAR instance 151 """ 152 153 # Get the current data pipe. 154 cdp = pipes.get_pipe() 155 156 # Initialise the spin specific data lists. 157 mol_name_list = [] 158 res_num_list = [] 159 res_name_list = [] 160 atom_name_list = [] 161 isotope_list = [] 162 element_list = [] 163 attached_atom_name_list = [] 164 attached_isotope_list = [] 165 attached_element_list = [] 166 167 # Store the spin specific data in lists for later use. 168 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True): 169 # Skip deselected spins. 170 if not spin.select: 171 continue 172 173 # Check the data for None (not allowed in BMRB!). 174 if res_num == None: 175 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id) 176 if res_name == None: 177 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id) 178 if spin.name == None: 179 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id) 180 if spin.heteronuc_type == None: 181 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id) 182 183 # The molecule/residue/spin info. 184 mol_name_list.append(mol_name) 185 res_num_list.append(str(res_num)) 186 res_name_list.append(str(res_name)) 187 atom_name_list.append(str(spin.name)) 188 189 # The attached atom info. 190 if hasattr(spin, 'attached_atom'): 191 attached_atom_name_list.append(str(spin.attached_atom)) 192 else: 193 attached_atom_name_list.append(str(spin.attached_proton)) 194 attached_element_list.append(element_from_isotope(spin.proton_type)) 195 attached_isotope_list.append(str(number_from_isotope(spin.proton_type))) 196 197 # Other info. 198 isotope_list.append(int(string.strip(spin.heteronuc_type, string.ascii_letters))) 199 element_list.append(spin.element) 200 201 # Convert the molecule names into the entity IDs. 202 entity_ids = zeros(len(mol_name_list), int32) 203 mol_names = get_molecule_names() 204 for i in range(len(mol_name_list)): 205 for j in range(len(mol_names)): 206 if mol_name_list[i] == mol_names[j]: 207 entity_ids[i] = j+1 208 209 # The tensor geometric shape. 210 geometric_shape = cdp.diff_tensor.type 211 if geometric_shape == 'spheroid': 212 geometric_shape = "%s %s" % (cdp.diff_tensor.spheroid_type, geometric_shape) 213 214 # The tensor symmetry. 215 shapes = ['sphere', 'oblate spheroid', 'prolate spheroid', 'ellipsoid'] 216 sym = ['isotropic', 'axial symmetry', 'axial symmetry', 'rhombic'] 217 for i in range(len(shapes)): 218 if geometric_shape == shapes[i]: 219 tensor_symmetry = sym[i] 220 221 # Axial symmetry axis. 222 theta = None 223 phi = None 224 if tensor_symmetry == 'axial symmetry': 225 theta = cdp.diff_tensor.theta 226 phi = cdp.diff_tensor.phi 227 228 # Euler angles. 229 alpha, beta, gamma = None, None, None 230 if tensor_symmetry == 'rhombic': 231 alpha = cdp.diff_tensor.alpha 232 beta = cdp.diff_tensor.beta 233 gamma = cdp.diff_tensor.gamma 234 235 # The tensor eigenvalues. 236 Diso = cdp.diff_tensor.Diso 237 Da = None 238 Dr = None 239 if tensor_symmetry == 'axial symmetry': 240 Da = cdp.diff_tensor.Da 241 elif tensor_symmetry == 'rhombic': 242 Dr = cdp.diff_tensor.Dr 243 244 # The full tensor. 245 tensor_11 = cdp.diff_tensor.tensor[0, 0] 246 tensor_12 = cdp.diff_tensor.tensor[0, 1] 247 tensor_13 = cdp.diff_tensor.tensor[0, 2] 248 tensor_21 = cdp.diff_tensor.tensor[1, 0] 249 tensor_22 = cdp.diff_tensor.tensor[1, 1] 250 tensor_23 = cdp.diff_tensor.tensor[1, 2] 251 tensor_31 = cdp.diff_tensor.tensor[2, 0] 252 tensor_32 = cdp.diff_tensor.tensor[2, 1] 253 tensor_33 = cdp.diff_tensor.tensor[2, 2] 254 255 256 # Add the diffusion tensor. 257 star.tensor.add(tensor_type='diffusion', euler_type='zyz', geometric_shape=geometric_shape, tensor_symmetry=tensor_symmetry, matrix_val_units='s-1', angle_units='rad', iso_val_formula='Diso = 1/(6.tm)', aniso_val_formula='Da = Dpar - Dper', rhomb_val_formula='Dr = (Dy - Dx)/2Da', entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, axial_sym_axis_polar_angle=theta, axial_sym_axis_azimuthal_angle=phi, iso_val=Diso, aniso_val=Da, rhombic_val=Dr, euler_alpha=alpha, euler_beta=beta, euler_gamma=gamma, tensor_11=tensor_11, tensor_12=tensor_12, tensor_13=tensor_13, tensor_21=tensor_21, tensor_22=tensor_22, tensor_23=tensor_23, tensor_31=tensor_31, tensor_32=tensor_32, tensor_33=tensor_33)

258 259 260

261 -def copy(pipe_from=None, pipe_to=None):

262 """Function for copying diffusion tensor data from one data pipe to another. 263 264 @param pipe_from: The data pipe to copy the diffusion tensor data from. This defaults to the 265 current data pipe. 266 @type pipe_from: str 267 @param pipe_to: The data pipe to copy the diffusion tensor data to. This defaults to the 268 current data pipe. 269 @type pipe_to: str 270 """ 271 272 # Defaults. 273 if pipe_from == None and pipe_to == None: 274 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 275 elif pipe_from == None: 276 pipe_from = pipes.cdp_name() 277 elif pipe_to == None: 278 pipe_to = pipes.cdp_name() 279 280 # Test if the pipe_from and pipe_to data pipes exist. 281 pipes.test(pipe_from) 282 pipes.test(pipe_to) 283 284 # Get the data pipes. 285 dp_from = pipes.get_pipe(pipe_from) 286 dp_to = pipes.get_pipe(pipe_to) 287 288 # Test if pipe_from contains diffusion tensor data. 289 if not diff_data_exists(pipe_from): 290 raise RelaxNoTensorError('diffusion') 291 292 # Test if pipe_to contains diffusion tensor data. 293 if diff_data_exists(pipe_to): 294 raise RelaxTensorError('diffusion') 295 296 # Copy the data. 297 dp_to.diff_tensor = deepcopy(dp_from.diff_tensor)

298 299

300 -def data_names():

301 """Function for returning a list of names of data structures associated with the sequence.""" 302 303 names = [ 'diff_type', 304 'diff_params' ] 305 306 return names

307 308

309 -def default_value(param):

310 """Return the default values for the diffusion tensor parameters. 311 312 @param param: The name of the parameter. 313 @type param: str 314 @return: The default value. 315 @rtype: float 316 """ 317 318 # tm. 319 if param == 'tm': 320 return 10.0 * 1e-9 321 322 # Diso, Dx, Dy, Dz, Dpar, Dper. 323 elif param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 324 return 1.666 * 1e7 325 326 # Da, Dr. 327 elif param == 'Da' or param == 'Dr': 328 return 0.0 329 330 # Dratio. 331 elif param == 'Dratio': 332 return 1.0 333 334 # All angles. 335 elif param == 'alpha' or param == 'beta' or param == 'gamma' or param == 'theta' or param == 'phi': 336 return 0.0

337 338 # User function documentation. 339 __default_value_prompt_doc__ = ["Diffusion tensor parameter default values", """ 340 ________________________________________________________________________ 341 | | | | 342 | Data type | Object name | Value | 343 |________________________|____________________|________________________| 344 | | | | 345 | tm | 'tm' | 10 * 1e-9 | 346 | | | | 347 | Diso | 'Diso' | 1.666 * 1e7 | 348 | | | | 349 | Da | 'Da' | 0.0 | 350 | | | | 351 | Dr | 'Dr' | 0.0 | 352 | | | | 353 | Dx | 'Dx' | 1.666 * 1e7 | 354 | | | | 355 | Dy | 'Dy' | 1.666 * 1e7 | 356 | | | | 357 | Dz | 'Dz' | 1.666 * 1e7 | 358 | | | | 359 | Dpar | 'Dpar' | 1.666 * 1e7 | 360 | | | | 361 | Dper | 'Dper' | 1.666 * 1e7 | 362 | | | | 363 | Dratio | 'Dratio' | 1.0 | 364 | | | | 365 | alpha | 'alpha' | 0.0 | 366 | | | | 367 | beta | 'beta' | 0.0 | 368 | | | | 369 | gamma | 'gamma' | 0.0 | 370 | | | | 371 | theta | 'theta' | 0.0 | 372 | | | | 373 | phi | 'phi' | 0.0 | 374 |________________________|____________________|________________________| 375 """] 376 377

378 -def delete():

379 """Function for deleting diffusion tensor data.""" 380 381 # Test if the current data pipe exists. 382 pipes.test() 383 384 # Test if diffusion tensor data exists. 385 if not diff_data_exists(): 386 raise RelaxNoTensorError('diffusion') 387 388 # Delete the diffusion data. 389 del(cdp.diff_tensor)

390 391

392 -def diff_data_exists(pipe=None):

393 """Function for determining if diffusion data exists in the current data pipe. 394 395 @param pipe: The data pipe to search for data in. 396 @type pipe: str 397 @return: The answer to the question. 398 @rtype: bool 399 """ 400 401 # The data pipe to check. 402 if pipe == None: 403 pipe = pipes.cdp_name() 404 405 # Get the data pipe. 406 dp = pipes.get_pipe(pipe) 407 408 # Test if the data structure exists. 409 if hasattr(dp, 'diff_tensor'): 410 return True 411 else: 412 return False

413 414

415 -def display():

416 """Function for displaying the diffusion tensor.""" 417 418 # Test if the current data pipe exists. 419 pipes.test() 420 421 # Test if diffusion tensor data exists. 422 if not diff_data_exists(): 423 raise RelaxNoTensorError('diffusion') 424 425 # Spherical diffusion. 426 if cdp.diff_tensor.type == 'sphere': 427 # Tensor type. 428 print("Type: Spherical diffusion") 429 430 # Parameters. 431 print("\nParameters {tm}.") 432 print(("tm (s): " + repr(cdp.diff_tensor.tm))) 433 434 # Alternate parameters. 435 print("\nAlternate parameters {Diso}.") 436 print(("Diso (1/s): " + repr(cdp.diff_tensor.Diso))) 437 438 # Fixed flag. 439 print(("\nFixed: " + repr(cdp.diff_tensor.fixed))) 440 441 # Spheroidal diffusion. 442 elif cdp.diff_tensor.type == 'spheroid': 443 # Tensor type. 444 print("Type: Spheroidal diffusion") 445 446 # Parameters. 447 print("\nParameters {tm, Da, theta, phi}.") 448 print(("tm (s): " + repr(cdp.diff_tensor.tm))) 449 print(("Da (1/s): " + repr(cdp.diff_tensor.Da))) 450 print(("theta (rad): " + repr(cdp.diff_tensor.theta))) 451 print(("phi (rad): " + repr(cdp.diff_tensor.phi))) 452 453 # Alternate parameters. 454 print("\nAlternate parameters {Diso, Da, theta, phi}.") 455 print(("Diso (1/s): " + repr(cdp.diff_tensor.Diso))) 456 print(("Da (1/s): " + repr(cdp.diff_tensor.Da))) 457 print(("theta (rad): " + repr(cdp.diff_tensor.theta))) 458 print(("phi (rad): " + repr(cdp.diff_tensor.phi))) 459 460 # Alternate parameters. 461 print("\nAlternate parameters {Dpar, Dper, theta, phi}.") 462 print(("Dpar (1/s): " + repr(cdp.diff_tensor.Dpar))) 463 print(("Dper (1/s): " + repr(cdp.diff_tensor.Dper))) 464 print(("theta (rad): " + repr(cdp.diff_tensor.theta))) 465 print(("phi (rad): " + repr(cdp.diff_tensor.phi))) 466 467 # Alternate parameters. 468 print("\nAlternate parameters {tm, Dratio, theta, phi}.") 469 print(("tm (s): " + repr(cdp.diff_tensor.tm))) 470 print(("Dratio: " + repr(cdp.diff_tensor.Dratio))) 471 print(("theta (rad): " + repr(cdp.diff_tensor.theta))) 472 print(("phi (rad): " + repr(cdp.diff_tensor.phi))) 473 474 # Fixed flag. 475 print(("\nFixed: " + repr(cdp.diff_tensor.fixed))) 476 477 # Ellipsoidal diffusion. 478 elif cdp.diff_tensor.type == 'ellipsoid': 479 # Tensor type. 480 print("Type: Ellipsoidal diffusion") 481 482 # Parameters. 483 print("\nParameters {tm, Da, Dr, alpha, beta, gamma}.") 484 print(("tm (s): " + repr(cdp.diff_tensor.tm))) 485 print(("Da (1/s): " + repr(cdp.diff_tensor.Da))) 486 print(("Dr: " + repr(cdp.diff_tensor.Dr))) 487 print(("alpha (rad): " + repr(cdp.diff_tensor.alpha))) 488 print(("beta (rad): " + repr(cdp.diff_tensor.beta))) 489 print(("gamma (rad): " + repr(cdp.diff_tensor.gamma))) 490 491 # Alternate parameters. 492 print("\nAlternate parameters {Diso, Da, Dr, alpha, beta, gamma}.") 493 print(("Diso (1/s): " + repr(cdp.diff_tensor.Diso))) 494 print(("Da (1/s): " + repr(cdp.diff_tensor.Da))) 495 print(("Dr: " + repr(cdp.diff_tensor.Dr))) 496 print(("alpha (rad): " + repr(cdp.diff_tensor.alpha))) 497 print(("beta (rad): " + repr(cdp.diff_tensor.beta))) 498 print(("gamma (rad): " + repr(cdp.diff_tensor.gamma))) 499 500 # Alternate parameters. 501 print("\nAlternate parameters {Dx, Dy, Dz, alpha, beta, gamma}.") 502 print(("Dx (1/s): " + repr(cdp.diff_tensor.Dx))) 503 print(("Dy (1/s): " + repr(cdp.diff_tensor.Dy))) 504 print(("Dz (1/s): " + repr(cdp.diff_tensor.Dz))) 505 print(("alpha (rad): " + repr(cdp.diff_tensor.alpha))) 506 print(("beta (rad): " + repr(cdp.diff_tensor.beta))) 507 print(("gamma (rad): " + repr(cdp.diff_tensor.gamma))) 508 509 # Fixed flag. 510 print(("\nFixed: " + repr(cdp.diff_tensor.fixed)))

511 512

513 -def ellipsoid(params=None, time_scale=None, d_scale=None, angle_units=None, param_types=None):

514 """Function for setting up a ellipsoidal diffusion tensor. 515 516 @param params: The diffusion tensor parameter. 517 @type params: float 518 @param time_scale: The correlation time scaling value. 519 @type time_scale: float 520 @param d_scale: The diffusion tensor eigenvalue scaling value. 521 @type d_scale: float 522 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 523 @type angle_units: str 524 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 525 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 526 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 527 @type param_types: int 528 """ 529 530 # The diffusion type. 531 cdp.diff_tensor.type = 'ellipsoid' 532 533 # (tm, Da, Dr, alpha, beta, gamma). 534 if param_types == 0: 535 # Unpack the tuple. 536 tm, Da, Dr, alpha, beta, gamma = params 537 538 # Scaling. 539 tm = tm * time_scale 540 Da = Da * d_scale 541 542 # Set the parameters. 543 set(value=[tm, Da, Dr], param=['tm', 'Da', 'Dr']) 544 545 # (Diso, Da, Dr, alpha, beta, gamma). 546 elif param_types == 1: 547 # Unpack the tuple. 548 Diso, Da, Dr, alpha, beta, gamma = params 549 550 # Scaling. 551 Diso = Diso * d_scale 552 Da = Da * d_scale 553 554 # Set the parameters. 555 set(value=[Diso, Da, Dr], param=['Diso', 'Da', 'Dr']) 556 557 # (Dx, Dy, Dz, alpha, beta, gamma). 558 elif param_types == 2: 559 # Unpack the tuple. 560 Dx, Dy, Dz, alpha, beta, gamma = params 561 562 # Scaling. 563 Dx = Dx * d_scale 564 Dy = Dy * d_scale 565 Dz = Dz * d_scale 566 567 # Set the parameters. 568 set(value=[Dx, Dy, Dz], param=['Dx', 'Dy', 'Dz']) 569 570 # (Dxx, Dyy, Dzz, Dxy, Dxz, Dyz). 571 elif param_types == 3: 572 # Unpack the tuple. 573 Dxx, Dyy, Dzz, Dxy, Dxz, Dyz = params 574 575 # Build the tensor. 576 tensor = zeros((3, 3), float64) 577 tensor[0, 0] = Dxx 578 tensor[1, 1] = Dyy 579 tensor[2, 2] = Dzz 580 tensor[0, 1] = tensor[1, 0] = Dxy 581 tensor[0, 2] = tensor[2, 0] = Dxz 582 tensor[1, 2] = tensor[2, 1] = Dyz 583 584 # Scaling. 585 tensor = tensor * d_scale 586 587 # Eigenvalues. 588 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 589 590 # Set the parameters. 591 set(value=[Di[0], Di[1], Di[2]], param=['Dx', 'Dy', 'Dz']) 592 593 # Change the angular units. 594 angle_units = 'rad' 595 596 # Unknown parameter combination. 597 else: 598 raise RelaxUnknownParamCombError('param_types', param_types) 599 600 # Convert the angles to radians. 601 if angle_units == 'deg': 602 print("Converting the angles from degrees to radian units.") 603 alpha = (alpha / 360.0) * 2.0 * pi 604 beta = (beta / 360.0) * 2.0 * pi 605 gamma = (gamma / 360.0) * 2.0 * pi 606 607 # Set the orientational parameters. 608 set(value=[alpha, beta, gamma], param=['alpha', 'beta', 'gamma'])

609 610

611 -def fold_angles(sim_index=None):

612 """Wrap the Euler or spherical angles and remove the glide reflection and translational symmetries. 613 614 Wrap the angles such that:: 615 616 0 <= theta <= pi, 617 0 <= phi <= 2pi, 618 619 and:: 620 621 0 <= alpha <= 2pi, 622 0 <= beta <= pi, 623 0 <= gamma <= 2pi. 624 625 626 For the simulated values, the angles are wrapped as:: 627 628 theta - pi/2 <= theta_sim <= theta + pi/2 629 phi - pi <= phi_sim <= phi + pi 630 631 and:: 632 633 alpha - pi <= alpha_sim <= alpha + pi 634 beta - pi/2 <= beta_sim <= beta + pi/2 635 gamma - pi <= gamma_sim <= gamma + pi 636 637 638 @param sim_index: The simulation index. If set to None then the actual values will be folded 639 rather than the simulation values. 640 @type sim_index: int or None 641 """ 642 643 644 # Wrap the angles. 645 ################## 646 647 # Spheroid. 648 if cdp.diff_tensor.type == 'spheroid': 649 # Get the current angles. 650 theta = cdp.diff_tensor.theta 651 phi = cdp.diff_tensor.phi 652 653 # Simulated values. 654 if sim_index != None: 655 theta_sim = cdp.diff_tensor.theta_sim[sim_index] 656 phi_sim = cdp.diff_tensor.phi_sim[sim_index] 657 658 # Normal value. 659 if sim_index == None: 660 cdp.diff_tensor.theta = wrap_angles(theta, 0.0, 2.0*pi) 661 cdp.diff_tensor.phi = wrap_angles(phi, 0.0, 2.0*pi) 662 663 # Simulated theta and phi values. 664 else: 665 cdp.diff_tensor.theta_sim[sim_index] = wrap_angles(theta_sim, theta - pi, theta + pi) 666 cdp.diff_tensor.phi_sim[sim_index] = wrap_angles(phi_sim, phi - pi, phi + pi) 667 668 # Ellipsoid. 669 elif cdp.diff_tensor.type == 'ellipsoid': 670 # Get the current angles. 671 alpha = cdp.diff_tensor.alpha 672 beta = cdp.diff_tensor.beta 673 gamma = cdp.diff_tensor.gamma 674 675 # Simulated values. 676 if sim_index != None: 677 alpha_sim = cdp.diff_tensor.alpha_sim[sim_index] 678 beta_sim = cdp.diff_tensor.beta_sim[sim_index] 679 gamma_sim = cdp.diff_tensor.gamma_sim[sim_index] 680 681 # Normal value. 682 if sim_index == None: 683 cdp.diff_tensor.alpha = wrap_angles(alpha, 0.0, 2.0*pi) 684 cdp.diff_tensor.beta = wrap_angles(beta, 0.0, 2.0*pi) 685 cdp.diff_tensor.gamma = wrap_angles(gamma, 0.0, 2.0*pi) 686 687 # Simulated alpha, beta, and gamma values. 688 else: 689 cdp.diff_tensor.alpha_sim[sim_index] = wrap_angles(alpha_sim, alpha - pi, alpha + pi) 690 cdp.diff_tensor.beta_sim[sim_index] = wrap_angles(beta_sim, beta - pi, beta + pi) 691 cdp.diff_tensor.gamma_sim[sim_index] = wrap_angles(gamma_sim, gamma - pi, gamma + pi) 692 693 694 # Remove the glide reflection and translational symmetries. 695 ########################################################### 696 697 # Spheroid. 698 if cdp.diff_tensor.type == 'spheroid': 699 # Normal value. 700 if sim_index == None: 701 # Fold phi inside 0 and pi. 702 if cdp.diff_tensor.phi >= pi: 703 theta, phi = fold_spherical_angles(cdp.diff_tensor.theta, cdp.diff_tensor.phi) 704 cdp.diff_tensor.theta = theta 705 cdp.diff_tensor.phi = phi 706 707 # Simulated theta and phi values. 708 else: 709 # Fold phi_sim inside phi-pi/2 and phi+pi/2. 710 if cdp.diff_tensor.phi_sim[sim_index] >= cdp.diff_tensor.phi + pi/2.0: 711 cdp.diff_tensor.theta_sim[sim_index] = pi - cdp.diff_tensor.theta_sim[sim_index] 712 cdp.diff_tensor.phi_sim[sim_index] = cdp.diff_tensor.phi_sim[sim_index] - pi 713 elif cdp.diff_tensor.phi_sim[sim_index] <= cdp.diff_tensor.phi - pi/2.0: 714 cdp.diff_tensor.theta_sim[sim_index] = pi - cdp.diff_tensor.theta_sim[sim_index] 715 cdp.diff_tensor.phi_sim[sim_index] = cdp.diff_tensor.phi_sim[sim_index] + pi 716 717 # Ellipsoid. 718 elif cdp.diff_tensor.type == 'ellipsoid': 719 # Normal value. 720 if sim_index == None: 721 # Fold alpha inside 0 and pi. 722 if cdp.diff_tensor.alpha >= pi: 723 cdp.diff_tensor.alpha = cdp.diff_tensor.alpha - pi 724 725 # Fold beta inside 0 and pi. 726 if cdp.diff_tensor.beta >= pi: 727 cdp.diff_tensor.alpha = pi - cdp.diff_tensor.alpha 728 cdp.diff_tensor.beta = cdp.diff_tensor.beta - pi 729 730 # Fold gamma inside 0 and pi. 731 if cdp.diff_tensor.gamma >= pi: 732 cdp.diff_tensor.alpha = pi - cdp.diff_tensor.alpha 733 cdp.diff_tensor.beta = pi - cdp.diff_tensor.beta 734 cdp.diff_tensor.gamma = cdp.diff_tensor.gamma - pi 735 736 # Simulated theta and phi values. 737 else: 738 # Fold alpha_sim inside alpha-pi/2 and alpha+pi/2. 739 if cdp.diff_tensor.alpha_sim[sim_index] >= cdp.diff_tensor.alpha + pi/2.0: 740 cdp.diff_tensor.alpha_sim[sim_index] = cdp.diff_tensor.alpha_sim[sim_index] - pi 741 elif cdp.diff_tensor.alpha_sim[sim_index] <= cdp.diff_tensor.alpha - pi/2.0: 742 cdp.diff_tensor.alpha_sim[sim_index] = cdp.diff_tensor.alpha_sim[sim_index] + pi 743 744 # Fold beta_sim inside beta-pi/2 and beta+pi/2. 745 if cdp.diff_tensor.beta_sim[sim_index] >= cdp.diff_tensor.beta + pi/2.0: 746 cdp.diff_tensor.alpha_sim[sim_index] = pi - cdp.diff_tensor.alpha_sim[sim_index] 747 cdp.diff_tensor.beta_sim[sim_index] = cdp.diff_tensor.beta_sim[sim_index] - pi 748 elif cdp.diff_tensor.beta_sim[sim_index] <= cdp.diff_tensor.beta - pi/2.0: 749 cdp.diff_tensor.alpha_sim[sim_index] = pi - cdp.diff_tensor.alpha_sim[sim_index] 750 cdp.diff_tensor.beta_sim[sim_index] = cdp.diff_tensor.beta_sim[sim_index] + pi 751 752 # Fold gamma_sim inside gamma-pi/2 and gamma+pi/2. 753 if cdp.diff_tensor.gamma_sim[sim_index] >= cdp.diff_tensor.gamma + pi/2.0: 754 cdp.diff_tensor.alpha_sim[sim_index] = pi - cdp.diff_tensor.alpha_sim[sim_index] 755 cdp.diff_tensor.beta_sim[sim_index] = pi - cdp.diff_tensor.beta_sim[sim_index] 756 cdp.diff_tensor.gamma_sim[sim_index] = cdp.diff_tensor.gamma_sim[sim_index] - pi 757 elif cdp.diff_tensor.gamma_sim[sim_index] <= cdp.diff_tensor.gamma - pi/2.0: 758 cdp.diff_tensor.alpha_sim[sim_index] = pi - cdp.diff_tensor.alpha_sim[sim_index] 759 cdp.diff_tensor.beta_sim[sim_index] = pi - cdp.diff_tensor.beta_sim[sim_index] 760 cdp.diff_tensor.gamma_sim[sim_index] = cdp.diff_tensor.gamma_sim[sim_index] + pi

761 762

763 -def init(params=None, time_scale=1.0, d_scale=1.0, angle_units='deg', param_types=0, spheroid_type=None, fixed=1):

764 """Function for initialising the diffusion tensor. 765 766 @param params: The diffusion tensor parameters. 767 @type params: float 768 @param time_scale: The correlation time scaling value. 769 @type time_scale: float 770 @param d_scale: The diffusion tensor eigenvalue scaling value. 771 @type d_scale: float 772 @param angle_units: The units for the angle parameters. 773 @type angle_units: str (either 'deg' or 'rad') 774 @param param_types: The type of parameters supplied. For spherical diffusion, the flag 775 values correspond to 0: tm, 1: Diso. For spheroidal diffusion, 0: {tm, 776 Da, theta, phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 777 3: {Dpar, Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. For 778 ellipsoidal diffusion, 0: {tm, Da, Dr, alpha, beta, gamma}, 1: {Diso, 779 Da, Dr, alpha, beta, gamma}, 2: {Dx, Dy, Dz, alpha, beta, gamma}. 780 @type param_types: int 781 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 782 restrict the tensor to either being 'oblate' or 'prolate'. 783 @type spheroid_type: str 784 @param fixed: A flag specifying whether the diffusion tensor is fixed or can be 785 optimised. 786 @type fixed: bin 787 """ 788 789 # Test if the current data pipe exists. 790 pipes.test() 791 792 # Test if diffusion tensor data already exists. 793 if diff_data_exists(): 794 raise RelaxTensorError('diffusion') 795 796 # Check the validity of the angle_units argument. 797 valid_types = ['deg', 'rad'] 798 if not angle_units in valid_types: 799 raise RelaxError("The diffusion tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.") 800 801 # Add the diff_tensor object to the data pipe. 802 cdp.diff_tensor = DiffTensorData() 803 804 # Set the fixed flag. 805 cdp.diff_tensor.fixed = fixed 806 807 # Spherical diffusion. 808 if isinstance(params, float): 809 num_params = 1 810 sphere(params, time_scale, param_types) 811 812 # Spheroidal diffusion. 813 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 4: 814 num_params = 4 815 spheroid(params, time_scale, d_scale, angle_units, param_types, spheroid_type) 816 817 # Ellipsoidal diffusion. 818 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 6: 819 num_params = 6 820 ellipsoid(params, time_scale, d_scale, angle_units, param_types) 821 822 # Unknown. 823 else: 824 raise RelaxError("The diffusion tensor parameters " + repr(params) + " are of an unknown type.") 825 826 # Test the validity of the parameters. 827 test_params(num_params)

828 829

830 -def map_bounds(param, spin_id=None):

831 """The function for creating bounds for the mapping function. 832 833 @param param: The name of the parameter to return the bounds for. 834 @type param: str 835 @keyword spin_id: The spin identification string. This arg is unused. 836 @type spin_id: None or str 837 @return: The bounds for the parameter. 838 @rtype: list of len 2 of floats 839 """ 840 841 # tm. 842 if param == 'tm': 843 return [0, 10.0 * 1e-9] 844 845 # {Diso, Dx, Dy, Dz, Dpar, Dper}. 846 if param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 847 return [1e6, 1e7] 848 849 # Da. 850 if param == 'Da': 851 return [-3.0/2.0 * 1e7, 3.0 * 1e7] 852 853 # Dr. 854 if param == 'Dr': 855 return [0, 1] 856 857 # Dratio. 858 if param == 'Dratio': 859 return [1.0/3.0, 3.0] 860 861 # theta. 862 if param == 'theta': 863 return [0, pi] 864 865 # phi. 866 if param == 'phi': 867 return [0, 2*pi] 868 869 # alpha. 870 if param == 'alpha': 871 return [0, 2*pi] 872 873 # beta. 874 if param == 'beta': 875 return [0, pi] 876 877 # gamma. 878 if param == 'gamma': 879 return [0, 2*pi]

880 881

882 -def map_labels(index, params, bounds, swap, inc):

883 """Function for creating labels, tick locations, and tick values for an OpenDX map.""" 884 885 # Initialise. 886 labels = "{" 887 tick_locations = [] 888 tick_values = [] 889 n = len(params) 890 axis_incs = 5 891 loc_inc = inc / axis_incs 892 893 # Increment over the model parameters. 894 for i in xrange(n): 895 # Parameter conversion factors. 896 factor = return_conversion_factor(params[swap[i]]) 897 898 # Parameter units. 899 units = return_units(params[swap[i]]) 900 901 # Labels. 902 if units: 903 labels = labels + "\"" + params[swap[i]] + " (" + units + ")\"" 904 else: 905 labels = labels + "\"" + params[swap[i]] + "\"" 906 907 # Tick values. 908 vals = bounds[swap[i], 0] / factor 909 val_inc = (bounds[swap[i], 1] - bounds[swap[i], 0]) / (axis_incs * factor) 910 911 if i < n - 1: 912 labels = labels + " " 913 else: 914 labels = labels + "}" 915 916 # Tick locations. 917 string = "{" 918 val = 0.0 919 for j in xrange(axis_incs + 1): 920 string = string + " " + repr(val) 921 val = val + loc_inc 922 string = string + " }" 923 tick_locations.append(string) 924 925 # Tick values. 926 string = "{" 927 for j in xrange(axis_incs + 1): 928 string = string + "\"" + "%.2f" % vals + "\" " 929 vals = vals + val_inc 930 string = string + "}" 931 tick_values.append(string) 932 933 return labels, tick_locations, tick_values

934 935

936 -def return_conversion_factor(param):

937 """Function for returning the factor of conversion between different parameter units. 938 939 For example, the internal representation of tm is in seconds, whereas the external 940 representation is in nanoseconds, therefore this function will return 1e-9 for tm. 941 942 943 @param param: The name of the parameter to return the conversion factor for. 944 @type param: str 945 @return: The conversion factor. 946 @rtype: float 947 """ 948 949 # Get the object name. 950 object_name = return_data_name(param) 951 952 # tm (nanoseconds). 953 if object_name == 'tm': 954 return 1e-9 955 956 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 957 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 958 return 1e6 959 960 # Angles. 961 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 962 return (2.0*pi) / 360.0 963 964 # No conversion factor. 965 return 1.0

966 967

968 -def return_data_name(name):

969 """Return the parameter name. 970 971 @param name: The name of the parameter to return the name of. 972 @type name: str 973 @return: The parameter name. 974 @rtype: str 975 """ 976 977 # Enforce that the name must be a string. 978 if not isinstance(name, str): 979 raise RelaxStrError('name', name) 980 981 # Local tm. 982 if search('^tm$', name): 983 return 'tm' 984 985 # Diso. 986 if search('[Dd]iso', name): 987 return 'Diso' 988 989 # Da. 990 if search('[Dd]a', name): 991 return 'Da' 992 993 # Dr. 994 if search('[Dd]r$', name): 995 return 'Dr' 996 997 # Dx. 998 if search('[Dd]x', name): 999 return 'Dx' 1000 1001 # Dy. 1002 if search('[Dd]y', name): 1003 return 'Dy' 1004 1005 # Dz. 1006 if search('[Dd]z', name): 1007 return 'Dz' 1008 1009 # Dpar. 1010 if search('[Dd]par', name): 1011 return 'Dpar' 1012 1013 # Dper. 1014 if search('[Dd]per', name): 1015 return 'Dper' 1016 1017 # Dratio. 1018 if search('[Dd]ratio', name): 1019 return 'Dratio' 1020 1021 # alpha. 1022 if search('^a$', name) or search('alpha', name): 1023 return 'alpha' 1024 1025 # beta. 1026 if search('^b$', name) or search('beta', name): 1027 return 'beta' 1028 1029 # gamma. 1030 if search('^g$', name) or search('gamma', name): 1031 return 'gamma' 1032 1033 # theta. 1034 if search('theta', name): 1035 return 'theta' 1036 1037 # phi. 1038 if search('phi', name): 1039 return 'phi'

1040 1041 # User function documentation. 1042 __return_data_name_prompt_doc__ = ["Diffusion tensor parameter string matching patterns", """ 1043 ____________________________________________________________________________________________ 1044 | | | | 1045 | Data type | Object name | Patterns | 1046 |________________________________________________________|______________|__________________| 1047 | | | | 1048 | Global correlation time - tm | 'tm' | '^tm$' | 1049 | | | | 1050 | Isotropic component of the diffusion tensor - Diso | 'Diso' | '[Dd]iso' | 1051 | | | | 1052 | Anisotropic component of the diffusion tensor - Da | 'Da' | '[Dd]a' | 1053 | | | | 1054 | Rhombic component of the diffusion tensor - Dr | 'Dr' | '[Dd]r$' | 1055 | | | | 1056 | Eigenvalue associated with the x-axis of the diffusion | 'Dx' | '[Dd]x' | 1057 | diffusion tensor - Dx | | | 1058 | | | | 1059 | Eigenvalue associated with the y-axis of the diffusion | 'Dy' | '[Dd]y' | 1060 | diffusion tensor - Dy | | | 1061 | | | | 1062 | Eigenvalue associated with the z-axis of the diffusion | 'Dz' | '[Dd]z' | 1063 | diffusion tensor - Dz | | | 1064 | | | | 1065 | Diffusion coefficient parallel to the major axis of | 'Dpar' | '[Dd]par' | 1066 | the spheroid diffusion tensor - Dpar | | | 1067 | | | | 1068 | Diffusion coefficient perpendicular to the major axis | 'Dper' | '[Dd]per' | 1069 | of the spheroid diffusion tensor - Dper | | | 1070 | | | | 1071 | Ratio of the parallel and perpendicular components of | 'Dratio' | '[Dd]ratio' | 1072 | the spheroid diffusion tensor - Dratio | | | 1073 | | | | 1074 | The first Euler angle of the ellipsoid diffusion | 'alpha' | '^a$' or 'alpha' | 1075 | tensor - alpha | | | 1076 | | | | 1077 | The second Euler angle of the ellipsoid diffusion | 'beta' | '^b$' or 'beta' | 1078 | tensor - beta | | | 1079 | | | | 1080 | The third Euler angle of the ellipsoid diffusion | 'gamma' | '^g$' or 'gamma' | 1081 | tensor - gamma | | | 1082 | | | | 1083 | The polar angle defining the major axis of the | 'theta' | 'theta' | 1084 | spheroid diffusion tensor - theta | | | 1085 | | | | 1086 | The azimuthal angle defining the major axis of the | 'phi' | 'phi' | 1087 | spheroid diffusion tensor - phi | | | 1088 |________________________________________________________|______________|__________________| 1089 """] 1090 1091

1092 -def return_eigenvalues():

1093 """Function for returning Dx, Dy, and Dz.""" 1094 1095 # Reassign the data. 1096 data = cdp.diff_tensor 1097 1098 # Diso. 1099 Diso = 1.0 / (6.0 * data.tm) 1100 1101 # Dx. 1102 Dx = Diso - 1.0/3.0 * data.Da * (1.0 + 3.0 * data.Dr) 1103 1104 # Dy. 1105 Dy = Diso - 1.0/3.0 * data.Da * (1.0 - 3.0 * data.Dr) 1106 1107 # Dz. 1108 Dz = Diso + 2.0/3.0 * data.Da 1109 1110 # Return the eigenvalues. 1111 return Dx, Dy, Dz

1112 1113

1114 -def return_units(param):

1115 """Function for returning a string representing the parameters units. 1116 1117 For example, the internal representation of tm is in seconds, whereas the external 1118 representation is in nanoseconds, therefore this function will return the string 1119 'nanoseconds' for tm. 1120 1121 1122 @param param: The name of the parameter to return the units for. 1123 @type param: str 1124 @return: The parameter units string. 1125 @rtype: str 1126 """ 1127 1128 # Get the object name. 1129 object_name = return_data_name(param) 1130 1131 # tm (nanoseconds). 1132 if object_name == 'tm': 1133 return 'ns' 1134 1135 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 1136 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 1137 return '1e6 1/s' 1138 1139 # Angles. 1140 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1141 return 'deg'

1142 1143

1144 -def set(value=None, param=None):

1145 """Set the diffusion tensor parameters. 1146 1147 @keyword tensor: The diffusion tensor object. 1148 @type tensor: DiffTensorData instance 1149 @keyword value: The list of values to set the parameters to. 1150 @type value: list of float 1151 @keyword param: The list of parameter names. 1152 @type param: list of str 1153 """ 1154 1155 # Set up the diffusion tensor data if it doesn't exist. 1156 if not diff_data_exists(): 1157 raise RelaxNoTensorError('diffusion') 1158 1159 # Initialise. 1160 geo_params = [] 1161 geo_values = [] 1162 orient_params = [] 1163 orient_values = [] 1164 1165 # Loop over the parameters. 1166 for i in xrange(len(param)): 1167 # Get the object name. 1168 param[i] = return_data_name(param[i]) 1169 1170 # Unknown parameter. 1171 if not param[i]: 1172 raise RelaxUnknownParamError("diffusion tensor", param[i]) 1173 1174 # Default value. 1175 if value[i] == None: 1176 value[i] = default_value(param[i]) 1177 1178 # Geometric parameter. 1179 if param[i] in ['tm', 'Diso', 'Da', 'Dratio', 'Dper', 'Dpar', 'Dr', 'Dx', 'Dy', 'Dz']: 1180 geo_params.append(param[i]) 1181 geo_values.append(value[i]) 1182 1183 # Orientational parameter. 1184 if param[i] in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1185 orient_params.append(param[i]) 1186 orient_values.append(value[i]) 1187 1188 1189 # Spherical diffusion. 1190 ###################### 1191 1192 if cdp.diff_tensor.type == 'sphere': 1193 # Geometric parameters. 1194 ####################### 1195 1196 # A single geometric parameter. 1197 if len(geo_params) == 1: 1198 # The single parameter tm. 1199 if geo_params[0] == 'tm': 1200 cdp.diff_tensor.tm = geo_values[0] 1201 1202 # The single parameter Diso. 1203 elif geo_params[0] == 'Diso': 1204 cdp.diff_tensor.tm = 1.0 / (6.0 * geo_values[0]) 1205 1206 # Cannot set the single parameter. 1207 else: 1208 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1209 1210 # More than one geometric parameters. 1211 elif len(geo_params) > 1: 1212 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1213 1214 1215 # Orientational parameters. 1216 ########################### 1217 1218 # ??? 1219 if len(orient_params): 1220 raise RelaxError("For spherical diffusion, the orientation parameters " + repr(orient_params) + " should not exist.") 1221 1222 1223 # Spheroidal diffusion. 1224 ####################### 1225 1226 elif cdp.diff_tensor.type == 'spheroid': 1227 # Geometric parameters. 1228 ####################### 1229 1230 # A single geometric parameter. 1231 if len(geo_params) == 1: 1232 # The single parameter tm. 1233 if geo_params[0] == 'tm': 1234 cdp.diff_tensor.tm = geo_values[0] 1235 1236 # The single parameter Diso. 1237 elif geo_params[0] == 'Diso': 1238 cdp.diff_tensor.tm = 1.0 / (6.0 * geo_values[0]) 1239 1240 # The single parameter Da. 1241 elif geo_params[0] == 'Da': 1242 cdp.diff_tensor.Da = geo_values[0] 1243 1244 # The single parameter Dratio. 1245 elif geo_params[0] == 'Dratio': 1246 Dratio = geo_values[0] 1247 cdp.diff_tensor.Da = (Dratio - 1.0) / (2.0 * cdp.diff_tensor.tm * (Dratio + 2.0)) 1248 1249 # Cannot set the single parameter. 1250 else: 1251 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1252 1253 # Two geometric parameters. 1254 elif len(geo_params) == 2: 1255 # The geometric parameter set {tm, Da}. 1256 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1257 # The parameters. 1258 tm = geo_values[geo_params.index('tm')] 1259 Da = geo_values[geo_params.index('Da')] 1260 1261 # Set the internal parameter values. 1262 cdp.diff_tensor.tm = tm 1263 cdp.diff_tensor.Da = Da 1264 1265 # The geometric parameter set {Diso, Da}. 1266 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1267 # The parameters. 1268 Diso = geo_values[geo_params.index('Diso')] 1269 Da = geo_values[geo_params.index('Da')] 1270 1271 # Set the internal parameter values. 1272 cdp.diff_tensor.tm = 1.0 / (6.0 * Diso) 1273 cdp.diff_tensor.Da = Da 1274 1275 # The geometric parameter set {tm, Dratio}. 1276 elif geo_params.count('tm') == 1 and geo_params.count('Dratio') == 1: 1277 # The parameters. 1278 tm = geo_values[geo_params.index('tm')] 1279 Dratio = geo_values[geo_params.index('Dratio')] 1280 1281 # Set the internal parameter values. 1282 cdp.diff_tensor.tm = tm 1283 cdp.diff_tensor.Da = (Dratio - 1.0) / (2.0 * tm * (Dratio + 2.0)) 1284 1285 # The geometric parameter set {Dpar, Dper}. 1286 elif geo_params.count('Dpar') == 1 and geo_params.count('Dper') == 1: 1287 # The parameters. 1288 Dpar = geo_values[geo_params.index('Dpar')] 1289 Dper = geo_values[geo_params.index('Dper')] 1290 1291 # Set the internal parameter values. 1292 cdp.diff_tensor.tm = 1.0 / (2.0 * (Dpar + 2.0*Dper)) 1293 cdp.diff_tensor.Da = Dpar - Dper 1294 1295 # The geometric parameter set {Diso, Dratio}. 1296 elif geo_params.count('Diso') == 1 and geo_params.count('Dratio') == 1: 1297 # The parameters. 1298 Diso = geo_values[geo_params.index('Diso')] 1299 Dratio = geo_values[geo_params.index('Dratio')] 1300 1301 # Set the internal parameter values. 1302 cdp.diff_tensor.tm = 1.0 / (6.0 * Diso) 1303 cdp.diff_tensor.Da = 3.0 * Diso * (Dratio - 1.0) / (Dratio + 2.0) 1304 1305 # Unknown parameter combination. 1306 else: 1307 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1308 1309 # More than two geometric parameters. 1310 elif len(geo_params) > 2: 1311 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1312 1313 1314 # Orientational parameters. 1315 ########################### 1316 1317 # A single orientational parameter. 1318 if len(orient_params) == 1: 1319 # The single parameter theta. 1320 if orient_params[0] == 'theta': 1321 cdp.diff_tensor.theta = orient_values[orient_params.index('theta')] 1322 1323 # The single parameter phi. 1324 elif orient_params[0] == 'phi': 1325 cdp.diff_tensor.phi = orient_values[orient_params.index('phi')] 1326 1327 # Disallowed parameter. 1328 else: 1329 raise RelaxError("For spheroidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1330 1331 # Two orientational parameters. 1332 elif len(orient_params) == 2: 1333 # The orientational parameter set {theta, phi}. 1334 if orient_params.count('theta') == 1 and orient_params.count('phi') == 1: 1335 cdp.diff_tensor.theta = orient_values[orient_params.index('theta')] 1336 cdp.diff_tensor.phi = orient_values[orient_params.index('phi')] 1337 1338 # Unknown parameter combination. 1339 else: 1340 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1341 1342 # More than two orientational parameters. 1343 elif len(orient_params) > 2: 1344 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1345 1346 1347 # Ellipsoidal diffusion. 1348 ######################## 1349 1350 elif cdp.diff_tensor.type == 'ellipsoid': 1351 # Geometric parameters. 1352 ####################### 1353 1354 # A single geometric parameter. 1355 if len(geo_params) == 1: 1356 # The single parameter tm. 1357 if geo_params[0] == 'tm': 1358 cdp.diff_tensor.tm = geo_values[0] 1359 1360 # The single parameter Diso. 1361 elif geo_params[0] == 'Diso': 1362 cdp.diff_tensor.tm = 1.0 / (6.0 * geo_values[0]) 1363 1364 # The single parameter Da. 1365 elif geo_params[0] == 'Da': 1366 cdp.diff_tensor.Da = geo_values[0] 1367 1368 # The single parameter Dr. 1369 elif geo_params[0] == 'Dr': 1370 cdp.diff_tensor.Dr = geo_values[0] 1371 1372 # Cannot set the single parameter. 1373 else: 1374 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1375 1376 # Two geometric parameters. 1377 elif len(geo_params) == 2: 1378 # The geometric parameter set {tm, Da}. 1379 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1380 # The parameters. 1381 tm = geo_values[geo_params.index('tm')] 1382 Da = geo_values[geo_params.index('Da')] 1383 1384 # Set the internal parameter values. 1385 cdp.diff_tensor.tm = tm 1386 cdp.diff_tensor.Da = Da 1387 1388 # The geometric parameter set {tm, Dr}. 1389 elif geo_params.count('tm') == 1 and geo_params.count('Dr') == 1: 1390 # The parameters. 1391 tm = geo_values[geo_params.index('tm')] 1392 Dr = geo_values[geo_params.index('Dr')] 1393 1394 # Set the internal parameter values. 1395 cdp.diff_tensor.tm = tm 1396 cdp.diff_tensor.Dr = Dr 1397 1398 # The geometric parameter set {Diso, Da}. 1399 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1400 # The parameters. 1401 Diso = geo_values[geo_params.index('Diso')] 1402 Da = geo_values[geo_params.index('Da')] 1403 1404 # Set the internal parameter values. 1405 cdp.diff_tensor.tm = 1.0 / (6.0 * Diso) 1406 cdp.diff_tensor.Da = Da 1407 1408 # The geometric parameter set {Diso, Dr}. 1409 elif geo_params.count('Diso') == 1 and geo_params.count('Dr') == 1: 1410 # The parameters. 1411 Diso = geo_values[geo_params.index('Diso')] 1412 Dr = geo_values[geo_params.index('Dr')] 1413 1414 # Set the internal parameter values. 1415 cdp.diff_tensor.tm = 1.0 / (6.0 * Diso) 1416 cdp.diff_tensor.Dr = Dr 1417 1418 # The geometric parameter set {Da, Dr}. 1419 elif geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1420 # The parameters. 1421 Da = geo_values[geo_params.index('Da')] 1422 Dr = geo_values[geo_params.index('Dr')] 1423 1424 # Set the internal parameter values. 1425 cdp.diff_tensor.Da = Da 1426 cdp.diff_tensor.Da = Dr 1427 1428 # Unknown parameter combination. 1429 else: 1430 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1431 1432 # Three geometric parameters. 1433 elif len(geo_params) == 3: 1434 # The geometric parameter set {tm, Da, Dr}. 1435 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1436 # The parameters. 1437 tm = geo_values[geo_params.index('tm')] 1438 Da = geo_values[geo_params.index('Da')] 1439 Dr = geo_values[geo_params.index('Dr')] 1440 1441 # Set the internal parameter values. 1442 cdp.diff_tensor.tm = tm 1443 cdp.diff_tensor.Da = Da 1444 cdp.diff_tensor.Dr = Dr 1445 1446 # The geometric parameter set {Diso, Da, Dr}. 1447 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1448 # The parameters. 1449 Diso = geo_values[geo_params.index('Diso')] 1450 Da = geo_values[geo_params.index('Da')] 1451 Dr = geo_values[geo_params.index('Dr')] 1452 1453 # Set the internal parameter values. 1454 cdp.diff_tensor.tm = 1.0 / (6.0 * Diso) 1455 cdp.diff_tensor.Da = Da 1456 cdp.diff_tensor.Dr = Dr 1457 1458 # The geometric parameter set {Dx, Dy, Dz}. 1459 elif geo_params.count('Dx') == 1 and geo_params.count('Dy') == 1 and geo_params.count('Dz') == 1: 1460 # The parameters. 1461 Dx = geo_values[geo_params.index('Dx')] 1462 Dy = geo_values[geo_params.index('Dy')] 1463 Dz = geo_values[geo_params.index('Dz')] 1464 1465 # Set the internal tm value. 1466 if Dx + Dy + Dz == 0.0: 1467 cdp.diff_tensor.tm = 1e99 1468 else: 1469 cdp.diff_tensor.tm = 0.5 / (Dx + Dy + Dz) 1470 1471 # Set the internal Da value. 1472 cdp.diff_tensor.Da = Dz - 0.5*(Dx + Dy) 1473 1474 # Set the internal Dr value. 1475 if cdp.diff_tensor.Da == 0.0: 1476 cdp.diff_tensor.Dr = (Dy - Dx) * 1e99 1477 else: 1478 cdp.diff_tensor.Dr = (Dy - Dx) / (2.0*cdp.diff_tensor.Da) 1479 1480 # Unknown parameter combination. 1481 else: 1482 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1483 1484 1485 # More than three geometric parameters. 1486 elif len(geo_params) > 3: 1487 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1488 1489 1490 # Orientational parameters. 1491 ########################### 1492 1493 # A single orientational parameter. 1494 if len(orient_params) == 1: 1495 # The single parameter alpha. 1496 if orient_params[0] == 'alpha': 1497 cdp.diff_tensor.alpha = orient_values[orient_params.index('alpha')] 1498 1499 # The single parameter beta. 1500 elif orient_params[0] == 'beta': 1501 cdp.diff_tensor.beta = orient_values[orient_params.index('beta')] 1502 1503 # The single parameter gamma. 1504 elif orient_params[0] == 'gamma': 1505 cdp.diff_tensor.gamma = orient_values[orient_params.index('gamma')] 1506 1507 # Disallowed parameter. 1508 else: 1509 raise RelaxError("For ellipsoidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1510 1511 # Two orientational parameters. 1512 elif len(orient_params) == 2: 1513 # The orientational parameter set {alpha, beta}. 1514 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1515 cdp.diff_tensor.alpha = orient_values[orient_params.index('alpha')] 1516 cdp.diff_tensor.beta = orient_values[orient_params.index('beta')] 1517 1518 # The orientational parameter set {alpha, gamma}. 1519 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1: 1520 cdp.diff_tensor.alpha = orient_values[orient_params.index('alpha')] 1521 cdp.diff_tensor.gamma = orient_values[orient_params.index('gamma')] 1522 1523 # The orientational parameter set {beta, gamma}. 1524 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1: 1525 cdp.diff_tensor.beta = orient_values[orient_params.index('beta')] 1526 cdp.diff_tensor.gamma = orient_values[orient_params.index('gamma')] 1527 1528 # Unknown parameter combination. 1529 else: 1530 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1531 1532 # Three orientational parameters. 1533 elif len(orient_params) == 3: 1534 # The orientational parameter set {alpha, beta, gamma}. 1535 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1536 cdp.diff_tensor.alpha = orient_values[orient_params.index('alpha')] 1537 cdp.diff_tensor.beta = orient_values[orient_params.index('beta')] 1538 cdp.diff_tensor.gamma = orient_values[orient_params.index('gamma')] 1539 1540 # Unknown parameter combination. 1541 else: 1542 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1543 1544 # More than three orientational parameters. 1545 elif len(orient_params) > 3: 1546 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1547 1548 1549 # Fold the angles in. 1550 ##################### 1551 1552 if orient_params: 1553 fold_angles()

1554 1555 # User function documentation. 1556 __set_prompt_doc__ = ["Diffusion tensor set details", """ 1557 If the diffusion tensor has not been setup, use the more powerful function 'diffusion_tensor.init' to initialise the tensor parameters. This function cannot be used to initialise a diffusion tensor. 1558 1559 The units of the parameters are: 1560 1561 Inverse seconds for tm. 1562 Seconds for Diso, Da, Dx, Dy, Dz, Dpar, Dper. 1563 Unitless for Dratio and Dr. 1564 Radians for all angles (alpha, beta, gamma, theta, phi). 1565 1566 When setting a diffusion tensor parameter, the residue number has no effect. As the internal parameters of spherical diffusion are {tm}, spheroidal diffusion are {tm, Da, theta, phi}, and ellipsoidal diffusion are {tm, Da, Dr, alpha, beta, gamma}, supplying geometric parameters must be done in the following way. If a single geometric parameter is supplied, it must be one of tm, Diso, Da, Dr, or Dratio. For the parameters Dpar, Dper, Dx, Dy, and Dx, it is not possible to determine how to use the currently set values together with the supplied value to calculate the new internal parameters. For spheroidal diffusion, when supplying multiple geometric parameters, the set must belong to one of 1567 1568 {tm, Da}, 1569 {Diso, Da}, 1570 {tm, Dratio}, 1571 {Dpar, Dper}, 1572 {Diso, Dratio}, 1573 1574 where either theta, phi, or both orientational parameters can be additionally supplied. For ellipsoidal diffusion, again when supplying multiple geometric parameters, the set must belong to one of 1575 1576 {tm, Da, Dr}, 1577 {Diso, Da, Dr}, 1578 {Dx, Dy, Dz}, 1579 1580 where any number of the orientational parameters, alpha, beta, or gamma can be additionally supplied. 1581 """] 1582 1583

1584 -def sphere(params=None, time_scale=None, param_types=None):

1585 """Function for setting up a spherical diffusion tensor. 1586 1587 @param params: The diffusion tensor parameter. 1588 @type params: float 1589 @param time_scale: The correlation time scaling value. 1590 @type time_scale: float 1591 @param param_types: The type of parameter supplied. If 0, then the parameter is tm. If 1, then 1592 the parameter is Diso. 1593 @type param_types: int 1594 """ 1595 1596 # The diffusion type. 1597 cdp.diff_tensor.type = 'sphere' 1598 1599 # tm. 1600 if param_types == 0: 1601 # Scaling. 1602 tm = params * time_scale 1603 1604 # Set the parameters. 1605 set(value=[tm], param=['tm']) 1606 1607 # Diso 1608 elif param_types == 1: 1609 # Scaling. 1610 Diso = params * d_scale 1611 1612 # Set the parameters. 1613 set(value=[Diso], param=['Diso']) 1614 1615 # Unknown parameter combination. 1616 else: 1617 raise RelaxUnknownParamCombError('param_types', param_types)

1618 1619

1620 -def spheroid(params=None, time_scale=None, d_scale=None, angle_units=None, param_types=None, spheroid_type=None):

1621 """Function for setting up a spheroidal diffusion tensor. 1622 1623 @param params: The diffusion tensor parameter. 1624 @type params: float 1625 @param time_scale: The correlation time scaling value. 1626 @type time_scale: float 1627 @param d_scale: The diffusion tensor eigenvalue scaling value. 1628 @type d_scale: float 1629 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 1630 @type angle_units: str 1631 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 1632 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 1633 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 1634 @type param_types: int 1635 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 1636 restrict the tensor to either being 'oblate' or 'prolate'. 1637 @type spheroid_type: str 1638 """ 1639 1640 # The diffusion type. 1641 cdp.diff_tensor.type = 'spheroid' 1642 1643 # Spheroid diffusion type. 1644 allowed_types = [None, 'oblate', 'prolate'] 1645 if spheroid_type not in allowed_types: 1646 raise RelaxError("The 'spheroid_type' argument " + repr(spheroid_type) + " should be 'oblate', 'prolate', or None.") 1647 cdp.diff_tensor.spheroid_type = spheroid_type 1648 1649 # (tm, Da, theta, phi). 1650 if param_types == 0: 1651 # Unpack the tuple. 1652 tm, Da, theta, phi = params 1653 1654 # Scaling. 1655 tm = tm * time_scale 1656 Da = Da * d_scale 1657 1658 # Set the parameters. 1659 set(value=[tm, Da], param=['tm', 'Da']) 1660 1661 # (Diso, Da, theta, phi). 1662 elif param_types == 1: 1663 # Unpack the tuple. 1664 Diso, Da, theta, phi = params 1665 1666 # Scaling. 1667 Diso = Diso * d_scale 1668 Da = Da * d_scale 1669 1670 # Set the parameters. 1671 set(value=[Diso, Da], param=['Diso', 'Da']) 1672 1673 # (tm, Dratio, theta, phi). 1674 elif param_types == 2: 1675 # Unpack the tuple. 1676 tm, Dratio, theta, phi = params 1677 1678 # Scaling. 1679 tm = tm * time_scale 1680 1681 # Set the parameters. 1682 set(value=[tm, Dratio], param=['tm', 'Dratio']) 1683 1684 # (Dpar, Dper, theta, phi). 1685 elif param_types == 3: 1686 # Unpack the tuple. 1687 Dpar, Dper, theta, phi = params 1688 1689 # Scaling. 1690 Dpar = Dpar * d_scale 1691 Dper = Dper * d_scale 1692 1693 # Set the parameters. 1694 set(value=[Dpar, Dper], param=['Dpar', 'Dper']) 1695 1696 # (Diso, Dratio, theta, phi). 1697 elif param_types == 4: 1698 # Unpack the tuple. 1699 Diso, Dratio, theta, phi = params 1700 1701 # Scaling. 1702 Diso = Diso * d_scale 1703 1704 # Set the parameters. 1705 set(value=[Diso, Dratio], param=['Diso', 'Dratio']) 1706 1707 # Unknown parameter combination. 1708 else: 1709 raise RelaxUnknownParamCombError('param_types', param_types) 1710 1711 # Convert the angles to radians. 1712 if angle_units == 'deg': 1713 print("Converting the angles from degrees to radian units.") 1714 theta = (theta / 360.0) * 2.0 * pi 1715 phi = (phi / 360.0) * 2.0 * pi 1716 1717 # Set the orientational parameters. 1718 set(value=[theta, phi], param=['theta', 'phi'])

1719 1720

1721 -def tensor_eigen_system(tensor):

1722 """Determine the eigenvalues and vectors for the tensor, sorting the entries. 1723 1724 @return: The eigenvalues, rotation matrix, and the Euler angles in zyz notation. 1725 @rtype: 3D rank-1 array, 3D rank-2 array, float, float, float 1726 """ 1727 1728 # Eigenvalues. 1729 R, Di, A = svd(tensor) 1730 D_diag = zeros((3, 3), float64) 1731 for i in range(3): 1732 D_diag[i, i] = Di[i] 1733 1734 # Reordering structure. 1735 tup_struct = [] 1736 for i in range(3): 1737 tup_struct.append((i, Di[i])) 1738 1739 # The indices. 1740 reorder_data = sorted(tup_struct, key=itemgetter(1)) 1741 reorder = zeros(3, int) 1742 Di_sort = zeros(3, float) 1743 for i in range(3): 1744 reorder[i], Di_sort[i] = reorder_data[i] 1745 1746 # Reorder columns. 1747 R_new = zeros((3, 3), float64) 1748 for i in range(3): 1749 R_new[:, i] = R[:, reorder[i]] 1750 1751 # Switch from the left handed to right handed universes (if needed). 1752 if norm(cross(R_new[:, 0], R_new[:, 1]) - R_new[:, 2]) > 1e-7: 1753 R_new[:, 2] = -R_new[:, 2] 1754 1755 # Reverse the rotation. 1756 R_new = transpose(R_new) 1757 1758 # Euler angles (reverse rotation in the rotated axis system). 1759 gamma, beta, alpha = R_to_euler_zyz(R_new) 1760 1761 # Collapse the pi axis rotation symmetries. 1762 if alpha >= pi: 1763 alpha = alpha - pi 1764 if gamma >= pi: 1765 alpha = pi - alpha 1766 beta = pi - beta 1767 gamma = gamma - pi 1768 if beta >= pi: 1769 alpha = pi - alpha 1770 beta = beta - pi 1771 1772 # Return the values. 1773 return Di_sort, R_new, alpha, beta, gamma

1774 1775

1776 -def test_params(num_params):

1777 """Function for testing the validity of the input parameters.""" 1778 1779 # An allowable error to account for machine precision, optimisation quality, etc. 1780 error = 1e-4 1781 1782 # tm. 1783 tm = cdp.diff_tensor.tm 1784 if tm <= 0.0 or tm > 1e-6: 1785 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1786 1787 # Spheroid. 1788 if num_params == 4: 1789 # Parameters. 1790 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1791 Da = cdp.diff_tensor.Da 1792 1793 # Da. 1794 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1795 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1796 1797 # Ellipsoid. 1798 if num_params == 6: 1799 # Parameters. 1800 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1801 Da = cdp.diff_tensor.Da 1802 Dr = cdp.diff_tensor.Dr 1803 1804 # Da. 1805 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1806 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1807 1808 # Dr. 1809 if Dr < (0.0 - error) or Dr > (1.0 + error): 1810 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1811 1812

1813 -def unit_axes():

1814 """Function for calculating the unit axes of the diffusion tensor. 1815 1816 Spheroid 1817 ======== 1818 1819 The unit Dpar vector is:: 1820 1821 | sin(theta) * cos(phi) | 1822 Dpar = | sin(theta) * sin(phi) | 1823 | cos(theta) | 1824 1825 1826 Ellipsoid 1827 ========= 1828 1829 The unit Dx vector is:: 1830 1831 | -sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma) | 1832 Dx = | -sin(alpha) * cos(gamma) - cos(alpha) * cos(beta) * sin(gamma) | 1833 | cos(alpha) * sin(beta) | 1834 1835 The unit Dy vector is:: 1836 1837 | cos(alpha) * sin(gamma) + sin(alpha) * cos(beta) * cos(gamma) | 1838 Dy = | cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma) | 1839 | sin(alpha) * sin(beta) | 1840 1841 The unit Dz vector is:: 1842 1843 | -sin(beta) * cos(gamma) | 1844 Dz = | sin(beta) * sin(gamma) | 1845 | cos(beta) | 1846 1847 """ 1848 1849 # Spheroid. 1850 if cdp.diff_tensor.type == 'spheroid': 1851 # Initialise. 1852 Dpar = zeros(3, float64) 1853 1854 # Trig. 1855 sin_theta = sin(cdp.diff_tensor.theta) 1856 cos_theta = cos(cdp.diff_tensor.theta) 1857 sin_phi = sin(cdp.diff_tensor.phi) 1858 cos_phi = cos(cdp.diff_tensor.phi) 1859 1860 # Unit Dpar axis. 1861 Dpar[0] = sin_theta * cos_phi 1862 Dpar[1] = sin_theta * sin_phi 1863 Dpar[2] = cos_theta 1864 1865 # Return the vector. 1866 return Dpar 1867 1868 # Ellipsoid. 1869 if cdp.diff_tensor.type == 'ellipsoid': 1870 # Initialise. 1871 Dx = zeros(3, float64) 1872 Dy = zeros(3, float64) 1873 Dz = zeros(3, float64) 1874 1875 # Trig. 1876 sin_alpha = sin(cdp.diff_tensor.alpha) 1877 cos_alpha = cos(cdp.diff_tensor.alpha) 1878 sin_beta = sin(cdp.diff_tensor.beta) 1879 cos_beta = cos(cdp.diff_tensor.beta) 1880 sin_gamma = sin(cdp.diff_tensor.gamma) 1881 cos_gamma = cos(cdp.diff_tensor.gamma) 1882 1883 # Unit Dx axis. 1884 Dx[0] = -sin_alpha * sin_gamma + cos_alpha * cos_beta * cos_gamma 1885 Dx[1] = -sin_alpha * cos_gamma - cos_alpha * cos_beta * sin_gamma 1886 Dx[2] = cos_alpha * sin_beta 1887 1888 # Unit Dy axis. 1889 Dx[0] = cos_alpha * sin_gamma + sin_alpha * cos_beta * cos_gamma 1890 Dx[1] = cos_alpha * cos_gamma - sin_alpha * cos_beta * sin_gamma 1891 Dx[2] = sin_alpha * sin_beta 1892 1893 # Unit Dz axis. 1894 Dx[0] = -sin_beta * cos_gamma 1895 Dx[1] = sin_beta * sin_gamma 1896 Dx[2] = cos_beta 1897 1898 # Return the vectors. 1899 return Dx, Dy, Dz

1900

Source Code for Module generic_fns.diffusion_tensor