maths_fns.correlation

1 ############################################################################### 2 # # 3 # Copyright (C) 2004 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 from math import sqrt 24 25 26 # Isotropic global correlation time equation. 27 ############################################# 28

29 -def calc_iso_ti(data, diff_data):

30 """Diffusional correlation times for isotropic diffusion. 31 32 t0 = tm 33 """ 34 35 data.ti[0] = diff_data.params[0]

36 37 38 39 # Isotropic global correlation time gradient. 40 ############################################# 41

42 -def calc_iso_dti(data, diff_data):

43 """Partial derivatives of the diffusional correlation times. 44 45 The tm partial derivatives are: 46 47 dt0 48 --- = 1 49 dtm 50 """ 51 52 data.dti[0] = 1.0

53 54 55 56 # Axially symmetric global correlation time equation. 57 ##################################################### 58

59 -def calc_axial_ti(data, diff_data):

60 """Diffusional correlation times. 61 62 The equations for the parameters {Dper, Dpar} are: 63 64 1 65 t0 = ----- 66 6Dper 67 68 1 69 t1 = ------------ 70 5Dper + Dpar 71 72 1 73 t2 = ------------- 74 2Dper + 4Dpar 75 76 77 The equations for the parameters {Diso, Da} are: 78 79 t0 = 1/6 (Diso - Da)**-1 80 81 t1 = 1/6 (Diso - Da/2)**-1 82 83 t2 = 1/6 (Diso + Da)**-1 84 85 The diffusion parameter set in data.diff_params is {tm, Da, theta, phi}. 86 """ 87 88 # Calculate Diso. 89 if diff_data.params[0] == 0.0: 90 data.Diso = 1e99 91 else: 92 data.Diso = 1.0 / (6.0 * diff_data.params[0]) 93 94 # Components. 95 data.t_0_comp = data.Diso - diff_data.params[1] 96 data.t_1_comp = data.Diso - 0.5 * diff_data.params[1] 97 data.t_2_comp = data.Diso + diff_data.params[1] 98 99 # t0. 100 data.ti[0] = 6.0 * data.t_0_comp 101 if data.ti[0] == 0.0: 102 data.ti[0] = 1e99 103 else: 104 data.ti[0] = 1.0 / data.ti[0] 105 106 # t1. 107 data.ti[1] = 6.0 * data.t_1_comp 108 if data.ti[1] == 0.0: 109 data.ti[1] = 1e99 110 else: 111 data.ti[1] = 1.0 / data.ti[1] 112 113 # t2. 114 data.ti[2] = 6.0 * data.t_2_comp 115 if data.ti[2] == 0.0: 116 data.ti[2] = 1e99 117 else: 118 data.ti[2] = 1.0 / data.ti[2]

119 120 121 122 123 # Axially symmetric global correlation time gradient. 124 ##################################################### 125

126 -def calc_axial_dti(data, diff_data):

127 """Diffusional correlation time gradients. 128 129 tm partial derivatives 130 ~~~~~~~~~~~~~~~~~~~~~~ 131 132 dt0 1 dDiso 133 --- = - - ----- (Diso - Da)**-2 134 dtm 6 dtm 135 136 dt1 1 dDiso 137 --- = - - ----- (Diso - Da/2)**-2 138 dtm 6 dtm 139 140 dt2 1 dDiso 141 --- = - - ----- (Diso + Da)**-2 142 dtm 6 dtm 143 144 145 dDiso 146 ----- = -1/6 * tm**-2 147 dtm 148 149 150 Da partial derivatives 151 ~~~~~~~~~~~~~~~~~~~~~~ 152 153 dt0 154 --- = 1/6 (Diso - Da)**-2 155 dDa 156 157 dt1 158 --- = 1/12 (Diso - Da/2)**-2 159 dDa 160 161 dt2 162 --- = -1/6 (Diso + Da)**-2 163 dDa 164 165 166 The diffusion parameter set in data.diff_params is {tm, Da, theta, phi}. 167 """ 168 169 # Components. 170 data.t_0_comp_sqrd = data.t_0_comp**2 171 data.t_1_comp_sqrd = data.t_1_comp**2 172 data.t_2_comp_sqrd = data.t_2_comp**2 173 174 175 # tm partial derivative. 176 ######################## 177 178 # Components. 179 data.inv_dDiso_dtm = -6.0 * diff_data.params[0]**2 180 181 # t0. 182 data.dti[0, 0] = -6.0 * data.inv_dDiso_dtm * data.t_0_comp_sqrd 183 if data.dti[0, 0] == 0.0: 184 data.dti[0, 0] = 1e99 185 else: 186 data.dti[0, 0] = 1.0 / data.dti[0, 0] 187 188 # t1. 189 data.dti[0, 1] = -6.0 * data.inv_dDiso_dtm * data.t_1_comp_sqrd 190 if data.dti[0, 1] == 0.0: 191 data.dti[0, 1] = 1e99 192 else: 193 data.dti[0, 1] = 1.0 / data.dti[0, 1] 194 195 # t2. 196 data.dti[0, 2] = -6.0 * data.inv_dDiso_dtm * data.t_2_comp_sqrd 197 if data.dti[0, 2] == 0.0: 198 data.dti[0, 2] = 1e99 199 else: 200 data.dti[0, 2] = 1.0 / data.dti[0, 2] 201 202 203 # Da partial derivative. 204 ######################## 205 206 # t0. 207 data.dti[1, 0] = 6.0 * data.t_0_comp_sqrd 208 if data.dti[1, 0] == 0.0: 209 if data.mu != 0.0: 210 data.dti[1, 0] = 1e99 211 else: 212 data.dti[1, 0] = 1.0 / data.dti[1, 0] 213 214 # t1. 215 data.dti[1, 1] = 12.0 * data.t_1_comp_sqrd 216 if data.dti[1, 1] == 0.0: 217 data.dti[1, 1] = 1e99 218 else: 219 data.dti[1, 1] = 1.0 / data.dti[1, 1] 220 221 # t2. 222 data.dti[1, 2] = -6.0 * data.t_2_comp_sqrd 223 if data.dti[1, 2] == 0.0: 224 if data.mu != 0.0: 225 data.dti[1, 2] = 1e99 226 else: 227 data.dti[1, 2] = 1.0 / data.dti[1, 2]

228 229 230 # Axially symmetric global correlation time Hessian. 231 #################################################### 232

233 -def calc_axial_d2ti(data, diff_data):

234 """Diffusional correlation time Hessians. 235 236 tm-tm partial derivatives 237 ~~~~~~~~~~~~~~~~~~~~~~~~~ 238 239 d2t0 1 / dDiso \ 2 1 d2Diso 240 ---- = - | ----- | (Diso - Da)**-3 - - ------ (Diso - Da)**-2 241 dtm2 3 \ dtm / 6 dtm2 242 243 d2t1 1 / dDiso \ 2 1 d2Diso 244 ---- = - | ----- | (Diso - Da/2)**-3 - - ------ (Diso - Da/2)**-2 245 dtm2 3 \ dtm / 6 dtm2 246 247 d2t2 1 / dDiso \ 2 1 d2Diso 248 ---- = - | ----- | (Diso + Da)**-3 - - ------ (Diso + Da)**-2 249 dtm2 3 \ dtm / 6 dtm2 250 251 252 d2Diso 253 ------ = 1/3 * tm**-3 254 dtm2 255 256 257 tm-Da partial derivatives 258 ~~~~~~~~~~~~~~~~~~~~~~~~~ 259 260 d2t0 1 dDiso 261 ------- = - - ----- (Diso - Da)**-3 262 dtm.dDa 3 dtm 263 264 d2t1 1 dDiso 265 ------- = - - ----- (Diso - Da/2)**-3 266 dtm.dDa 6 dtm 267 268 d2t2 1 dDiso 269 ------- = - ----- (Diso + Da)**-3 270 dtm.dDa 3 dtm 271 272 273 Da-Da partial derivatives 274 ~~~~~~~~~~~~~~~~~~~~~~~~~ 275 276 d2t0 277 ---- = 1/3 (Diso - Da)**-3 278 dDa2 279 280 d2t1 281 ---- = 1/12 (Diso - Da/2)**-3 282 dDa2 283 284 d2t2 285 ---- = 1/3 (Diso + Da)**-3 286 dDa2 287 288 289 The diffusion parameter set in data.diff_params is {tm, Da, theta, phi}. 290 """ 291 292 # Components. 293 data.t_0_comp_cubed = data.t_0_comp**3 294 data.t_1_comp_cubed = data.t_1_comp**3 295 data.t_2_comp_cubed = data.t_2_comp**3 296 297 298 # tm-tm partial derivative. 299 ########################### 300 301 # Components. 302 data.inv_d2Diso_dtm2 = 3.0 * diff_data.params[0]**3 303 304 # t0. 305 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_0_comp_cubed 306 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_0_comp_sqrd 307 if a == 0.0 or b == 0.0: 308 data.d2ti[0, 0, 0] = 1e99 309 else: 310 data.d2ti[0, 0, 0] = 1.0 / a + 1.0 / b 311 312 # t1. 313 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_1_comp_cubed 314 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_1_comp_sqrd 315 if a == 0.0 or b == 0.0: 316 data.d2ti[0, 0, 1] = 1e99 317 else: 318 data.d2ti[0, 0, 1] = 1.0 / a + 1.0 / b 319 320 # t2. 321 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_2_comp_cubed 322 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_2_comp_sqrd 323 if a == 0.0 or b == 0.0: 324 data.d2ti[0, 0, 2] = 1e99 325 else: 326 data.d2ti[0, 0, 2] = 1.0 / a + 1.0 / b 327 328 329 # tm-Da partial derivative. 330 ########################### 331 332 # t0. 333 a = -3.0 * data.inv_dDiso_dtm * data.t_0_comp_cubed 334 if a == 0.0: 335 data.d2ti[0, 1, 0] = data.d2ti[1, 0, 0] = 1e99 336 else: 337 data.d2ti[0, 1, 0] = data.d2ti[1, 0, 0] = 1.0 / a 338 339 # t1. 340 a = -6.0 * data.inv_dDiso_dtm * data.t_1_comp_cubed 341 if a == 0.0: 342 data.d2ti[0, 1, 1] = data.d2ti[1, 0, 1] = 1e99 343 else: 344 data.d2ti[0, 1, 1] = data.d2ti[1, 0, 1] = 1.0 / a 345 346 # t2. 347 a = 3.0 * data.inv_dDiso_dtm * data.t_2_comp_cubed 348 if a == 0.0: 349 data.d2ti[0, 1, 2] = data.d2ti[1, 0, 2] = 1e99 350 else: 351 data.d2ti[0, 1, 2] = data.d2ti[1, 0, 2] = 1.0 / a 352 353 354 # Da-Da partial derivative. 355 ########################### 356 357 # t0. 358 a = 3.0 * data.t_0_comp_cubed 359 if a == 0.0: 360 data.d2ti[1, 1, 0] = 1e99 361 else: 362 data.d2ti[1, 1, 0] = 1.0 / a 363 364 # t1. 365 a = 12.0 * data.t_1_comp_cubed 366 if a == 0.0: 367 data.d2ti[1, 1, 1] = 1e99 368 else: 369 data.d2ti[1, 1, 1] = 1.0 / a 370 371 # t2. 372 a = 3.0 * data.t_2_comp_cubed 373 if a == 0.0: 374 data.d2ti[1, 1, 2] = 1e99 375 else: 376 data.d2ti[1, 1, 2] = 1.0 / a

377 378 379 380 381 # Anisotropic global correlation time equation. 382 ############################################### 383

384 -def calc_aniso_ti(data, diff_data):

385 """Diffusional correlation times. 386 387 The equations for the parameters {Diso, Da, Dr} are: 388 389 t-2 = 1/6 (Diso + Da)**-1 390 391 t-1 = 1/6 (Diso - (Da + Dr)/2)**-1 392 393 t0 = 1/6 (Diso - mu)**-1 394 395 t1 = 1/6 (Diso - (Da - Dr)/2)**-1 396 397 t2 = 1/6 (Diso + mu)**-1 398 399 where: 400 __________________ 401 mu = \/ Da**2 + Dr**2 / 3 402 403 The diffusion parameter set in data.diff_params is {tm, Da, Dr, alpha, beta, gamma}. 404 """ 405 406 # Calculate Diso. 407 if diff_data.params[0] == 0.0: 408 data.Diso = 1e99 409 else: 410 data.Diso = 1.0 / (6.0 * diff_data.params[0]) 411 412 # Calculate mu. 413 data.mu = sqrt(diff_data.params[1]**2 + diff_data.params[2]**2 / 3.0) 414 415 # Components. 416 data.t_m2_comp = data.Diso + diff_data.params[1] 417 data.t_m1_comp = data.Diso - 0.5 * (diff_data.params[1] + diff_data.params[2]) 418 data.t_0_comp = data.Diso - data.mu 419 data.t_1_comp = data.Diso - 0.5 * (diff_data.params[1] - diff_data.params[2]) 420 data.t_2_comp = data.Diso + data.mu 421 422 # t-2. 423 data.ti[0] = 6.0 * data.t_m2_comp 424 if data.ti[0] == 0.0: 425 data.ti[0] = 1e99 426 else: 427 data.ti[0] = 1.0 / data.ti[0] 428 429 # t-1. 430 data.ti[1] = 6.0 * data.t_m1_comp 431 if data.ti[1] == 0.0: 432 data.ti[1] = 1e99 433 else: 434 data.ti[1] = 1.0 / data.ti[1] 435 436 # t0. 437 data.ti[2] = 6.0 * data.t_0_comp 438 if data.ti[2] == 0.0: 439 data.ti[2] = 1e99 440 else: 441 data.ti[2] = 1.0 / data.ti[2] 442 443 # t1. 444 data.ti[3] = 6.0 * data.t_1_comp 445 if data.ti[3] == 0.0: 446 data.ti[3] = 1e99 447 else: 448 data.ti[3] = 1.0 / data.ti[3] 449 450 # t2. 451 data.ti[4] = 6.0 * data.t_2_comp 452 if data.ti[4] == 0.0: 453 data.ti[4] = 1e99 454 else: 455 data.ti[4] = 1.0 / data.ti[4]

456 457 458 459 # Anisotropic global correlation time gradient. 460 ############################################### 461

462 -def calc_aniso_dti(data, diff_data):

463 """Diffusional correlation time gradients. 464 465 tm partial derivatives 466 ~~~~~~~~~~~~~~~~~~~~~~ 467 468 dt-2 1 dDiso 469 ---- = - - ----- (Diso + Da)**-2 470 dtm 6 dtm 471 472 dt-1 1 dDiso 473 ---- = - - ----- (Diso - (Da + Dr)/2)**-2 474 dtm 6 dtm 475 476 dt0 1 dDiso 477 --- = - - ----- (Diso - mu)**-2 478 dtm 6 dtm 479 480 dt1 1 dDiso 481 --- = - - ----- (Diso - (Da - Dr)/2)**-2 482 dtm 6 dtm 483 484 dt2 1 dDiso 485 --- = - - ----- (Diso + mu)**-2 486 dtm 6 dtm 487 488 489 dDiso 490 ----- = -1/6 * tm**-2 491 dtm 492 493 494 Da partial derivatives 495 ~~~~~~~~~~~~~~~~~~~~~~ 496 497 dt-2 498 ---- = -1/6 (Diso + Da)**-2 499 dDa 500 501 dt-1 502 ---- = 1/12 (Diso - (Da + Dr)/2)**-2 503 dDa 504 505 dt0 506 --- = 1/6 Da/mu (Diso - mu)**-2 507 dDa 508 509 dt1 510 --- = 1/12 (Diso - (Da - Dr)/2)**-2 511 dDa 512 513 dt2 514 --- = -1/6 Da/mu (Diso + mu)**-2 515 dDa 516 517 518 Dr partial derivatives 519 ~~~~~~~~~~~~~~~~~~~~~~ 520 521 dt-2 522 ---- = 0 523 dDr 524 525 dt-1 526 ---- = 1/12 (Diso - (Da + Dr)/2)**-2 527 dDr 528 529 dt0 530 --- = 1/18 Dr/mu (Diso - mu)**-2 531 dDr 532 533 dt1 534 --- = -1/12 (Diso - (Da - Dr)/2)**-2 535 dDr 536 537 dt2 538 --- = -1/18 Dr/mu (Diso + mu)**-2 539 dDr 540 541 The diffusion parameter set in data.diff_params is {tm, Da, Dr, alpha, beta, gamma}. 542 """ 543 544 # Components. 545 data.t_m2_comp_sqrd = data.t_m2_comp**2 546 data.t_m1_comp_sqrd = data.t_m1_comp**2 547 data.t_0_comp_sqrd = data.t_0_comp**2 548 data.t_1_comp_sqrd = data.t_1_comp**2 549 data.t_2_comp_sqrd = data.t_2_comp**2 550 551 552 # tm partial derivative. 553 ######################## 554 555 # Components. 556 data.inv_dDiso_dtm = -6.0 * diff_data.params[0]**2 557 558 # t-2. 559 data.dti[0, 0] = -6.0 * data.inv_dDiso_dtm * data.t_m2_comp_sqrd 560 if data.dti[0, 0] == 0.0: 561 data.dti[0, 0] = 1e99 562 else: 563 data.dti[0, 0] = 1.0 / data.dti[0, 0] 564 565 # t-1. 566 data.dti[0, 1] = -6.0 * data.inv_dDiso_dtm * data.t_m1_comp_sqrd 567 if data.dti[0, 1] == 0.0: 568 data.dti[0, 1] = 1e99 569 else: 570 data.dti[0, 1] = 1.0 / data.dti[0, 1] 571 572 # t0. 573 data.dti[0, 2] = -6.0 * data.inv_dDiso_dtm * data.t_0_comp_sqrd 574 if data.dti[0, 2] == 0.0: 575 data.dti[0, 2] = 1e99 576 else: 577 data.dti[0, 2] = 1.0 / data.dti[0, 2] 578 579 # t1. 580 data.dti[0, 3] = -6.0 * data.inv_dDiso_dtm * data.t_1_comp_sqrd 581 if data.dti[0, 3] == 0.0: 582 data.dti[0, 3] = 1e99 583 else: 584 data.dti[0, 3] = 1.0 / data.dti[0, 3] 585 586 # t2. 587 data.dti[0, 4] = -6.0 * data.inv_dDiso_dtm * data.t_2_comp_sqrd 588 if data.dti[0, 4] == 0.0: 589 data.dti[0, 4] = 1e99 590 else: 591 data.dti[0, 4] = 1.0 / data.dti[0, 4] 592 593 594 # Da partial derivative. 595 ######################## 596 597 # t-2. 598 data.dti[1, 0] = -6.0 * data.t_m2_comp_sqrd 599 if data.dti[1, 0] == 0.0: 600 data.dti[1, 0] = 1e99 601 else: 602 data.dti[1, 0] = 1.0 / data.dti[1, 0] 603 604 # t-1. 605 data.dti[1, 1] = 12.0 * data.t_m1_comp_sqrd 606 if data.dti[1, 1] == 0.0: 607 data.dti[1, 1] = 1e99 608 else: 609 data.dti[1, 1] = 1.0 / data.dti[1, 1] 610 611 # t0. 612 data.dti[1, 2] = 6.0 * data.mu * data.t_0_comp_sqrd 613 if data.dti[1, 2] == 0.0: 614 if data.mu != 0.0: 615 data.dti[1, 2] = 1e99 616 else: 617 data.dti[1, 2] = 0.0 618 else: 619 data.dti[1, 2] = diff_data.params[1] / data.dti[1, 2] 620 621 # t1. 622 data.dti[1, 3] = 12.0 * data.t_1_comp_sqrd 623 if data.dti[1, 3] == 0.0: 624 data.dti[1, 3] = 1e99 625 else: 626 data.dti[1, 3] = 1.0 / data.dti[1, 3] 627 628 # t2. 629 data.dti[1, 4] = -6.0 * data.mu * data.t_2_comp_sqrd 630 if data.dti[1, 4] == 0.0: 631 if data.mu != 0.0: 632 data.dti[1, 4] = 1e99 633 else: 634 data.dti[1, 4] = 0.0 635 else: 636 data.dti[1, 4] = diff_data.params[1] / data.dti[1, 4] 637 638 639 # Dr partial derivative. 640 ######################## 641 642 # t-1. 643 data.dti[2, 1] = 12.0 * data.t_m1_comp_sqrd 644 if data.dti[2, 1] == 0.0: 645 data.dti[2, 1] = 1e99 646 else: 647 data.dti[2, 1] = 1.0 / data.dti[2, 1] 648 649 # t0. 650 data.dti[2, 2] = 18.0 * data.mu * data.t_0_comp_sqrd 651 if data.dti[2, 2] == 0.0: 652 if data.mu != 0.0: 653 data.dti[2, 2] = 1e99 654 else: 655 data.dti[2, 2] = 0.0 656 else: 657 data.dti[2, 2] = diff_data.params[2] / data.dti[2, 2] 658 659 # t1. 660 data.dti[2, 3] = -12.0 * data.t_1_comp_sqrd 661 if data.dti[2, 3] == 0.0: 662 data.dti[2, 3] = 1e99 663 else: 664 data.dti[2, 3] = 1.0 / data.dti[2, 3] 665 666 # t2. 667 data.dti[2, 4] = -18.0 * data.mu * data.t_2_comp_sqrd 668 if data.dti[2, 4] == 0.0: 669 if data.mu != 0.0: 670 data.dti[2, 4] = 1e99 671 else: 672 data.dti[2, 4] = 0.0 673 else: 674 data.dti[2, 4] = diff_data.params[2] / data.dti[2, 4]

675 676 677 678 # Anisotropic global correlation time Hessians. 679 ############################################### 680

681 -def calc_aniso_d2ti(data, diff_data):

682 """Diffusional correlation time Hessians. 683 684 tm-tm partial derivatives 685 ~~~~~~~~~~~~~~~~~~~~~~~~~ 686 687 d2t-2 1 / dDiso \ 2 1 d2Diso 688 ----- = - | ----- | (Diso + Da)**-3 - - ------ (Diso + Da)**-2 689 dtm2 3 \ dtm / 6 dtm2 690 691 d2t-1 1 / dDiso \ 2 1 d2Diso 692 ----- = - | ----- | (Diso - (Da + Dr)/2)**-3 - - ------ (Diso - (Da + Dr)/2)**-2 693 dtm2 3 \ dtm / 6 dtm2 694 695 d2t0 1 / dDiso \ 2 1 d2Diso 696 ---- = - | ----- | (Diso - mu)**-3 - - ------ (Diso - mu)**-2 697 dtm2 3 \ dtm / 6 dtm2 698 699 d2t1 1 / dDiso \ 2 1 d2Diso 700 ---- = - | ----- | (Diso - (Da - Dr)/2)**-3 - - ------ (Diso - (Da - Dr)/2)**-2 701 dtm2 3 \ dtm / 6 dtm2 702 703 d2t2 1 / dDiso \ 2 1 d2Diso 704 ---- = - | ----- | (Diso + mu)**-3 - - ------ (Diso + mu)**-2 705 dtm2 3 \ dtm / 6 dtm2 706 707 708 d2Diso 709 ------ = 1/3 * tm**-3 710 dtm2 711 712 713 tm-Da partial derivatives 714 ~~~~~~~~~~~~~~~~~~~~~~~~~ 715 716 d2t-2 1 dDiso 717 ------- = - ----- (Diso + Da)**-3 718 dtm.dDa 3 dtm 719 720 d2t-1 1 dDiso 721 ------- = - - ----- (Diso - (Da + Dr)/2)**-3 722 dtm.dDa 6 dtm 723 724 d2t0 1 dDiso Da 725 ------- = - - ----- -- (Diso - mu)**-3 726 dtm.dDa 3 dtm mu 727 728 d2t1 1 dDiso 729 ------- = - - ----- (Diso - (Da - Dr)/2)**-3 730 dtm.dDa 6 dtm 731 732 d2t2 1 dDiso Da 733 ------- = - ----- -- (Diso + mu)**-3 734 dtm.dDa 3 dtm mu 735 736 737 tm-Dr partial derivatives 738 ~~~~~~~~~~~~~~~~~~~~~~~~~ 739 740 d2t-2 741 ------- = 0 742 dtm.dDr 743 744 d2t-1 1 dDiso 745 ------- = - - ----- (Diso - (Da + Dr)/2)**-3 746 dtm.dDr 6 dtm 747 748 d2t0 1 dDiso Dr 749 ------- = - - ----- -- (Diso - mu)**-3 750 dtm.dDr 9 dtm mu 751 752 d2t1 1 dDiso 753 ------- = - ----- (Diso - (Da - Dr)/2)**-3 754 dtm.dDr 6 dtm 755 756 d2t2 1 dDiso Dr 757 ------- = - ----- -- (Diso + mu)**-3 758 dtm.dDr 9 dtm mu 759 760 761 Da-Da partial derivatives 762 ~~~~~~~~~~~~~~~~~~~~~~~~~ 763 764 d2t-2 765 ----- = 1/3 (Diso + Da)**-3 766 dDa2 767 768 d2t-1 769 ----- = 1/12 (Diso - (Da + Dr)/2)**-3 770 dDa2 771 772 d2t0 1 Da**2 1 / Da**2 \ 773 ---- = - ----- (Diso - mu)**-3 - --- | ----- - 1 | (Diso - mu)**-2 774 dDa2 3 mu**2 6mu \ mu**2 / 775 776 d2t1 777 ---- = 1/12 (Diso - (Da - Dr)/2)**-3 778 dDa2 779 780 d2t2 1 Da**2 1 / Da**2 \ 781 ---- = - ----- (Diso + mu)**-3 + --- | ----- - 1 | (Diso + mu)**-2 782 dDa2 3 mu**2 6mu \ mu**2 / 783 784 785 Da-Dr partial derivatives 786 ~~~~~~~~~~~~~~~~~~~~~~~~~ 787 788 d2t-2 789 ------- = 0 790 dDa.dDr 791 792 d2t-1 793 ------- = 1/12 (Diso - (Da + Dr)/2)**-3 794 dDa.dDr 795 796 d2t0 1 Da.Dr 1 Da.Dr 797 ------- = - ----- (Diso - mu)**-3 - -- ----- (Diso - mu)**-2 798 dDa.dDr 9 mu**2 18 mu**3 799 800 d2t1 801 ------- = -1/12 (Diso - (Da - Dr)/2)**-3 802 dDa.dDr 803 804 d2t2 1 Da.Dr 1 Da.Dr 805 ------- = - ----- (Diso + mu)**-3 + -- ----- (Diso + mu)**-2 806 dDa.dDr 9 mu**2 18 mu**3 807 808 809 Dr-Dr partial derivatives 810 ~~~~~~~~~~~~~~~~~~~~~~~~~ 811 812 d2t-2 813 ----- = 0 814 dDr2 815 816 d2t-1 817 ----- = 1/12 (Diso - (Da + Dr)/2)**-3 818 dDr2 819 820 d2t0 1 Dr**2 1 / Dr**2 \ 821 ---- = -- ----- (Diso - mu)**-3 - ---- | ------ - 1 | (Diso - mu)**-2 822 dDr2 27 mu**2 18mu \ 3mu**2 / 823 824 d2t1 825 ---- = 1/12 (Diso - (Da - Dr)/2)**-3 826 dDr2 827 828 d2t2 1 Dr**2 1 / Dr**2 \ 829 ---- = -- ----- (Diso + mu)**-3 + ---- | ------ - 1 | (Diso + mu)**-2 830 dDr2 27 mu**2 18mu \ 3mu**2 / 831 832 833 The diffusion parameter set in data.diff_params is {tm, Da, Dr, alpha, beta, gamma}. 834 """ 835 836 # Components. 837 data.t_m2_comp_cubed = data.t_m2_comp**3 838 data.t_m1_comp_cubed = data.t_m1_comp**3 839 data.t_0_comp_cubed = data.t_0_comp**3 840 data.t_1_comp_cubed = data.t_1_comp**3 841 data.t_2_comp_cubed = data.t_2_comp**3 842 843 844 # tm-tm partial derivative. 845 ########################### 846 847 # Components. 848 data.inv_d2Diso_dtm2 = 3.0 * diff_data.params[0]**3 849 850 # t-2. 851 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_m2_comp_cubed 852 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_m2_comp_sqrd 853 if a == 0.0 or b == 0.0: 854 data.d2ti[0, 0, 0] = 1e99 855 else: 856 data.d2ti[0, 0, 0] = 1.0 / a + 1.0 / b 857 858 # t-1. 859 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_m1_comp_cubed 860 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_m1_comp_sqrd 861 if a == 0.0 or b == 0.0: 862 data.d2ti[0, 0, 1] = 1e99 863 else: 864 data.d2ti[0, 0, 1] = 1.0 / a + 1.0 / b 865 866 # t0. 867 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_0_comp_cubed 868 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_0_comp_sqrd 869 if a == 0.0 or b == 0.0: 870 data.d2ti[0, 0, 2] = 1e99 871 else: 872 data.d2ti[0, 0, 2] = 1.0 / a + 1.0 / b 873 874 # t1. 875 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_1_comp_cubed 876 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_1_comp_sqrd 877 if a == 0.0 or b == 0.0: 878 data.d2ti[0, 0, 3] = 1e99 879 else: 880 data.d2ti[0, 0, 3] = 1.0 / a + 1.0 / b 881 882 # t2. 883 a = 3.0 * data.inv_dDiso_dtm**2 * data.t_2_comp_cubed 884 b = -6.0 * data.inv_d2Diso_dtm2 * data.t_2_comp_sqrd 885 if a == 0.0 or b == 0.0: 886 data.d2ti[0, 0, 4] = 1e99 887 else: 888 data.d2ti[0, 0, 4] = 1.0 / a + 1.0 / b 889 890 891 # tm-Da partial derivative. 892 ########################### 893 894 # t-2. 895 a = 3.0 * data.inv_dDiso_dtm * data.t_m2_comp_cubed 896 if a == 0.0: 897 data.d2ti[0, 1, 0] = data.d2ti[1, 0, 0] = 1e99 898 else: 899 data.d2ti[0, 1, 0] = data.d2ti[1, 0, 0] = 1.0 / a 900 901 # t-1. 902 a = -6.0 * data.inv_dDiso_dtm * data.t_m1_comp_cubed 903 if a == 0.0: 904 data.d2ti[0, 1, 1] = data.d2ti[1, 0, 1] = 1e99 905 else: 906 data.d2ti[0, 1, 1] = data.d2ti[1, 0, 1] = 1.0 / a 907 908 # t0. 909 a = -3.0 * data.inv_dDiso_dtm * data.mu * data.t_0_comp_cubed 910 if a == 0.0: 911 data.d2ti[0, 1, 2] = data.d2ti[1, 0, 2] = 1e99 912 else: 913 data.d2ti[0, 1, 2] = data.d2ti[1, 0, 2] = diff_data.params[1] / a 914 915 # t1. 916 a = -6.0 * data.inv_dDiso_dtm * data.t_1_comp_cubed 917 if a == 0.0: 918 data.d2ti[0, 1, 3] = data.d2ti[1, 0, 3] = 1e99 919 else: 920 data.d2ti[0, 1, 3] = data.d2ti[1, 0, 3] = 1.0 / a 921 922 # t2. 923 a = 3.0 * data.inv_dDiso_dtm * data.mu * data.t_2_comp_cubed 924 if a == 0.0: 925 data.d2ti[0, 1, 4] = data.d2ti[1, 0, 4] = 1e99 926 else: 927 data.d2ti[0, 1, 4] = data.d2ti[1, 0, 4] = diff_data.params[1] / a 928 929 930 # tm-Dr partial derivative. 931 ########################### 932 933 # t-1. 934 a = -6.0 * data.inv_dDiso_dtm * data.t_m1_comp_cubed 935 if a == 0.0: 936 data.d2ti[0, 2, 1] = data.d2ti[2, 0, 1] = 1e99 937 else: 938 data.d2ti[0, 2, 1] = data.d2ti[2, 0, 1] = 1.0 / a 939 940 # t0. 941 a = -9.0 * data.inv_dDiso_dtm * data.mu * data.t_0_comp_cubed 942 if a == 0.0: 943 data.d2ti[0, 2, 2] = data.d2ti[2, 0, 2] = 1e99 944 else: 945 data.d2ti[0, 2, 2] = data.d2ti[2, 0, 2] = diff_data.params[2] / a 946 947 # t1. 948 a = 6.0 * data.inv_dDiso_dtm * data.t_1_comp_cubed 949 if a == 0.0: 950 data.d2ti[0, 2, 3] = data.d2ti[2, 0, 3] = 1e99 951 else: 952 data.d2ti[0, 2, 3] = data.d2ti[2, 0, 3] = 1.0 / a 953 954 # t2. 955 a = 9.0 * data.inv_dDiso_dtm * data.mu * data.t_2_comp_cubed 956 if a == 0.0: 957 data.d2ti[0, 2, 4] = data.d2ti[2, 0, 4] = 1e99 958 else: 959 data.d2ti[0, 2, 4] = data.d2ti[2, 0, 4] = diff_data.params[2] / a 960 961 962 # Da-Da partial derivative. 963 ########################### 964 965 # t-2. 966 a = 3.0 * data.t_m2_comp_cubed 967 if a == 0.0: 968 data.d2ti[1, 1, 0] = 1e99 969 else: 970 data.d2ti[1, 1, 0] = 1.0 / a 971 972 # t-1. 973 a = 12.0 * data.t_m1_comp_cubed 974 if a == 0.0: 975 data.d2ti[1, 1, 1] = 1e99 976 else: 977 data.d2ti[1, 1, 1] = 1.0 / a 978 979 # t0. 980 a = 3.0 * data.mu**2 * data.t_0_comp_cubed 981 b = -6.0 * data.mu * data.t_0_comp_sqrd 982 if a == 0.0 or b == 0.0: 983 data.d2ti[1, 1, 2] = 1e99 984 else: 985 data.d2ti[1, 1, 2] = diff_data.params[1]**2 / a + (diff_data.params[1]**2 / data.mu**2 - 1.0) / b 986 987 # t1. 988 a = 12.0 * data.t_1_comp_cubed 989 if a == 0.0: 990 data.d2ti[1, 1, 3] = 1e99 991 else: 992 data.d2ti[1, 1, 3] = 1.0 / a 993 994 # t2. 995 a = 3.0 * data.mu**2 * data.t_2_comp_cubed 996 b = 6.0 * data.mu * data.t_2_comp_sqrd 997 if a == 0.0 or b == 0.0: 998 data.d2ti[1, 1, 4] = 1e99 999 else: 1000 data.d2ti[1, 1, 4] = diff_data.params[1]**2 / a + (diff_data.params[1]**2 / data.mu**2 - 1.0) / b 1001 1002 1003 # Da-Dr partial derivative. 1004 ########################### 1005 1006 # t-1. 1007 a = 12.0 * data.t_m1_comp_cubed 1008 if a == 0.0: 1009 data.d2ti[1, 2, 1] = data.d2ti[2, 1, 1] = 1e99 1010 else: 1011 data.d2ti[1, 2, 1] = data.d2ti[2, 1, 1] = 1.0 / a 1012 1013 # t0. 1014 a = 9.0 * data.mu**2 * data.t_0_comp_cubed 1015 b = -18.0 * data.mu**3 * data.t_0_comp_sqrd 1016 if a == 0.0 or b == 0.0: 1017 data.d2ti[1, 2, 2] = data.d2ti[2, 1, 2] = 1e99 1018 else: 1019 data.d2ti[1, 2, 2] = data.d2ti[2, 1, 2] = diff_data.params[1] * diff_data.params[2] / a + diff_data.params[1] * diff_data.params[2] / b 1020 1021 # t1. 1022 a = -12.0 * data.t_1_comp_cubed 1023 if a == 0.0: 1024 data.d2ti[1, 2, 3] = data.d2ti[2, 1, 3] = 1e99 1025 else: 1026 data.d2ti[1, 2, 3] = data.d2ti[2, 1, 3] = 1.0 / a 1027 1028 # t2. 1029 a = 9.0 * data.mu**2 * data.t_2_comp_cubed 1030 b = 18.0 * data.mu**3 * data.t_2_comp_sqrd 1031 if a == 0.0 or b == 0.0: 1032 data.d2ti[1, 2, 4] = data.d2ti[2, 1, 4] = 1e99 1033 else: 1034 data.d2ti[1, 2, 4] = data.d2ti[2, 1, 4] = diff_data.params[1] * diff_data.params[2] / a + diff_data.params[1] * diff_data.params[2] / b 1035 1036 1037 # Dr-Dr partial derivative. 1038 ########################### 1039 1040 # t-1. 1041 a = 12.0 * data.t_m1_comp_cubed 1042 if a == 0.0: 1043 data.d2ti[1, 1, 1] = 1e99 1044 else: 1045 data.d2ti[1, 1, 1] = 1.0 / a 1046 1047 # t0. 1048 a = 27.0 * data.mu**2 * data.t_0_comp_cubed 1049 b = -18.0 * data.mu * data.t_0_comp_sqrd 1050 if a == 0.0 or b == 0.0: 1051 data.d2ti[1, 1, 2] = 1e99 1052 else: 1053 data.d2ti[1, 1, 2] = diff_data.params[2]**2 / a + (diff_data.params[2]**2 / (3.0 * data.mu**2) - 1.0) / b 1054 1055 # t1. 1056 a = 12.0 * data.t_1_comp_cubed 1057 if a == 0.0: 1058 data.d2ti[1, 1, 3] = 1e99 1059 else: 1060 data.d2ti[1, 1, 3] = 1.0 / a 1061 1062 # t2. 1063 a = 27.0 * data.mu**2 * data.t_2_comp_cubed 1064 b = 18.0 * data.mu * data.t_2_comp_sqrd 1065 if a == 0.0 or b == 0.0: 1066 data.d2ti[1, 1, 4] = 1e99 1067 else: 1068 data.d2ti[1, 1, 4] = diff_data.params[2]**2 / a + (diff_data.params[2]**2 / (3.0 * data.mu**2) - 1.0) / b

1069

Source Code for Module maths_fns.correlation_time