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 from user_functions.data import Uf_tables; uf_tables = Uf_tables() 46 from user_functions.objects import Desc_container 47 48

49 -def bmrb_read(star):

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

146 147

148 -def bmrb_write(star):

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

263 264 265

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

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

303 304

305 -def data_names():

306 """Function for returning a list of names of data structures associated with the sequence.""" 307 308 names = [ 'diff_type', 309 'diff_params' ] 310 311 return names

312 313

314 -def default_value(param):

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

342 343 # User function documentation. 344 __default_value_doc__ = Desc_container("Diffusion tensor parameter default values") 345 table = uf_tables.add_table(label="table: diff default values", caption="Diffusion tensor parameter default values.") 346 table.add_headings(["Data type", "Object name", "Value"]) 347 table.add_row(["tm", "'tm'", "10 * 1e-9"]) 348 table.add_row(["Diso", "'Diso'", "1.666 * 1e7"]) 349 table.add_row(["Da", "'Da'", "0.0"]) 350 table.add_row(["Dr", "'Dr'", "0.0"]) 351 table.add_row(["Dx", "'Dx'", "1.666 * 1e7"]) 352 table.add_row(["Dy", "'Dy'", "1.666 * 1e7"]) 353 table.add_row(["Dz", "'Dz'", "1.666 * 1e7"]) 354 table.add_row(["Dpar", "'Dpar'", "1.666 * 1e7"]) 355 table.add_row(["Dper", "'Dper'", "1.666 * 1e7"]) 356 table.add_row(["Dratio", "'Dratio'", "1.0"]) 357 table.add_row(["alpha", "'alpha'", "0.0"]) 358 table.add_row(["beta", "'beta'", "0.0"]) 359 table.add_row(["gamma", "'gamma'", "0.0"]) 360 table.add_row(["theta", "'theta'", "0.0"]) 361 table.add_row(["phi", "'phi'", "0.0"]) 362 __default_value_doc__.add_table(table.label) 363 364

365 -def delete():

366 """Function for deleting diffusion tensor data.""" 367 368 # Test if the current data pipe exists. 369 pipes.test() 370 371 # Test if diffusion tensor data exists. 372 if not diff_data_exists(): 373 raise RelaxNoTensorError('diffusion') 374 375 # Delete the diffusion data. 376 del(cdp.diff_tensor)

377 378

379 -def diff_data_exists(pipe=None):

380 """Function for determining if diffusion data exists in the current data pipe. 381 382 @param pipe: The data pipe to search for data in. 383 @type pipe: str 384 @return: The answer to the question. 385 @rtype: bool 386 """ 387 388 # The data pipe to check. 389 if pipe == None: 390 pipe = pipes.cdp_name() 391 392 # Get the data pipe. 393 dp = pipes.get_pipe(pipe) 394 395 # Test if the data structure exists. 396 if hasattr(dp, 'diff_tensor'): 397 return True 398 else: 399 return False

400 401

402 -def display():

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

498 499

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

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

596 597

598 -def fold_angles(sim_index=None):

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

748 749

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

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

815 816

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

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

867 868

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

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

921 922

923 -def return_conversion_factor(param):

924 """Function for returning the factor of conversion between different parameter units. 925 926 For example, the internal representation of tm is in seconds, whereas the external 927 representation is in nanoseconds, therefore this function will return 1e-9 for tm. 928 929 930 @param param: The name of the parameter to return the conversion factor for. 931 @type param: str 932 @return: The conversion factor. 933 @rtype: float 934 """ 935 936 # Get the object name. 937 object_name = return_data_name(param) 938 939 # tm (nanoseconds). 940 if object_name == 'tm': 941 return 1e-9 942 943 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 944 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 945 return 1e6 946 947 # Angles. 948 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 949 return (2.0*pi) / 360.0 950 951 # No conversion factor. 952 return 1.0

953 954

955 -def return_data_name(name):

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

1027 1028 # User function documentation. 1029 __return_data_name_doc__ = Desc_container("Diffusion tensor parameter string matching patterns") 1030 table = uf_tables.add_table(label="table: diff data type patterns", caption="Diffusion tensor parameter string matching patterns.") 1031 table.add_headings(["Data type", "Object name", "Patterns"]) 1032 table.add_row(["Global correlation time - tm", "'tm'", "'^tm$'"]) 1033 table.add_row(["Isotropic component of the diffusion tensor - Diso", "'Diso'", "'[Dd]iso'"]) 1034 table.add_row(["Anisotropic component of the diffusion tensor - Da", "'Da'", "'[Dd]a'"]) 1035 table.add_row(["Rhombic component of the diffusion tensor - Dr", "'Dr'", "'[Dd]r$'"]) 1036 table.add_row(["Eigenvalue associated with the x-axis of the diffusion tensor - Dx", "'Dx'", "'[Dd]x'"]) 1037 table.add_row(["Eigenvalue associated with the y-axis of the diffusion tensor - Dy", "'Dy'", "'[Dd]y'"]) 1038 table.add_row(["Eigenvalue associated with the z-axis of the diffusion tensor - Dz", "'Dz'", "'[Dd]z'"]) 1039 table.add_row(["Diffusion coefficient parallel to the major axis of the spheroid diffusion tensor - Dpar", "'Dpar'", "'[Dd]par'"]) 1040 table.add_row(["Diffusion coefficient perpendicular to the major axis of the spheroid diffusion tensor - Dper", "'Dper'", "'[Dd]per'"]) 1041 table.add_row(["Ratio of the parallel and perpendicular components of the spheroid diffusion tensor - Dratio", "'Dratio'", "'[Dd]ratio'"]) 1042 table.add_row(["The first Euler angle of the ellipsoid diffusion tensor - alpha", "'alpha'", "'^a$' or 'alpha'"]) 1043 table.add_row(["The second Euler angle of the ellipsoid diffusion tensor - beta", "'beta'", "'^b$' or 'beta'"]) 1044 table.add_row(["The third Euler angle of the ellipsoid diffusion tensor - gamma", "'gamma'", "'^g$' or 'gamma'"]) 1045 table.add_row(["The polar angle defining the major axis of the spheroid diffusion tensor - theta", "'theta'", "'theta'"]) 1046 table.add_row(["The azimuthal angle defining the major axis of the spheroid diffusion tensor - phi", "'phi'", "'phi'"]) 1047 __return_data_name_doc__.add_table(table.label) 1048 1049

1050 -def return_eigenvalues():

1051 """Function for returning Dx, Dy, and Dz.""" 1052 1053 # Reassign the data. 1054 data = cdp.diff_tensor 1055 1056 # Diso. 1057 Diso = 1.0 / (6.0 * data.tm) 1058 1059 # Dx. 1060 Dx = Diso - 1.0/3.0 * data.Da * (1.0 + 3.0 * data.Dr) 1061 1062 # Dy. 1063 Dy = Diso - 1.0/3.0 * data.Da * (1.0 - 3.0 * data.Dr) 1064 1065 # Dz. 1066 Dz = Diso + 2.0/3.0 * data.Da 1067 1068 # Return the eigenvalues. 1069 return Dx, Dy, Dz

1070 1071

1072 -def return_units(param):

1073 """Function for returning a string representing the parameters units. 1074 1075 For example, the internal representation of tm is in seconds, whereas the external 1076 representation is in nanoseconds, therefore this function will return the string 1077 'nanoseconds' for tm. 1078 1079 1080 @param param: The name of the parameter to return the units for. 1081 @type param: str 1082 @return: The parameter units string. 1083 @rtype: str 1084 """ 1085 1086 # Get the object name. 1087 object_name = return_data_name(param) 1088 1089 # tm (nanoseconds). 1090 if object_name == 'tm': 1091 return 'ns' 1092 1093 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 1094 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 1095 return '1e6 1/s' 1096 1097 # Angles. 1098 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1099 return 'deg'

1100 1101

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

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

1512 1513 # User function documentation. 1514 __set_doc__ = Desc_container("Diffusion tensor set details") 1515 __set_doc__.add_paragraph("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.") 1516 __set_doc__.add_paragraph("The units of the parameters are:") 1517 __set_doc__.add_list_element("Inverse seconds for tm.") 1518 __set_doc__.add_list_element("Seconds for Diso, Da, Dx, Dy, Dz, Dpar, Dper.") 1519 __set_doc__.add_list_element("Unitless for Dratio and Dr.") 1520 __set_doc__.add_list_element("Radians for all angles (alpha, beta, gamma, theta, phi).") 1521 __set_doc__.add_paragraph("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") 1522 __set_doc__.add_list_element("{tm, Da},") 1523 __set_doc__.add_list_element("{Diso, Da},") 1524 __set_doc__.add_list_element("{tm, Dratio},") 1525 __set_doc__.add_list_element("{Dpar, Dper},") 1526 __set_doc__.add_list_element("{Diso, Dratio},") 1527 __set_doc__.add_paragraph("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") 1528 __set_doc__.add_list_element("{tm, Da, Dr},") 1529 __set_doc__.add_list_element("{Diso, Da, Dr},") 1530 __set_doc__.add_list_element("{Dx, Dy, Dz},") 1531 __set_doc__.add_paragraph("where any number of the orientational parameters, alpha, beta, or gamma can be additionally supplied.") 1532 1533

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

1535 """Function for setting up a spherical diffusion tensor. 1536 1537 @param params: The diffusion tensor parameter. 1538 @type params: float 1539 @param time_scale: The correlation time scaling value. 1540 @type time_scale: float 1541 @param param_types: The type of parameter supplied. If 0, then the parameter is tm. If 1, then 1542 the parameter is Diso. 1543 @type param_types: int 1544 """ 1545 1546 # The diffusion type. 1547 cdp.diff_tensor.type = 'sphere' 1548 1549 # tm. 1550 if param_types == 0: 1551 # Scaling. 1552 tm = params * time_scale 1553 1554 # Set the parameters. 1555 set(value=[tm], param=['tm']) 1556 1557 # Diso 1558 elif param_types == 1: 1559 # Scaling. 1560 Diso = params * d_scale 1561 1562 # Set the parameters. 1563 set(value=[Diso], param=['Diso']) 1564 1565 # Unknown parameter combination. 1566 else: 1567 raise RelaxUnknownParamCombError('param_types', param_types)

1568 1569

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

1571 """Function for setting up a spheroidal diffusion tensor. 1572 1573 @param params: The diffusion tensor parameter. 1574 @type params: float 1575 @param time_scale: The correlation time scaling value. 1576 @type time_scale: float 1577 @param d_scale: The diffusion tensor eigenvalue scaling value. 1578 @type d_scale: float 1579 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 1580 @type angle_units: str 1581 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 1582 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 1583 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 1584 @type param_types: int 1585 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 1586 restrict the tensor to either being 'oblate' or 'prolate'. 1587 @type spheroid_type: str 1588 """ 1589 1590 # The diffusion type. 1591 cdp.diff_tensor.type = 'spheroid' 1592 1593 # Spheroid diffusion type. 1594 allowed_types = [None, 'oblate', 'prolate'] 1595 if spheroid_type not in allowed_types: 1596 raise RelaxError("The 'spheroid_type' argument " + repr(spheroid_type) + " should be 'oblate', 'prolate', or None.") 1597 cdp.diff_tensor.spheroid_type = spheroid_type 1598 1599 # (tm, Da, theta, phi). 1600 if param_types == 0: 1601 # Unpack the tuple. 1602 tm, Da, theta, phi = params 1603 1604 # Scaling. 1605 tm = tm * time_scale 1606 Da = Da * d_scale 1607 1608 # Set the parameters. 1609 set(value=[tm, Da], param=['tm', 'Da']) 1610 1611 # (Diso, Da, theta, phi). 1612 elif param_types == 1: 1613 # Unpack the tuple. 1614 Diso, Da, theta, phi = params 1615 1616 # Scaling. 1617 Diso = Diso * d_scale 1618 Da = Da * d_scale 1619 1620 # Set the parameters. 1621 set(value=[Diso, Da], param=['Diso', 'Da']) 1622 1623 # (tm, Dratio, theta, phi). 1624 elif param_types == 2: 1625 # Unpack the tuple. 1626 tm, Dratio, theta, phi = params 1627 1628 # Scaling. 1629 tm = tm * time_scale 1630 1631 # Set the parameters. 1632 set(value=[tm, Dratio], param=['tm', 'Dratio']) 1633 1634 # (Dpar, Dper, theta, phi). 1635 elif param_types == 3: 1636 # Unpack the tuple. 1637 Dpar, Dper, theta, phi = params 1638 1639 # Scaling. 1640 Dpar = Dpar * d_scale 1641 Dper = Dper * d_scale 1642 1643 # Set the parameters. 1644 set(value=[Dpar, Dper], param=['Dpar', 'Dper']) 1645 1646 # (Diso, Dratio, theta, phi). 1647 elif param_types == 4: 1648 # Unpack the tuple. 1649 Diso, Dratio, theta, phi = params 1650 1651 # Scaling. 1652 Diso = Diso * d_scale 1653 1654 # Set the parameters. 1655 set(value=[Diso, Dratio], param=['Diso', 'Dratio']) 1656 1657 # Unknown parameter combination. 1658 else: 1659 raise RelaxUnknownParamCombError('param_types', param_types) 1660 1661 # Convert the angles to radians. 1662 if angle_units == 'deg': 1663 print("Converting the angles from degrees to radian units.") 1664 theta = (theta / 360.0) * 2.0 * pi 1665 phi = (phi / 360.0) * 2.0 * pi 1666 1667 # Set the orientational parameters. 1668 set(value=[theta, phi], param=['theta', 'phi'])

1669 1670

1671 -def tensor_eigen_system(tensor):

1672 """Determine the eigenvalues and vectors for the tensor, sorting the entries. 1673 1674 @return: The eigenvalues, rotation matrix, and the Euler angles in zyz notation. 1675 @rtype: 3D rank-1 array, 3D rank-2 array, float, float, float 1676 """ 1677 1678 # Eigenvalues. 1679 R, Di, A = svd(tensor) 1680 D_diag = zeros((3, 3), float64) 1681 for i in range(3): 1682 D_diag[i, i] = Di[i] 1683 1684 # Reordering structure. 1685 tup_struct = [] 1686 for i in range(3): 1687 tup_struct.append((i, Di[i])) 1688 1689 # The indices. 1690 reorder_data = sorted(tup_struct, key=itemgetter(1)) 1691 reorder = zeros(3, int) 1692 Di_sort = zeros(3, float) 1693 for i in range(3): 1694 reorder[i], Di_sort[i] = reorder_data[i] 1695 1696 # Reorder columns. 1697 R_new = zeros((3, 3), float64) 1698 for i in range(3): 1699 R_new[:, i] = R[:, reorder[i]] 1700 1701 # Switch from the left handed to right handed universes (if needed). 1702 if norm(cross(R_new[:, 0], R_new[:, 1]) - R_new[:, 2]) > 1e-7: 1703 R_new[:, 2] = -R_new[:, 2] 1704 1705 # Reverse the rotation. 1706 R_new = transpose(R_new) 1707 1708 # Euler angles (reverse rotation in the rotated axis system). 1709 gamma, beta, alpha = R_to_euler_zyz(R_new) 1710 1711 # Collapse the pi axis rotation symmetries. 1712 if alpha >= pi: 1713 alpha = alpha - pi 1714 if gamma >= pi: 1715 alpha = pi - alpha 1716 beta = pi - beta 1717 gamma = gamma - pi 1718 if beta >= pi: 1719 alpha = pi - alpha 1720 beta = beta - pi 1721 1722 # Return the values. 1723 return Di_sort, R_new, alpha, beta, gamma

1724 1725

1726 -def test_params(num_params):

1727 """Function for testing the validity of the input parameters.""" 1728 1729 # An allowable error to account for machine precision, optimisation quality, etc. 1730 error = 1e-4 1731 1732 # tm. 1733 tm = cdp.diff_tensor.tm 1734 if tm <= 0.0 or tm > 1e-6: 1735 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1736 1737 # Spheroid. 1738 if num_params == 4: 1739 # Parameters. 1740 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1741 Da = cdp.diff_tensor.Da 1742 1743 # Da. 1744 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1745 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1746 1747 # Ellipsoid. 1748 if num_params == 6: 1749 # Parameters. 1750 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1751 Da = cdp.diff_tensor.Da 1752 Dr = cdp.diff_tensor.Dr 1753 1754 # Da. 1755 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1756 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1757 1758 # Dr. 1759 if Dr < (0.0 - error) or Dr > (1.0 + error): 1760 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1761 1762

1763 -def unit_axes():

1764 """Function for calculating the unit axes of the diffusion tensor. 1765 1766 Spheroid 1767 ======== 1768 1769 The unit Dpar vector is:: 1770 1771 | sin(theta) * cos(phi) | 1772 Dpar = | sin(theta) * sin(phi) | 1773 | cos(theta) | 1774 1775 1776 Ellipsoid 1777 ========= 1778 1779 The unit Dx vector is:: 1780 1781 | -sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma) | 1782 Dx = | -sin(alpha) * cos(gamma) - cos(alpha) * cos(beta) * sin(gamma) | 1783 | cos(alpha) * sin(beta) | 1784 1785 The unit Dy vector is:: 1786 1787 | cos(alpha) * sin(gamma) + sin(alpha) * cos(beta) * cos(gamma) | 1788 Dy = | cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma) | 1789 | sin(alpha) * sin(beta) | 1790 1791 The unit Dz vector is:: 1792 1793 | -sin(beta) * cos(gamma) | 1794 Dz = | sin(beta) * sin(gamma) | 1795 | cos(beta) | 1796 1797 """ 1798 1799 # Spheroid. 1800 if cdp.diff_tensor.type == 'spheroid': 1801 # Initialise. 1802 Dpar = zeros(3, float64) 1803 1804 # Trig. 1805 sin_theta = sin(cdp.diff_tensor.theta) 1806 cos_theta = cos(cdp.diff_tensor.theta) 1807 sin_phi = sin(cdp.diff_tensor.phi) 1808 cos_phi = cos(cdp.diff_tensor.phi) 1809 1810 # Unit Dpar axis. 1811 Dpar[0] = sin_theta * cos_phi 1812 Dpar[1] = sin_theta * sin_phi 1813 Dpar[2] = cos_theta 1814 1815 # Return the vector. 1816 return Dpar 1817 1818 # Ellipsoid. 1819 if cdp.diff_tensor.type == 'ellipsoid': 1820 # Initialise. 1821 Dx = zeros(3, float64) 1822 Dy = zeros(3, float64) 1823 Dz = zeros(3, float64) 1824 1825 # Trig. 1826 sin_alpha = sin(cdp.diff_tensor.alpha) 1827 cos_alpha = cos(cdp.diff_tensor.alpha) 1828 sin_beta = sin(cdp.diff_tensor.beta) 1829 cos_beta = cos(cdp.diff_tensor.beta) 1830 sin_gamma = sin(cdp.diff_tensor.gamma) 1831 cos_gamma = cos(cdp.diff_tensor.gamma) 1832 1833 # Unit Dx axis. 1834 Dx[0] = -sin_alpha * sin_gamma + cos_alpha * cos_beta * cos_gamma 1835 Dx[1] = -sin_alpha * cos_gamma - cos_alpha * cos_beta * sin_gamma 1836 Dx[2] = cos_alpha * sin_beta 1837 1838 # Unit Dy axis. 1839 Dx[0] = cos_alpha * sin_gamma + sin_alpha * cos_beta * cos_gamma 1840 Dx[1] = cos_alpha * cos_gamma - sin_alpha * cos_beta * sin_gamma 1841 Dx[2] = sin_alpha * sin_beta 1842 1843 # Unit Dz axis. 1844 Dx[0] = -sin_beta * cos_gamma 1845 Dx[1] = sin_beta * sin_gamma 1846 Dx[2] = cos_beta 1847 1848 # Return the vectors. 1849 return Dx, Dy, Dz

1850

Source Code for Module generic_fns.diffusion_tensor