maths_fns.weights

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 from Numeric import outerproduct 25 26 27 # Isotropic weight equation. 28 ############################ 29

30 -def calc_iso_ci(data, diff_data):

31 """Weight equations for isotropic diffusion. 32 33 c0 = 1 34 """ 35 36 data.ci[0] = 1.0

37 38 39 40 # Axially symmetric weight equation. 41 #################################### 42

43 -def calc_axial_ci(data, diff_data):

44 """Weight equations for axially symmetric diffusion. 45 46 The equations are: 47 48 c0 = 1/4 (3delta**2 - 1)**2 49 c1 = 3delta**2 (1 - delta**2) 50 c2 = 3/4 (1 - delta**2)**2 51 52 where delta is the dot product of the unit bond vector and the unit vector along Dpar. 53 """ 54 55 data.ci[0] = 0.25 * (3.0 * data.delta**2 - 1.0)**2 56 data.ci[1] = 3.0 * data.delta**2 * (1.0 - data.delta**2) 57 data.ci[2] = 0.75 * (1.0 - data.delta**2)**2

58 59 60 61 # Axially symmetric weight gradient. 62 #################################### 63

64 -def calc_axial_dci(data, diff_data):

65 """Weight gradient for axially symmetric diffusion. 66 67 The equations are: 68 69 dc0 ddelta 70 ----- = 3delta (3delta**2 - 1) ------ 71 dpsii dpsii 72 73 dc1 ddelta 74 ----- = 6delta (1 - 2delta**2) ------ 75 dpsii dpsii 76 77 dc2 ddelta 78 ----- = 3delta (delta**2 - 1) ------ 79 dpsii dpsii 80 81 where psi = {theta, phi} 82 """ 83 84 data.dci[2:, 0] = 3.0 * data.delta * (3.0 * data.delta**2 - 1.0) * data.ddelta_dpsi 85 data.dci[2:, 1] = 6.0 * data.delta * (1.0 - 2.0 * data.delta**2) * data.ddelta_dpsi 86 data.dci[2:, 2] = 3.0 * data.delta * (data.delta**2 - 1.0) * data.ddelta_dpsi

87 88 89 90 # Axially symmetric weight Hessian. 91 ################################### 92

93 -def calc_axial_d2ci(data, diff_data):

94 """Weight Hessian for axially symmetric diffusion. 95 96 The equations are: 97 98 d2c0 / ddelta ddelta d2delta \ 99 ----------- = 3 |(9delta**2 - 1) ------ . ------ + delta (3delta**2 - 1) ----------- | 100 dpsii.dpsij \ dpsii dpsij dpsii.dpsij / 101 102 d2c1 / ddelta ddelta d2delta \ 103 ----------- = 6 |(1 - 6delta**2) ------ . ------ + delta (1 - 2delta**2) ----------- | 104 dpsii.dpsij \ dpsii dpsij dpsii.dpsij / 105 106 d2c2 / ddelta ddelta d2delta \ 107 ----------- = 3 |(3delta**2 - 1) ------ . ------ + delta (delta**2 - 1) ----------- | 108 dpsii.dpsij \ dpsii dpsij dpsii.dpsij / 109 110 where psi = {theta, phi} 111 """ 112 113 # Outer product. 114 op = outerproduct(data.ddelta_dpsi, data.ddelta_dpsi) 115 116 # Hessian. 117 data.d2ci[2:, 2:, 0] = 3.0 * ((9.0 * data.delta**2 - 1.0) * op + data.delta * (3.0 * data.delta**2 - 1.0) * data.d2delta_dpsi2) 118 data.d2ci[2:, 2:, 1] = 6.0 * ((1.0 - 6.0 * data.delta**2) * op + data.delta * (1.0 - 2.0 * data.delta**2) * data.d2delta_dpsi2) 119 data.d2ci[2:, 2:, 2] = 3.0 * ((3.0 * data.delta**2 - 1.0) * op + data.delta * (data.delta**2 - 1.0) * data.d2delta_dpsi2)

120 121 122 123 # Anisotropic weight equation. 124 ############################## 125

126 -def calc_aniso_ci(data, diff_data):

127 """Weight equations for anisotropic diffusion. 128 129 In the following equations, the following short-hand notation will be used: 130 131 da = delta_alpha 132 133 db = delta_beta 134 135 dg = delta_gamma 136 137 138 The equations are: 139 140 c-2 = 3da**2.db**2 141 142 143 c-1 = 3da**2.dg**2 144 145 146 c0 = 1/4 (3(da**4 + db**4 + dg**4) - 1 - e) 147 148 149 c1 = 3db**2.dg**2 150 151 152 c2 = 1/4 (3(da**4 + db**4 + dg**4) - 1 + e) 153 154 155 Da - Dr Da + Dr 2Da 156 e = ------- (da**4 + 2db**2.dg**2) + ------- (db**4 + 2da**2.dg**2) - --- (dg**4 + 2da**2.db**2) 157 mu mu mu 158 159 160 where: 161 __________________ 162 mu = V Da**2 + Dr**2 / 3 163 164 delta_alpha is the dot product of the unit bond vector and the unit vector along Dx. 165 166 delta_beta is the dot product of the unit bond vector and the unit vector along Dy. 167 168 delta_gamma is the dot product of the unit bond vector and the unit vector along Dz. 169 170 alpha (in delta_alpha) is the directional cosine along Dx. 171 172 beta (in delta_beta) is the directional cosine along Dy. 173 174 gamma (in delta_gamma) is the directional cosine along Dz. 175 """ 176 177 # Calculate mu. 178 data.mu = sqrt(diff_data.params[1]**2 + diff_data.params[2]**2 / 3.0) 179 180 # Components. 181 data.c_alpha = data.delta_alpha**4 + 2.0 * data.delta_beta**2 * data.delta_gamma**2 182 data.c_beta = data.delta_beta**4 + 2.0 * data.delta_alpha**2 * data.delta_gamma**2 183 data.c_gamma = data.delta_gamma**4 + 2.0 * data.delta_alpha**2 * data.delta_beta**2 184 185 if data.mu == 0.0: 186 data.e1 = 0.0 187 data.e2 = 0.0 188 data.e3 = 0.0 189 else: 190 data.e1 = (diff_data.params[1] - diff_data.params[2]) / data.mu 191 data.e2 = (diff_data.params[1] + diff_data.params[2]) / data.mu 192 data.e3 = 2.0 * diff_data.params[1] / data.mu 193 194 # Calculate d. 195 d = 3.0 * (data.delta_alpha**4 + data.delta_beta**4 + data.delta_gamma**4) - 1.0 196 197 # Calculate e. 198 e = data.e1 * data.c_alpha + data.e2 * data.c_beta - data.e3 * data.c_gamma 199 200 # Weight c-2. 201 data.ci[0] = 3.0 * data.delta_alpha**2 * data.delta_beta**2 202 203 # Weight c-1. 204 data.ci[1] = 3.0 * data.delta_alpha**2 * data.delta_gamma**2 205 206 # Weight c0. 207 data.ci[2] = 0.25 * (d - e) 208 209 # Weight c1. 210 data.ci[3] = 3.0 * data.delta_beta**2 * data.delta_gamma**2 211 212 # Weight c2. 213 data.ci[4] = 0.25 * (d + e)

214 215 216 # Anisotropic weight gradient. 217 ############################## 218

219 -def calc_aniso_dci(data, diff_data):

220 """Weight gradient for anisotropic diffusion. 221 222 psii partial derivatives 223 ~~~~~~~~~~~~~~~~~~~~~~~~ 224 225 dc-2 / dda ddb \ 226 ----- = 6da.db | db ----- + da ----- | 227 dpsii \ dpsii dpsii / 228 229 230 dc-1 / dda ddg \ 231 ----- = 6da.dg | dg ----- + da ----- | 232 dpsii \ dpsii dpsii / 233 234 235 dc0 / dda ddb ddg \ de 236 ----- = 3 | da**3 ----- + db**3 ----- + dg**3 ----- | - ----- 237 dpsii \ dpsii dpsii dpsii / dpsii 238 239 240 dc1 / ddb ddg \ 241 ----- = 6db.dg | dg ----- + db ----- | 242 dpsii \ dpsii dpsii / 243 244 245 dc2 / dda ddb ddg \ de 246 ----- = 3 | da**3 ----- + db**3 ----- + dg**3 ----- | + ----- 247 dpsii \ dpsii dpsii dpsii / dpsii 248 249 250 251 de Da - Dr / dda / ddb ddg \ \ 252 ----- = ------- | da**3 ----- + db.dg | dg ----- + db ----- | | 253 dpsii mu \ dpsii \ dpsii dpsii / / 254 255 Da + Dr / ddb / dda ddg \ \ 256 + ------- | db**3 ----- + da.dg | dg ----- + da ----- | | 257 mu \ dpsii \ dpsii dpsii / / 258 259 2Da / ddg / dda ddb \ \ 260 - --- | dg**3 ----- + da.db | db ----- + da ----- | | 261 mu \ dpsii \ dpsii dpsii / / 262 263 264 where psi = {alpha, beta, gamma}. 265 266 267 Da partial derivatives 268 ~~~~~~~~~~~~~~~~~~~~~~ 269 270 dc-2 271 ---- = 0 272 dDa 273 274 dc-1 275 ---- = 0 276 dDa 277 278 dc0 1 de 279 --- = - - --- 280 dDa 4 dDa 281 282 dc1 283 --- = 0 284 dDa 285 286 dc2 1 de 287 --- = - --- 288 dDa 4 dDa 289 290 291 de 1 / (3Da + Dr)Dr (3Da - Dr)Dr 2Dr**2 \ 292 --- = - | ------------ (da**4 + 2db**2.dg**2) - ------------ (db**4 + 2da**2.dg**2) - ------ (dg**4 + 2da**2.db**2) | 293 dDa 3 \ mu**3 mu**3 mu**3 / 294 295 296 297 Dr partial derivatives 298 ~~~~~~~~~~~~~~~~~~~~~~ 299 300 dc-2 301 ---- = 0 302 dDr 303 304 dc-1 305 ---- = 0 306 dDr 307 308 dc0 1 de 309 --- = - - --- 310 dDr 4 dDr 311 312 dc1 313 --- = 0 314 dDr 315 316 dc2 1 de 317 --- = - --- 318 dDr 4 dDr 319 320 321 de 1 / (3Da + Dr)Da (3Da - Dr)Da 2Da.Dr \ 322 --- = - - | ------------ (da**4 + 2db**2.dg**2) - ------------ (db**4 + 2da**2.dg**2) - ------ (dg**4 + 2da**2.db**2) | 323 dDr 3 \ mu**3 mu**3 mu**3 / 324 325 326 tm partial derivatives 327 ~~~~~~~~~~~~~~~~~~~~~~ 328 329 dc-2 330 ---- = 0 331 dtm 332 333 dc-1 334 ---- = 0 335 dtm 336 337 dc0 338 --- = 0 339 dtm 340 341 dc1 342 --- = 0 343 dtm 344 345 dc2 346 --- = 0 347 dtm 348 349 """ 350 351 # psii partial derivative. 352 ########################## 353 354 # Components. 355 data.dc_dpsii_alpha_beta = data.delta_alpha * data.delta_beta * (data.delta_beta * data.ddelta_alpha_dpsi + data.delta_alpha * data.ddelta_beta_dpsi) 356 data.dc_dpsii_alpha_gamma = data.delta_alpha * data.delta_gamma * (data.delta_gamma * data.ddelta_alpha_dpsi + data.delta_alpha * data.ddelta_gamma_dpsi) 357 data.dc_dpsii_beta_gamma = data.delta_beta * data.delta_gamma * (data.delta_gamma * data.ddelta_beta_dpsi + data.delta_beta * data.ddelta_gamma_dpsi) 358 359 data.dc_dpsii_alpha = data.delta_alpha**3 * data.ddelta_alpha_dpsi 360 data.dc_dpsii_beta = data.delta_beta**3 * data.ddelta_beta_dpsi 361 data.dc_dpsii_gamma = data.delta_gamma**3 * data.ddelta_gamma_dpsi 362 363 data.dc_dpsii_alpha_ext = data.dc_dpsii_alpha + data.dc_dpsii_beta_gamma 364 data.dc_dpsii_beta_ext = data.dc_dpsii_beta + data.dc_dpsii_alpha_gamma 365 data.dc_dpsii_gamma_ext = data.dc_dpsii_gamma + data.dc_dpsii_alpha_beta 366 367 # Calculate dd_dpsii. 368 dd_dpsii = 3.0 * (data.dc_dpsii_alpha + data.dc_dpsii_beta + data.dc_dpsii_gamma) 369 370 # Calculate de_dpsii. 371 de_dpsii = data.e1 * data.dc_dpsii_alpha_ext + data.e2 * data.dc_dpsii_beta_ext - data.e3 * data.dc_dpsii_gamma_ext 372 373 # Weight c-2. 374 data.dci[3:, 0] = 6.0 * data.dc_dpsii_alpha_beta 375 376 # Weight c-1. 377 data.dci[3:, 1] = 6.0 * data.dc_dpsii_alpha_gamma 378 379 # Weight c0. 380 data.dci[3:, 2] = dd_dpsii - de_dpsii 381 382 # Weight c1. 383 data.dci[3:, 3] = 6.0 * data.dc_dpsii_beta_gamma 384 385 # Weight c2. 386 data.dci[3:, 4] = dd_dpsii + de_dpsii 387 388 389 # Da partial derivative. 390 ######################## 391 392 # Components. 393 data.mu_cubed = data.mu**3 394 395 if data.mu == 0.0: 396 data.de1_dDa = 0.0 397 data.de2_dDa = 0.0 398 data.de3_dDa = 0.0 399 else: 400 data.de1_dDa = (3.0 * diff_data.params[1] + diff_data.params[2]) * diff_data.params[2] / data.mu_cubed 401 data.de2_dDa = (3.0 * diff_data.params[1] - diff_data.params[2]) * diff_data.params[2] / data.mu_cubed 402 data.de3_dDa = 2.0 * diff_data.params[2]**2 / data.mu_cubed 403 404 # Calculate de_dDa. 405 de_dDa = (data.de1_dDa * data.c_alpha - data.de2_dDa * data.c_beta - data.de3_dDa * data.c_gamma) / 3.0 406 407 # Weight c0. 408 data.dci[1, 2] = -0.25 * de_dDa 409 410 # Weight c2. 411 data.dci[1, 4] = 0.25 * de_dDa 412 413 414 # Dr partial derivative. 415 ######################## 416 417 # Components. 418 if data.mu == 0.0: 419 data.de1_dDr = 0.0 420 data.de2_dDr = 0.0 421 data.de3_dDr = 0.0 422 else: 423 data.de1_dDr = (3.0 * diff_data.params[1] + diff_data.params[2]) * diff_data.params[1] / data.mu_cubed 424 data.de2_dDr = (3.0 * diff_data.params[1] - diff_data.params[2]) * diff_data.params[1] / data.mu_cubed 425 data.de3_dDr = 2.0 * diff_data.params[1] * diff_data.params[2] / data.mu_cubed 426 427 # Calculate de_dDr. 428 de_dDr = - (data.de1_dDr * data.c_alpha - data.de2_dDr * data.c_beta - data.de3_dDr * data.c_gamma) / 3.0 429 430 # Weight c0. 431 data.dci[2, 2] = -0.25 * de_dDr 432 433 # Weight c2. 434 data.dci[2, 4] = 0.25 * de_dDr

435 436 437 # Anisotropic weight Hessian. 438 ############################# 439

440 -def calc_aniso_d2ci(data, diff_data):

441 """Weight Hessian for anisotropic diffusion. 442 443 psii-psij partial derivatives 444 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 445 446 d2c-2 / / dda dda d2da \ / dda ddb ddb dda \ / ddb ddb d2db \ \ 447 ----------- = 6 | db**2 | ----- . ----- + da ----------- | + 2da.db | ----- . ----- + ----- . ----- | + da**2 | ----- . ----- + db ----------- | | 448 dpsii.dpsij \ \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 449 450 451 d2c-1 / / dda dda d2da \ / dda ddg ddg dda \ / ddg ddg d2dg \ \ 452 ----------- = 6 | dg**2 | ----- . ----- + da ----------- | + 2da.dg | ----- . ----- + ----- . ----- | + da**2 | ----- . ----- + dg ----------- | | 453 dpsii.dpsij \ \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 454 455 456 d2c0 / / dda dda d2da \ / ddb ddb d2db \ / ddg ddg d2dg \ \ d2e 457 ----------- = 3 | da**2 | ----- . ----- + da ----------- | + db**2 | ----- . ----- + db ----------- | + dg**2 | ----- . ----- + dg ----------- | | - ----------- 458 dpsii.dpsij \ \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii.dpsij / / dpsii.dpsij 459 460 461 d2c1 / / ddb ddb d2db \ / ddb ddg ddg ddb \ / ddg ddg d2dg \ \ 462 ----------- = 6 | dg**2 | ----- . ----- + db ----------- | + 2db.dg | ----- . ----- + ----- . ----- | + db**2 | ----- . ----- + dg ----------- | | 463 dpsii.dpsij \ \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 464 465 466 d2c2 / / dda dda d2da \ / ddb ddb d2db \ / ddg ddg d2dg \ \ d2e 467 ----------- = 3 | da**2 | ----- . ----- + da ----------- | + db**2 | ----- . ----- + db ----------- | + dg**2 | ----- . ----- + dg ----------- | | + ----------- 468 dpsii.dpsij \ \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii.dpsij / / dpsii.dpsij 469 470 471 472 d2e Da - Dr / / dda dda d2da \ 473 ----------- = ------- | da**2 | 3 ----- . ----- + da ----------- | 474 dpsii.dpsij mu \ \ dpsii dpsij dpsii.dpsij / 475 476 / ddb ddb d2db \ / ddb ddg ddg ddb \ / ddg ddg d2dg \ \ 477 + dg**2 | ----- . ----- + db ----------- | + 2db.dg | ----- . ----- + ----- . ----- | + db**2 | ----- . ----- + dg ----------- | | 478 \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 479 480 Da + Dr / / ddb ddb d2db \ 481 + ------- | db**2 | 3 ----- . ----- + db ----------- | 482 mu \ \ dpsii dpsij dpsii.dpsij / 483 484 / dda dda d2da \ / dda ddg ddg dda \ / ddg ddg d2dg \ \ 485 + dg**2 | ----- . ----- + da ----------- | + 2da.dg | ----- . ----- + ----- . ----- | + da**2 | ----- . ----- + dg ----------- | | 486 \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 487 488 2Da / / ddg ddg d2dg \ 489 - --- | dg**2 | 3 ----- . ----- + da ----------- | 490 mu \ \ dpsii dpsij dpsii.dpsij / 491 492 / dda dda d2da \ / dda ddb ddb dda \ / ddb ddb d2db \ \ 493 + db**2 | ----- . ----- + da ----------- | + 2da.db | ----- . ----- + ----- . ----- | + da**2 | ----- . ----- + db ----------- | | 494 \ dpsii dpsij dpsii.dpsij / \ dpsii dpsij dpsii dpsij / \ dpsii dpsij dpsii.dpsij / / 495 496 497 psii-Da partial derivatives 498 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 499 500 d2c-2 501 --------- = 0 502 dpsii.dDa 503 504 505 d2c-1 506 --------- = 0 507 dpsii.dDa 508 509 510 d2c0 d2e 511 --------- = - --------- 512 dpsii.dDa dpsii.dDa 513 514 515 d2c1 516 --------- = 0 517 dpsii.dDa 518 519 520 d2c2 d2e 521 --------- = --------- 522 dpsii.dDa dpsii.dDa 523 524 525 526 d2e 1 (3Da + Dr)Dr / dda / ddb ddg \ \ 527 --------- = - ------------ | da**3 ----- + db.dg | dg ----- + db ----- | | 528 dpsii.dDa 3 mu**3 \ dpsii \ dpsii dpsii / / 529 530 1 (3Da - Dr)Dr / ddb / dda ddg \ \ 531 - - ------------ | db**3 ----- + da.dg | dg ----- + da ----- | | 532 3 mu**3 \ dpsii \ dpsii dpsii / / 533 534 2 Dr**2 / ddg / dda ddb \ \ 535 - - ----- | dg**3 ----- + da.db | db ----- + da ----- | | 536 3 mu**3 \ dpsii \ dpsii dpsii / / 537 538 539 psii-Dr partial derivatives 540 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 541 542 d2c-2 543 --------- = 0 544 dpsii.dDr 545 546 547 d2c-1 548 --------- = 0 549 dpsii.dDr 550 551 552 d2c0 d2e 553 --------- = - --------- 554 dpsii.dDr dpsii.dDr 555 556 557 d2c1 558 --------- = 0 559 dpsii.dDr 560 561 562 d2c2 d2e 563 --------- = --------- 564 dpsii.dDr dpsii.dDr 565 566 567 568 d2e 1 (3Da + Dr)Da / dda / ddb ddg \ \ 569 --------- = - - ------------ | da**3 ----- + db.dg | dg ----- + db ----- | | 570 dpsii.dDr 3 mu**3 \ dpsii \ dpsii dpsii / / 571 572 1 (3Da - Dr)Da / ddb / dda ddg \ \ 573 - - ------------ | db**3 ----- + da.dg | dg ----- + da ----- | | 574 3 mu**3 \ dpsii \ dpsii dpsii / / 575 576 2 Dr**2 / ddg / dda ddb \ \ 577 - - ----- | dg**3 ----- + da.db | db ----- + da ----- | | 578 3 mu**3 \ dpsii \ dpsii dpsii / / 579 580 581 psii-tm partial derivatives 582 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 583 584 d2c-2 585 --------- = 0 586 dpsii.dtm 587 588 589 d2c-1 590 --------- = 0 591 dpsii.dtm 592 593 594 d2c0 595 --------- = 0 596 dpsii.dtm 597 598 599 d2c1 600 --------- = 0 601 dpsii.dtm 602 603 604 d2c2 605 --------- = 0 606 dpsii.dtm 607 608 609 Da-Da partial derivatives 610 ~~~~~~~~~~~~~~~~~~~~~~~~~ 611 612 d2c-2 613 ------ = 0 614 dDa**2 615 616 617 d2c-1 618 ------ = 0 619 dDa**2 620 621 622 d2c0 1 d2e 623 ------ = - - ------ 624 dDa**2 4 dDa**2 625 626 627 d2c1 628 ------ = 0 629 dDa**2 630 631 632 d2c2 1 d2e 633 ------ = - ------ 634 dDa**2 4 dDa**2 635 636 637 638 d2e 1 / (6Da**2 + 3Da.Dr - Dr**2)Dr (6Da**2 - 3Da.Dr - Dr**2)Dr 6Da.Dr**2 \ 639 ------ = - - | --------------------------- (da**4 + 2db**2.dg**2) - --------------------------- (db**4 + 2da**2.dg**2) - --------- (dg**4 + 2da**2.db**2) | 640 dDa**2 3 \ mu**5 mu**5 mu**5 / 641 642 643 Da-Dr partial derivatives 644 ~~~~~~~~~~~~~~~~~~~~~~~~~ 645 646 d2c-2 647 ------- = 0 648 dDa.dDr 649 650 651 d2c-1 652 ------- = 0 653 dDa.dDr 654 655 656 d2c0 1 d2e 657 ------- = - - ------- 658 dDa.dDr 4 dDa.dDr 659 660 661 d2c1 662 ------- = 0 663 dDa.dDr 664 665 666 d2c2 1 d2e 667 ------- = - ------- 668 dDa.dDr 4 dDa.dDr 669 670 671 672 d2e 1 / 9Da**3 + 6Da**2.Dr - 6Da.Dr**2 - Dr**3 9Da**3 - 6Da**2.Dr - 6Da.Dr**2 + Dr**3 673 ------- = - | -------------------------------------- (da**4 + 2db**2.dg**2) - -------------------------------------- (db**4 + 2da**2.dg**2) 674 dDa.dDr 9 \ mu**5 mu**5 675 676 2(Da**2 - Dr**2)Dr \ 677 - ------------------ (dg**4 + 2da**2.db**2) | 678 mu**5 / 679 680 681 Da-tm partial derivatives 682 ~~~~~~~~~~~~~~~~~~~~~~~~~ 683 684 d2c-2 685 ------- = 0 686 dDa.dtm 687 688 689 d2c-1 690 ------- = 0 691 dDa.dtm 692 693 694 d2c0 695 ------- = 0 696 dDa.dtm 697 698 699 d2c1 700 ------- = 0 701 dDa.dtm 702 703 704 d2c2 705 ------- = 0 706 dDa.dtm 707 708 709 710 Dr-Dr partial derivatives 711 ~~~~~~~~~~~~~~~~~~~~~~~~~ 712 713 d2c-2 714 ------ = 0 715 dDr**2 716 717 718 d2c-1 719 ------ = 0 720 dDr**2 721 722 723 d2c0 1 d2e 724 ------ = - - ------ 725 dDr**2 4 dDr**2 726 727 728 d2c1 729 ------ = 0 730 dDr**2 731 732 733 d2c2 1 d2e 734 ------ = - ------ 735 dDr**2 4 dDr**2 736 737 738 739 d2e 1 / (3Da**2 - 9Da.Dr - 2Dr**2)Da (3Da**2 + 9Da.Dr - 2Dr**2)Da 2(3Da**2 - 2Dr**2)Da \ 740 ------ = - - | ---------------------------- (da**4 + 2db**2.dg**2) + ---------------------------- (db**4 + 2da**2.dg**2) - -------------------- (dg**4 + 2da**2.db**2) | 741 dDr**2 9 \ mu**5 mu**5 mu**5 / 742 743 744 Dr-tm partial derivatives 745 ~~~~~~~~~~~~~~~~~~~~~~~~~ 746 747 d2c-2 748 ------- = 0 749 dDr.dtm 750 751 752 d2c-1 753 ------- = 0 754 dDr.dtm 755 756 757 d2c0 758 ------- = 0 759 dDr.dtm 760 761 762 d2c1 763 ------- = 0 764 dDr.dtm 765 766 767 d2c2 768 ------- = 0 769 dDr.dtm 770 771 """ 772 773 # psii-psij partial derivative. 774 ############################### 775 776 # Outer products. 777 op_alpha = outerproduct(data.ddelta_alpha_dpsi, data.ddelta_alpha_dpsi) 778 op_beta = outerproduct(data.ddelta_beta_dpsi, data.ddelta_beta_dpsi) 779 op_gamma = outerproduct(data.ddelta_gamma_dpsi, data.ddelta_gamma_dpsi) 780 781 op_alpha_beta = outerproduct(data.ddelta_alpha_dpsi, data.ddelta_beta_dpsi) 782 op_beta_alpha = outerproduct(data.ddelta_beta_dpsi, data.ddelta_alpha_dpsi) 783 op_alpha_gamma = outerproduct(data.ddelta_alpha_dpsi, data.ddelta_gamma_dpsi) 784 op_gamma_alpha = outerproduct(data.ddelta_gamma_dpsi, data.ddelta_alpha_dpsi) 785 op_beta_gamma = outerproduct(data.ddelta_beta_dpsi, data.ddelta_gamma_dpsi) 786 op_gamma_beta = outerproduct(data.ddelta_gamma_dpsi, data.ddelta_beta_dpsi) 787 788 # Components. 789 alpha_comp = op_alpha + data.delta_alpha * data.d2delta_alpha_dpsi2 790 beta_comp = op_beta + data.delta_beta * data.d2delta_beta_dpsi2 791 gamma_comp = op_gamma + data.delta_gamma * data.d2delta_gamma_dpsi2 792 793 alpha3_comp = data.delta_alpha**2 * (3.0 * op_alpha + data.delta_alpha * data.d2delta_alpha_dpsi2) 794 beta3_comp = data.delta_beta**2 * (3.0 * op_beta + data.delta_beta * data.d2delta_beta_dpsi2) 795 gamma3_comp = data.delta_gamma**2 * (3.0 * op_gamma + data.delta_gamma * data.d2delta_gamma_dpsi2) 796 797 alpha_beta_comp = data.delta_beta**2 * alpha_comp + 2.0 * data.delta_alpha * data.delta_beta * (op_alpha_beta + op_beta_alpha) + data.delta_alpha**2 * beta_comp 798 alpha_gamma_comp = data.delta_gamma**2 * alpha_comp + 2.0 * data.delta_alpha * data.delta_gamma * (op_alpha_gamma + op_gamma_alpha) + data.delta_alpha**2 * gamma_comp 799 beta_gamma_comp = data.delta_gamma**2 * beta_comp + 2.0 * data.delta_beta * data.delta_gamma * (op_beta_gamma + op_gamma_beta) + data.delta_beta**2 * gamma_comp 800 801 # Calculate d2d_dpsii_dpsij. 802 d2d_dpsii_dpsij = 3.0 * (alpha3_comp + beta3_comp + gamma3_comp) 803 804 # Calculate d2e_dpsii_dpsij. 805 d2e_dpsii_dpsij = data.e1 * (alpha3_comp + alpha_beta_comp) + data.e2 * (beta3_comp + alpha_gamma_comp) - data.e3 * (alpha3_comp + beta_gamma_comp) 806 807 # Weight c-2. 808 data.d2ci[3:, 3:, 0] = 6.0 * alpha_beta_comp 809 810 # Weight c-1. 811 data.d2ci[3:, 3:, 1] = 6.0 * alpha_gamma_comp 812 813 # Weight c0. 814 data.d2ci[3:, 3:, 2] = d2d_dpsii_dpsij - d2e_dpsii_dpsij 815 816 # Weight c1. 817 data.d2ci[3:, 3:, 3] = 6.0 * beta_gamma_comp 818 819 # Weight c2. 820 data.d2ci[3:, 3:, 4] = d2d_dpsii_dpsij + d2e_dpsii_dpsij 821 822 823 # psii-Da partial derivative. 824 ############################# 825 826 # Calculate d2e_dpsii_dDa. 827 d2e_dpsii_dDa = 1.0 / 3.0 * (data.de1_dDa * data.dc_dpsii_alpha - data.de2_dDa * data.dc_dpsii_beta - data.de2_dDa * data.dc_dpsii_gamma) 828 829 # Weight c0. 830 data.d2ci[3:, 1, 2] = data.d2ci[1, 3:, 2] = -d2e_dpsii_dDa 831 832 # Weight c2. 833 data.d2ci[3:, 1, 4] = data.d2ci[1, 3:, 4] = d2e_dpsii_dDa 834 835 836 # psii-Dr partial derivative. 837 ############################# 838 839 # Calculate d2e_dpsii_dDr. 840 d2e_dpsii_dDr = -1.0 / 3.0 * (data.de1_dDr * data.dc_dpsii_alpha - data.de2_dDr * data.dc_dpsii_beta - data.de2_dDr * data.dc_dpsii_gamma) 841 842 # Weight c0. 843 data.d2ci[3:, 2, 2] = data.d2ci[2, 3:, 2] = -d2e_dpsii_dDr 844 845 # Weight c2. 846 data.d2ci[3:, 2, 4] = data.d2ci[2, 3:, 4] = d2e_dpsii_dDr 847 848 849 # Da-Da partial derivative. 850 ########################### 851 852 # Components. 853 mu_five = data.mu**5 854 855 if data.mu == 0.0: 856 d2e1_dDa2 = 0.0 857 d2e2_dDa2 = 0.0 858 d2e3_dDa2 = 0.0 859 else: 860 d2e1_dDa2 = (6.0 * diff_data.params[1]**2 + 3.0 * diff_data.params[1] * diff_data.params[2] - diff_data.params[2]**2) * diff_data.params[2] / mu_five 861 d2e2_dDa2 = (6.0 * diff_data.params[1]**2 - 3.0 * diff_data.params[1] * diff_data.params[2] - diff_data.params[2]**2) * diff_data.params[2] / mu_five 862 d2e3_dDa2 = 6.0 * diff_data.params[1] * diff_data.params[2]**2 / mu_five 863 864 # Calculate d2e_dDa2. 865 d2e_dDa2 = -1.0 / 3.0 * (d2e1_dDa2 * data.c_alpha - d2e2_dDa2 * data.c_beta - d2e2_dDa2 * data.c_gamma) 866 867 # Weight c0. 868 data.d2ci[1, 1, 2] = -0.25 * d2e_dDa2 869 870 # Weight c2. 871 data.d2ci[1, 1, 4] = 0.25 * d2e_dDa2 872 873 874 # Da-Dr partial derivative. 875 ########################### 876 877 # Components. 878 if data.mu == 0.0: 879 d2e1_dDa_dDr = 0.0 880 d2e2_dDa_dDr = 0.0 881 d2e3_dDa_dDr = 0.0 882 else: 883 d2e1_dDa_dDr = (9.0 * diff_data.params[1]**3 + 6.0 * diff_data.params[1]**2 * diff_data.params[2] - 6.0 * diff_data.params[1] * diff_data.params[2]**2 - diff_data.params[2]**3) / mu_five 884 d2e2_dDa_dDr = (9.0 * diff_data.params[1]**3 - 6.0 * diff_data.params[1]**2 * diff_data.params[2] - 6.0 * diff_data.params[1] * diff_data.params[2]**2 + diff_data.params[2]**3) / mu_five 885 d2e3_dDa_dDr = 2.0 * (diff_data.params[1]**2 - diff_data.params[2]**2) * diff_data.params[2] / mu_five 886 887 # Calculate d2e_dDa2. 888 d2e_dDa_dDr = 1.0 / 9.0 * (d2e1_dDa_dDr * data.c_alpha - d2e2_dDa_dDr * data.c_beta - d2e2_dDa_dDr * data.c_gamma) 889 890 # Weight c0. 891 data.d2ci[1, 2, 2] = data.d2ci[2, 1, 2] = -0.25 * d2e_dDa_dDr 892 893 # Weight c2. 894 data.d2ci[1, 2, 4] = data.d2ci[2, 1, 4] = 0.25 * d2e_dDa_dDr 895 896 897 # Dr-Dr partial derivative. 898 ########################### 899 900 # Components. 901 if data.mu == 0.0: 902 d2e1_dDr2 = 0.0 903 d2e2_dDr2 = 0.0 904 d2e3_dDr2 = 0.0 905 else: 906 d2e1_dDr2 = (3.0 * diff_data.params[1]**2 - 9.0 * diff_data.params[1] * diff_data.params[2] - 2.0 * diff_data.params[2]**2) * diff_data.params[1] / mu_five 907 d2e2_dDr2 = (3.0 * diff_data.params[1]**2 + 9.0 * diff_data.params[1] * diff_data.params[2] - 2.0 * diff_data.params[2]**2) * diff_data.params[1] / mu_five 908 d2e3_dDr2 = 2.0 * (3.0 * diff_data.params[1]**2 - 2.0 * diff_data.params[2]**2) * diff_data.params[1] / mu_five 909 910 # Calculate d2e_dDr2. 911 d2e_dDr2 = -1.0 / 9.0 * (d2e1_dDr2 * data.c_alpha - d2e2_dDr2 * data.c_beta - d2e2_dDr2 * data.c_gamma) 912 913 # Weight c0. 914 data.d2ci[2, 2, 2] = -0.25 * d2e_dDr2 915 916 # Weight c2. 917 data.d2ci[2, 2, 4] = 0.25 * d2e_dDr2

918

Source Code for Module maths_fns.weights