1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2c Copyright (C) Bruno Pincon 3c 4c Copyright (C) 2012 - 2016 - Scilab Enterprises 5c 6c This file is hereby licensed under the terms of the GNU GPL v2.0, 7c pursuant to article 5.3.4 of the CeCILL v.2.1. 8c This file was originally licensed under the terms of the CeCILL v2.1, 9c and continues to be available under such terms. 10c For more information, see the COPYING file which you should have received 11c along with this program. 12 13 subroutine SplineCub(x, y, d, n, type, A_d, A_sd, qdy, lll) 14* 15* PURPOSE 16* computes a cubic spline interpolation function 17* in Hermite form (ie computes the derivatives d(i) of the 18* spline in each interpolation point (x(i), y(i))) 19* 20* ARGUMENTS 21* inputs : 22* n number of interpolation points (n >= 3) 23* x, y the n interpolation points, x must be in strict increasing order 24* type type of the spline : currently 25* type = 0 correspond to a NOT_A_KNOT spline where it is 26* imposed the conditions : 27* s'''(x(2)-) = s'''(x(2)+) 28* and s'''(x(n-1)-) = s'''(x(n-1)+) 29* type = 1 correspond to a NATURAL spline with the conditions : 30* s''(x1) = 0 31* and s''(xn) = 0 32* type = 2 correspond to a CLAMPED spline (d(1) and d(n) are given) 33* 34* type = 3 correspond to a PERIODIC spline 35* outputs : 36* d the derivatives in each x(i) i = 1..n 37* 38* work arrays : 39* A_d(1..n), A_sd(1..n-1), qdy(1..n-1) 40* lll(1..n-1) (used only in the periodic case) 41* 42* NOTES 43* this routine requires (i) n >= 3 (for natural) n >=4 (for not_a_knot) 44* (ii) strict increasing abscissae x(i) 45* (iii) y(1) = y(n) in the periodic case 46* THESE CONDITIONS MUST BE TESTED IN THE CALLING CODE 47* 48* AUTHOR 49* Bruno Pincon 50* 51* July 22 2004 : correction of the case not_a_knot which worked only 52* for equidistant abscissae 53* 54 implicit none 55 56 integer n, type 57 double precision x(n), y(n), d(n), A_d(n), A_sd(n-1), qdy(n-1), 58 $ lll(n-1) 59 60 include 'constinterp.h.f' 61 integer i 62 double precision r 63 64 if (n .eq. 2) then 65 if (type .ne. CLAMPED) then 66 d(1) = (y(2) - y(1)) / (x(2) - x(1)) 67 d(2) = d(1) 68 endif 69 return 70 endif 71 72 if (n .eq. 3 .and. type .eq. NOT_A_KNOT) then 73 call derivd(x, y, d, n, 1, FAST) 74 return 75 endif 76 77 do i = 1, n-1 78 A_sd(i) = 1.d0 / (x(i+1) - x(i)) 79 qdy(i) = (y(i+1) - y(i)) * A_sd(i)**2 80 enddo 81 82* compute the coef matrix and r.h.s. for rows 2..n-1 83* (which don't relies on the type) 84 do i = 2, n-1 85 A_d(i) = 2.d0*( A_sd(i-1) + A_sd(i) ) 86 d(i) = 3.d0 * ( qdy(i-1) + qdy(i) ) 87 enddo 88 89* compute equ 1 and n in function of the type 90 if (type .eq. NATURAL) then 91 A_d(1) = 2.d0*A_sd(1) 92 d(1) = 3.d0 * qdy(1) 93 A_d(n) = 2.d0*A_sd(n-1) 94 d(n) = 3.d0 * qdy(n-1) 95 call TriDiagLDLtSolve(A_d, A_sd, d, n) 96 97 else if (type .eq. NOT_A_KNOT) then 98 ! s'''(x(2)-) = s'''(x(2)+) 99 r = A_sd(2)/A_sd(1) 100 A_d(1) = A_sd(1)/(1.d0 + r) 101 d(1) = ((3.d0*r+2.d0)*qdy(1)+r*qdy(2))/(1.d0+r)**2 102 ! s'''(x(n-1)-) = s'''(x(n-1)+) 103 r = A_sd(n-2)/A_sd(n-1) 104 A_d(n) = A_sd(n-1)/(1.d0 + r) 105 d(n) = ((3.d0*r+2.d0)*qdy(n-1)+r*qdy(n-2))/(1.d0+r)**2 106 call TriDiagLDLtSolve(A_d, A_sd, d, n) 107 108 else if (type .eq. CLAMPED) then 109 ! d(1) and d(n) are already known 110 d(2) = d(2) - d(1)*A_sd(1) 111 d(n-1) = d(n-1) - d(n)*A_sd(n-1) 112 call TriDiagLDLtSolve(A_d(2), A_sd(2), d(2), n-2) 113 114 else if (type .eq. PERIODIC) then 115 A_d(1) = 2.d0*( A_sd(1) + A_sd(n-1) ) 116 d(1) = 3.d0 * ( qdy(1) + qdy(n-1) ) 117 lll(1) = A_sd(n-1) 118 call dset(n-2, 0.d0, lll(2),1) ! mise a zero 119 lll(n-2) = A_sd(n-2) 120 call CyclicTriDiagLDLtSolve(A_d, A_sd, lll, d, n-1) 121 d(n) = d(1) 122 123 endif 124 125 end ! subroutine SplineCub 126 127 128 subroutine TriDiagLDLtSolve(d, l, b, n) 129* 130* PURPOSE 131* solve a linear system A x = b with a symmetric tridiagonal positive definite 132* matrix A by using an LDL^t factorization 133* 134* PARAMETERS 135* d(1..n) : on input the diagonal of A 136* on output the diagonal of the (diagonal) matrix D 137* l(1..n-1) : on input the sub-diagonal of A 138* on output the sub-diagonal of L 139* b(1..n) : on input contains the r.h.s. b 140* on output the solution x 141* n : the dimension 142* 143* CAUTION 144* no zero pivot detection 145* 146 implicit none 147 integer n 148 integer flag 149 double precision d(n), l(n-1), b(n) 150 151 integer i 152 double precision temp 153 154 do i = 2, n 155 temp = l(i-1) 156 l(i-1) = l(i-1) / d(i-1) 157 d(i) = d(i) - temp * l(i-1) 158 b(i) = b(i) - l(i-1)*b(i-1) 159 enddo 160 161 b(n) = b(n) / d(n) 162 do i = n-1, 1, -1 163 b(i) = b(i)/d(i) - l(i)*b(i+1) 164 enddo 165 166 end ! subroutine TriDiagLDLtSolve 167 168 subroutine CyclicTriDiagLDLtSolve(d, lsd, lll, b, n) 169* 170* PURPOSE 171* solve a linear system A x = b with a symmetric "nearly" tridiagonal 172* positive definite matrix A by using an LDL^t factorization, 173* the matrix A has the form : 174* 175* |x x x| |1 | 176* |x x x | |x 1 | 177* | x x x | | x 1 | 178* | x x x | and so the L is like | x 1 | 179* | x x x | | x 1 | 180* | x x x| | x 1 | 181* |x x x| |x x x x x x 1| 182* 183* PARAMETERS 184* d(1..n) : on input the diagonal of A 185* on output the diagonal of the (diagonal) matrix D 186* lsd(1..n-2) : on input the sub-diagonal of A (without A(n,n-1)) 187* on output the sub-diagonal of L (without L(n,n-1)) 188* lll(1..n-1) : on input the last line of A (without A(n,n)) 189* on output the last line of L (without L(n,n)) 190* b(1..n) : on input contains the r.h.s. b 191* on output the solution x 192* n : the dimension 193* 194* CAUTION 195* no zero pivot detection 196* 197 implicit none 198 integer n 199 double precision d(n), lsd(n-2), lll(n-1), b(n) 200 201 integer i, j 202 double precision temp1, temp2 203 204* compute the LDL^t factorization 205 do i = 1, n-2 206 temp1 = lsd(i) 207 temp2 = lll(i) 208 lsd(i) = lsd(i) / d(i) ! elimination coef L(i,i-1) 209 lll(i) = lll(i) / d(i) ! elimination coef L(n,i-1) 210 211 d(i+1) = d(i+1) - lsd(i) * temp1 ! elimination on line i+1 212 lll(i+1) = lll(i+1) - lll(i) * temp1 ! elimination on line n 213 d(n) = d(n) - lll(i) * temp2 ! elimination on line n 214 enddo 215 temp2 = lll(n-1) 216 lll(n-1) = lll(n-1) / d(n-1) 217 d(n) = d(n) - lll(n-1) * temp2 218 219* solve LDL^t x = b (but use b for x and for the intermediary vectors...) 220 do i = 2, n-1 221 b(i) = b(i) - lsd(i-1)*b(i-1) 222 enddo 223 do j = 1, n-1 224 b(n) = b(n) - lll(j)*b(j) 225 enddo 226 227 do i = 1, n 228 b(i) = b(i) / d(i) 229 enddo 230 231 b(n-1) = b(n-1) - lll(n-1)*b(n) 232 do i = n-2, 1, -1 233 b(i) = b(i) - lsd(i)*b(i+1) - lll(i)*b(n) 234 enddo 235 236 end ! subroutine CyclicTriDiagLDLtSolve 237 238 239 integer function isearch(t, x, n) 240* 241* PURPOSE 242* x(1..n) being an array (with strict increasing order and n >=2) 243* representing intervals, this routine return i such that : 244* 245* x(i) <= t <= x(i+1) 246* 247* and 0 if t is not in [x(1), x(n)] 248* 249 implicit none 250 integer n 251 double precision t, x(n) 252 253 integer i1, i2, i 254 255 if ( x(1) .le. t .and. t .le. x(n)) then 256* dichotomic search 257 i1 = 1 258 i2 = n 259 do while ( i2 - i1 .gt. 1 ) 260 i = (i1 + i2)/2 261 if ( t .le. x(i) ) then 262 i2 = i 263 else 264 i1 = i 265 endif 266 enddo 267 isearch = i1 268 else 269 isearch = 0 270 endif 271 272 end 273 274 subroutine EvalPWHermite(t, st, dst, d2st, d3st, m, x, y, d, n, 275 $ outmode) 276* 277* PURPOSE 278* evaluation at the abscissae t(1..m) of the piecewise hermite function 279* define by x(1..n), y(1..n), d(1..n) (d being the derivatives at the 280* x(i)) together with its derivative, second derivative and third derivative 281* 282* PARAMETERS 283* 284* outmode : define what return in case t(j) not in [x(1), x(n)] 285 286 implicit none 287 integer m, n, outmode 288 double precision t(m), st(m), dst(m), d2st(m), d3st(m), 289 $ x(n), y(n), d(n) 290 291 include 'constinterp.h.f' 292 integer i, j 293 integer isearch, isanan 294 external isearch, isanan 295 double precision tt 296 external returnananfortran 297 logical new_call 298 common /INFO_HERMITE/new_call 299 300 new_call = .true. 301 i = 0 302 do j = 1, m 303 tt = t(j) 304 call fast_int_search(tt, x, n, i) ! recompute i only if necessary 305 306 if ( i .ne. 0 ) then 307 call EvalHermite(tt, x(i), x(i+1), y(i), y(i+1), d(i), 308 $ d(i+1), st(j), dst(j), d2st(j), d3st(j), i) 309 310 else ! t(j) is outside [x(1), x(n)] evaluation depend upon outmode 311 312 if (outmode .eq. BY_NAN .or. isanan(tt) .eq. 1) then 313 CALL returnananfortran(st(j)) 314 dst(j) = st(j) 315 d2st(j) = st(j) 316 d3st(j) = st(j) 317 318 elseif (outmode .eq. BY_ZERO) then 319 st(j) = 0.d0 320 dst(j) = 0.d0 321 d2st(j) = 0.d0 322 d3st(j) = 0.d0 323 324 elseif (outmode .eq. C0) then 325 dst(j) = 0.d0 326 d2st(j) = 0.d0 327 d3st(j) = 0.d0 328 if ( tt .lt. x(1) ) then 329 st(j) = y(1) 330 else 331 st(j) = y(n) 332 endif 333 334 elseif (outmode .eq. LINEAR) then 335 d2st(j) = 0.d0 336 d3st(j) = 0.d0 337 if ( tt .lt. x(1) ) then 338 dst(j) = d(1) 339 st(j) = y(1) + (tt - x(1))*d(1) 340 else 341 dst(j) = d(n) 342 st(j) = y(n) + (tt - x(n))*d(n) 343 endif 344 345 else 346 if (outmode .eq. NATURAL) then 347 call near_interval(tt, x, n, i) 348 elseif (outmode .eq. PERIODIC) then 349 call coord_by_periodicity(tt, x, n, i) 350 endif 351 call EvalHermite(tt, x(i), x(i+1), y(i), y(i+1), d(i), 352 $ d(i+1), st(j), dst(j), d2st(j), d3st(j), i) 353 endif 354 endif 355 enddo 356 357 end ! subroutine EvalPWHermite 358 359 subroutine EvalHermite(t, xa, xb, ya, yb, da, db, h, dh, ddh, 360 $ dddh, i) 361 362 implicit none 363 double precision t, xa, xb, ya, yb, da, db, h, dh, ddh, dddh 364 integer i 365 366 logical new_call 367 common /INFO_HERMITE/new_call 368 369 double precision tmxa, dx, p, c2, c3 370 integer old_i 371 save old_i, c2, c3 372 data old_i/0/ 373 374 if (old_i .ne. i .or. new_call) then 375* compute the following Newton form : 376* h(t) = ya + da*(t-xa) + c2*(t-xa)^2 + c3*(t-xa)^2*(t-xb) 377 dx = 1.d0/(xb - xa) 378 p = (yb - ya) * dx 379 c2 = (p - da) * dx 380 c3 = ( (db - p) + (da - p) ) * (dx*dx) 381 new_call = .false. 382 endif 383 old_i = i 384 385* eval h(t), h'(t), h"(t) and h"'(t), by a generalized Horner 's scheme 386 tmxa = t - xa 387 h = c2 + c3*(t - xb) 388 dh = h + c3*tmxa 389 ddh = 2.d0*( dh + c3*tmxa ) 390 dddh = 6.d0*c3 391 h = da + h*tmxa 392 dh = h + dh*tmxa 393 h = ya + h*tmxa 394 395 end ! subroutine EvalHermite 396 397 subroutine fast_int_search(xx, x, nx, i) 398 399 implicit none 400 integer nx, i 401 double precision xx, x(nx) 402 403 integer isearch 404 external isearch 405 406 if ( i .eq. 0 ) then 407 i = isearch(xx, x, nx) 408 elseif ( .not. (x(i) .le. xx .and. xx .le. x(i+1))) then 409 i = isearch(xx, x, nx) 410 endif 411 412 end 413 414 subroutine coord_by_periodicity(t, x, n, i) 415* 416* PURPOSE 417* recompute t such that t in [x(1), x(n)] by periodicity : 418* and then the interval i of this new t 419* 420 implicit none 421 integer n, i 422 double precision t, x(n) 423 integer isearch 424 external isearch 425 426 double precision r, dx 427 428 dx = x(n) - x(1) 429 r = (t - x(1)) / dx 430 431 if (r .ge. 0.d0) then 432 t = x(1) + (r - aint(r))*dx 433 else 434 r = abs(r) 435 t = x(n) - (r - aint(r))*dx 436 endif 437 438 ! some cautions in case of roundoff errors (is necessary ?) 439 if (t .lt. x(1)) then 440 t = x(1) 441 i = 1 442 else if (t .gt. x(n)) then 443 t = x(n) 444 i = n-1 445 else 446 i = isearch(t, x, n) 447 endif 448 449 end ! subroutine coord_by_periodicity 450 451 452 subroutine near_grid_point(xx, x, nx, i) 453* 454* calcule le point de la grille le plus proche ... a detailler 455* 456 implicit none 457 integer i, nx 458 double precision xx, x(nx) 459 460 if (xx .lt. x(1)) then 461 i = 1 462 xx = x(1) 463 else ! xx > x(nx) 464 i = nx-1 465 xx = x(nx) 466 endif 467 468 end 469 470 subroutine near_interval(xx, x, nx, i) 471* 472* idem sans modifier xx 473* 474 implicit none 475 integer i, nx 476 double precision xx, x(nx) 477 478 if (xx .lt. x(1)) then 479 i = 1 480 else ! xx > x(nx) 481 i = nx-1 482 endif 483 484 end 485 486 subroutine proj_by_per(t, xmin, xmax) 487* 488* PURPOSE 489* recompute t such that t in [xmin, xmax] by periodicity. 490* 491 implicit none 492 double precision t, xmin, xmax 493 double precision r, dx 494 495 dx = xmax - xmin 496 r = (t - xmin) / dx 497 498 if (r .ge. 0.d0) then 499 t = xmin + (r - aint(r))*dx 500 else 501 r = abs(r) 502 t = xmax - (r - aint(r))*dx 503 endif 504 505 ! some cautions in case of roundoff errors (is necessary ?) 506 if (t .lt. xmin) then 507 t = xmin 508 else if (t .gt. xmax) then 509 t = xmax 510 endif 511 512 end ! subroutine proj_by_per 513 514 515 subroutine proj_on_grid(xx, xmin, xmax) 516* 517 implicit none 518 double precision xx, xmin, xmax 519 520 if (xx .lt. xmin) then 521 xx = xmin 522 else 523 xx = xmax 524 endif 525 526 end 527 528 subroutine BiCubicSubSpline(x, y, u, nx, ny, C, p, q, r, type) 529* 530* PURPOSE 531* compute bicubic subsplines 532* 533 implicit none 534 integer nx, ny, type 535 double precision x(nx), y(ny), u(nx, ny), C(4,4,nx-1,ny-1), 536 $ p(nx, ny), q(nx, ny), r(nx, ny) 537 integer i, j 538 include 'constinterp.h.f' 539 540 if (type .eq. MONOTONE) then 541* approximation des derivees par SUBROUTINE DPCHIM(N,X,F,D,INCFD) 542 ! p = du/dx 543 do j = 1, ny 544 call dpchim(nx, x, u(1,j), p(1,j), 1) 545 enddo 546 ! q = du/dy 547 do i = 1, nx 548 call dpchim(ny, y, u(i,1), q(i,1), nx) 549 enddo 550 ! r = d2 u/ dx dy approchee via dq / dx 551 do j = 1, ny 552 call dpchim(nx, x, q(1,j), r(1,j), 1) 553 enddo 554 555 elseif (type .eq. FAST .or. type .eq. FAST_PERIODIC) then 556* approximation des derivees partielles par methode simple 557 558 ! p = du/dx 559 do j = 1, ny 560 call derivd(x, u(1,j), p(1,j), nx, 1, type) 561 enddo 562 ! q = du/dy 563 do i = 1, nx 564 call derivd(y, u(i,1), q(i,1), ny, nx, type) 565 enddo 566 ! r = d2 u/ dx dy approchee via dq / dx 567 do j = 1, ny 568 call derivd(x, q(1,j), r(1,j), nx, 1, type) 569 enddo 570 571 endif 572 573* calculs des coefficients dans les bases (x-x(i))^k (y-y(j))^l 0<= k,l <= 3 574* pour evaluation rapide via Horner par la suite 575 call coef_bicubic(u, p, q, r, x, y, nx, ny, C) 576 577 end ! subroutine BiCubicSubSpline 578 579 subroutine BiCubicSpline(x, y, u, nx, ny, C, p, q, r, 580 $ A_d, A_sd, d, ll, qdu, u_temp, type) 581* 582* PURPOSE 583* compute bicubic splines 584* 585 implicit none 586 integer nx, ny, type 587 double precision x(nx), y(ny), u(nx, ny), C(4,4,nx-1,ny-1), 588 $ p(nx, ny), q(nx, ny), r(nx, ny), A_d(*), 589 $ A_sd(*), d(ny), ll(*), qdu(*), u_temp(ny) 590 include 'constinterp.h.f' 591 integer i, j 592 593 ! compute du/dx 594 do j = 1, ny 595 call SplineCub(x, u(1,j), p(1,j), nx, type, A_d, A_sd, qdu, ll) 596 enddo 597 598 ! compute du/dy 599 do i = 1, nx 600 call dcopy(ny, u(i,1), nx, u_temp, 1) 601 call SplineCub(y, u_temp, d, ny, type, A_d, A_sd, qdu, ll) 602 call dcopy(ny, d, 1, q(i,1), nx) 603 enddo 604 605 ! compute ddu/dxdy 606 call SplineCub(x, q(1,1), r(1,1), nx, type, A_d, A_sd, qdu, ll) 607 call SplineCub(x, q(1,ny), r(1,ny), nx, type, A_d, A_sd, qdu, ll) 608 609 do i = 1, nx 610 call dcopy(ny, p(i,1), nx, u_temp, 1) 611 d(1) = r(i,1) 612 d(ny) = r(i,ny) 613 call SplineCub(y, u_temp, d, ny, CLAMPED, A_d, A_sd, qdu, ll) 614 call dcopy(ny-2, d(2), 1, r(i,2), nx) 615 enddo 616 617* calculs des coefficients dans les bases (x-x(i))^k (y-y(j))^l 0<= k,l <= 3 618* pour evaluation rapide via Horner par la suite 619 call coef_bicubic(u, p, q, r, x, y, nx, ny, C) 620 621 end ! subroutine BiCubicSpline 622 623 624 subroutine derivd(x, u, du, n, inc, type) 625* 626* PURPOSE 627* given functions values u(i) at points x(i), i = 1, ..., n 628* this subroutine computes approximations du(i) of the derivative 629* at the points x(i). 630* 631* METHOD 632* For i in [2,n-1], the "centered" formula of order 2 is used : 633* d(i) = derivative at x(i) of the interpolation polynomial 634* of the points {(x(j),u(j)), j in [i-1,i+1]} 635* 636* For i=1 and n, if type = FAST_PERIODIC (in which case u(n)=u(1)) then 637* the previus "centered" formula is also used else (type = FAST), d(1) 638* is the derivative at x(1) of the interpolation polynomial of 639* {(x(j),u(j)), j in [1,3]} and the same method is used for d(n) 640* 641* ARGUMENTS 642* inputs : 643* n integer : number of point (n >= 2) 644* x, u double precision : the n points, x must be in strict increasing order 645* type integer : FAST (the function is non periodic) or FAST_PERIODIC 646* (the function is periodic), in this last case u(n) must be equal to u(1)) 647* inc integer : to deal easily with 2d applications, u(i) is in fact 648* u(1,i) with u declared as u(inc,*) to avoid the direct management of 649* the increment inc (the i th value given with u(1 + inc*(i-1) ...) 650* outputs : 651* d the derivatives in each x(i) i = 1..n 652* 653* NOTES 654* this routine requires (i) n >= 2 655* (ii) strict increasing abscissae x(i) 656* (iii) u(1)=u(n) if type = FAST_PERIODIC 657* ALL THESE CONDITIONS MUST BE TESTED IN THE CALLING CODE 658* 659* AUTHOR 660* Bruno Pincon 661* 662 663 implicit none 664 integer n, inc, type 665 double precision x(n), u(inc,*), du(inc,*) 666 667 include 'constinterp.h.f' 668 double precision dx_l, du_l, dx_r, du_r, w_l, w_r 669 integer i, k 670 671 672 if (n .eq. 2) then ! special case used linear interp 673 du(1,1) = (u(1,2) - u(1,1))/(x(2)-x(1)) 674 du(1,2) = du(1,1) 675 return 676 endif 677 678 if (type .eq. FAST_PERIODIC) then 679 680 dx_r = x(n)-x(n-1) 681 du_r = (u(1,1) - u(1,n-1)) / dx_r 682 do i = 1, n-1 683 dx_l = dx_r 684 du_l = du_r 685 dx_r = x(i+1) - x(i) 686 du_r = (u(1,i+1) - u(1,i))/dx_r 687 w_l = dx_r/(dx_l + dx_r) 688 w_r = 1.d0 - w_l 689 du(1,i) = w_l*du_l + w_r*du_r 690 enddo 691 du(1,n) = du(1,1) 692 693 else if (type .eq. FAST) then 694 695 dx_l = x(2) - x(1) 696 du_l = (u(1,2) - u(1,1))/dx_l 697 dx_r = x(3) - x(2) 698 du_r = (u(1,3) - u(1,2))/dx_r 699 w_l = dx_r/(dx_l + dx_r) 700 w_r = 1.d0 - w_l 701 du(1,1) = (1.d0 + w_r)*du_l - w_r*du_r 702 du(1,2) = w_l*du_l + w_r*du_r 703 do i = 3, n-1 704 dx_l = dx_r 705 du_l = du_r 706 dx_r = x(i+1) - x(i) 707 du_r = (u(1,i+1) - u(1,i))/dx_r 708 w_l = dx_r/(dx_l + dx_r) 709 w_r = 1.d0 - w_l 710 du(1,i) = w_l*du_l + w_r*du_r 711 enddo 712 du(1,n) = (1.d0 + w_l)*du_r - w_l*du_l 713 714 endif 715 716 end 717 718 719 subroutine coef_bicubic(u, p, q, r, x, y, nx, ny, C) 720* 721* PURPOSE 722* compute for each polynomial (i,j)-patch (defined on 723* [x(i),x(i+1)]x[y(i),y(i+1)]) the following base 724* representation : 725* i,j _4_ _4_ i,j k-1 l-1 726* u (x,y) = >__ >__ C (x-x(i)) (y-y(j)) 727* k=1 l=1 k,l 728* 729* from the "Hermite" representation (values of u, p = du/dx, 730* q = du/dy, r = ddu/dxdy at the 4 vertices (x(i),y(j)), 731* (x(i+1),y(j)), (x(i+1),y(j+1)), (x(i),y(j+1)). 732* 733 implicit none 734 735 integer nx, ny 736 double precision u(nx, ny), p(nx, ny), q(nx, ny), r(nx, ny), 737 $ x(nx), y(ny), C(4,4,nx-1,ny-1) 738 739 integer i, j 740 double precision a, b, cc, d, dx, dy 741 742 do j = 1, ny-1 743 dy = 1.d0/(y(j+1) - y(j)) 744 745 do i = 1, nx-1 746 dx = 1.d0/(x(i+1) - x(i)) 747 748 C(1,1,i,j) = u(i,j) 749 C(2,1,i,j) = p(i,j) 750 C(1,2,i,j) = q(i,j) 751 C(2,2,i,j) = r(i,j) 752 753 a = (u(i+1,j) - u(i,j))*dx 754 C(3,1,i,j) = (3.d0*a -2.d0*p(i,j) - p(i+1,j))*dx 755 C(4,1,i,j) = (p(i+1,j) + p(i,j) - 2.d0*a)*(dx*dx) 756 757 a = (u(i,j+1) - u(i,j))*dy 758 C(1,3,i,j) = (3.d0*a -2.d0*q(i,j) - q(i,j+1))*dy 759 C(1,4,i,j) = (q(i,j+1) + q(i,j) - 2.d0*a)*(dy*dy) 760 761 a = (q(i+1,j) - q(i,j))*dx 762 C(3,2,i,j) = (3.d0*a - r(i+1,j) - 2.d0*r(i,j))*dx 763 C(4,2,i,j) = (r(i+1,j) + r(i,j) - 2.d0*a)*(dx*dx) 764 765 a = (p(i,j+1) - p(i,j))*dy 766 C(2,3,i,j) = (3.d0*a - r(i,j+1) - 2.d0*r(i,j))*dy 767 C(2,4,i,j) = (r(i,j+1) + r(i,j) - 2.d0*a)*(dy*dy) 768 769 a = (u(i+1,j+1)+u(i,j)-u(i+1,j)-u(i,j+1))*(dx*dx*dy*dy) 770 $ - (p(i,j+1)-p(i,j))*(dx*dy*dy) 771 $ - (q(i+1,j)-q(i,j))*(dx*dx*dy) 772 $ + r(i,j)*(dx*dy) 773 b = (p(i+1,j+1)+p(i,j)-p(i+1,j)-p(i,j+1))*(dx*dy*dy) 774 $ - (r(i+1,j)-r(i,j))*(dx*dy) 775 cc = (q(i+1,j+1)+q(i,j)-q(i+1,j)-q(i,j+1))*(dx*dx*dy) 776 $ - (r(i,j+1)-r(i,j))*(dx*dy) 777 d = (r(i+1,j+1)+r(i,j)-r(i+1,j)-r(i,j+1))*(dx*dy) 778 779 780 C(3,3,i,j) = 9.d0*a - 3.d0*b - 3.d0*cc + d 781 C(3,4,i,j) =(-6.d0*a + 2.d0*b + 3.d0*cc - d)*dy 782 C(4,3,i,j) =(-6.d0*a + 3.d0*b + 2.d0*cc - d)*dx 783 C(4,4,i,j) =( 4.d0*a - 2.d0*b - 2.d0*cc + d)*dx*dy 784 785 enddo 786 enddo 787 end 788 789 double precision function EvalBicubic(xx, yy, xk, yk, Ck) 790 791 implicit none 792 double precision xx, yy, xk, yk, Ck(4,4) 793 794 double precision dx, dy, u 795 integer i 796 797 dx = xx - xk 798 dy = yy - yk 799 800 u = 0.d0 801 do i = 4, 1, -1 802 u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx 803 enddo 804 EvalBicubic = u 805 806 end ! function EvalBicubic 807 808 subroutine EvalBicubic_with_grad(xx, yy, xk, yk, Ck, 809 $ u, dudx, dudy) 810 811 implicit none 812 double precision xx, yy, xk, yk, Ck(4,4), u, dudx, dudy 813 814 double precision dx, dy 815 integer i 816 817 dx = xx - xk 818 dy = yy - yk 819 u = 0.d0 820 dudx = 0.d0 821 dudy = 0.d0 822 do i = 4, 1, -1 823 u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx 824 dudx = Ck(2,i)+dx*(2.d0*Ck(3,i) + dx*3.d0*Ck(4,i)) + dudx*dy 825 dudy = Ck(i,2)+dy*(2.d0*Ck(i,3) + dy*3.d0*Ck(i,4)) + dudy*dx 826 enddo 827 828 end ! subroutine EvalBicubic_with_grad 829 830 831 subroutine EvalBicubic_with_grad_and_hes(xx, yy, xk, yk, Ck, 832 $ u, dudx, dudy, 833 $ d2udx2, d2udxy, d2udy2) 834 835 implicit none 836 double precision xx, yy, xk, yk, Ck(4,4), u, dudx, dudy, d2udx2, 837 $ d2udxy, d2udy2 838 839 double precision dx, dy 840 integer i 841 842 dx = xx - xk 843 dy = yy - yk 844 u = 0.d0 845 dudx = 0.d0 846 dudy = 0.d0 847 d2udx2 = 0.d0 848 d2udy2 = 0.d0 849 d2udxy = 0.d0 850 do i = 4, 1, -1 851 u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx 852 dudx = Ck(2,i)+dx*(2.d0*Ck(3,i) + dx*3.d0*Ck(4,i)) + dudx*dy 853 dudy = Ck(i,2)+dy*(2.d0*Ck(i,3) + dy*3.d0*Ck(i,4)) + dudy*dx 854 d2udx2 = 2.d0*Ck(3,i) + dx*6.d0*Ck(4,i) + d2udx2*dy 855 d2udy2 = 2.d0*Ck(i,3) + dy*6.d0*Ck(i,4) + d2udy2*dx 856 enddo 857 d2udxy = Ck(2,2)+dy*(2.d0*Ck(2,3)+dy*3.d0*Ck(2,4)) 858 . + dx*(2.d0*(Ck(3,2)+dy*(2.d0*Ck(3,3)+dy*3.d0*Ck(3,4))) 859 . + dx*(3.d0*(Ck(4,2)+dy*(2.d0*Ck(4,3)+dy*3.d0*Ck(4,4))))) 860 end ! subroutine EvalBicubic_with_grad_and_hes 861 862 863 subroutine BiCubicInterp(x, y, C, nx, ny, x_eval, y_eval, z_eval, 864 $ m, outmode) 865* 866* PURPOSE 867* bicubic interpolation : 868* the grid is defined by x(1..nx), y(1..ny) 869* the known values are z(1..nx,1..ny), (z(i,j) being the value 870* at point (x(i),y(j))) 871* the interpolation is done on the points x_eval,y_eval(1..m) 872* z_eval(k) is the result of the bicubic interpolation of 873* (x_eval(k), y_eval(k)) 874* 875 implicit none 876 integer nx, ny, m, outmode 877 double precision x(nx), y(ny), C(4,4,nx-1,ny-1), 878 $ x_eval(m), y_eval(m), z_eval(m) 879 880 double precision xx, yy 881 integer i, j, k 882 include 'constinterp.h.f' 883 integer isanan 884 double precision EvalBicubic 885 external isanan, returnananfortran, EvalBicubic 886 887 i = 0 888 j = 0 889 do k = 1, m 890 xx = x_eval(k) 891 call fast_int_search(xx, x, nx, i) 892 yy = y_eval(k) 893 call fast_int_search(yy, y, ny, j) 894 895 if (i .ne. 0 .and. j .ne. 0) then 896 z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j)) 897 898 elseif (outmode .eq. BY_NAN .or. isanan(xx) .eq. 1 899 $ .or. isanan(yy) .eq. 1) then 900 CALL returnananfortran(z_eval(k)) 901 902 elseif (outmode .eq. BY_ZERO) then 903 z_eval(k) = 0.d0 904 905 elseif (outmode .eq. PERIODIC) then 906 if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i) 907 if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j) 908 z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j)) 909 910 elseif (outmode .eq. C0) then 911 if (i .eq. 0) call near_grid_point(xx, x, nx, i) 912 if (j .eq. 0) call near_grid_point(yy, y, ny, j) 913 z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j)) 914 915 elseif (outmode .eq. NATURAL) then 916 if (i .eq. 0) call near_interval(xx, x, nx, i) 917 if (j .eq. 0) call near_interval(yy, y, ny, j) 918 z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j)) 919 endif 920 921 enddo 922 923 end 924 925 subroutine BiCubicInterpWithGrad(x, y, C, nx, ny, x_eval, y_eval, 926 $ z_eval, dzdx_eval, dzdy_eval,m, 927 $ outmode) 928* 929* PURPOSE 930* bicubic interpolation : 931* the grid is defined by x(1..nx), y(1..ny) 932* the known values are z(1..nx,1..ny), (z(i,j) being the value 933* at point (x(i),y(j))) 934* the interpolation is done on the points x_eval,y_eval(1..m) 935* z_eval(k) is the result of the bicubic interpolation of 936* (x_eval(k), y_eval(k)) and dzdx_eval(k), dzdy_eval(k) is the gradient 937* 938 implicit none 939 integer nx, ny, m, outmode 940 double precision x(nx), y(ny), C(4,4,nx-1,ny-1), 941 $ x_eval(m), y_eval(m), z_eval(m), 942 $ dzdx_eval(m), dzdy_eval(m) 943 944 double precision xx, yy 945 integer k, i, j 946 include 'constinterp.h.f' 947 integer isanan 948 double precision EvalBicubic 949 external isanan, returnananfortran, EvalBicubic 950 logical change_dzdx, change_dzdy 951 952 i = 0 953 j = 0 954 do k = 1, m 955 xx = x_eval(k) 956 call fast_int_search(xx, x, nx, i) 957 yy = y_eval(k) 958 call fast_int_search(yy, y, ny, j) 959 960 if (i .ne. 0 .and. j .ne. 0) then 961 call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j), 962 $ z_eval(k), 963 $ dzdx_eval(k), dzdy_eval(k)) 964 965 elseif ( outmode .eq. BY_NAN .or. isanan(xx) .eq. 1 966 $ .or. isanan(yy) .eq. 1) then 967 CALL returnananfortran(z_eval(k)) 968 dzdx_eval(k) = z_eval(k) 969 dzdy_eval(k) = z_eval(k) 970 971 elseif (outmode .eq. BY_ZERO ) then 972 z_eval(k) = 0.d0 973 dzdx_eval(k) = 0.d0 974 dzdy_eval(k) = 0.d0 975 976 elseif (outmode .eq. PERIODIC) then 977 if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i) 978 if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j) 979 call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j), 980 $ z_eval(k), 981 $ dzdx_eval(k), dzdy_eval(k)) 982 983 elseif (outmode .eq. C0) then 984 if (i .eq. 0) then 985 call near_grid_point(xx, x, nx, i) 986 change_dzdx = .true. 987 else 988 change_dzdx = .false. 989 endif 990 if (j .eq. 0) then 991 call near_grid_point(yy, y, ny, j) 992 change_dzdy = .true. 993 else 994 change_dzdy = .false. 995 endif 996 call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j), 997 $ z_eval(k), 998 $ dzdx_eval(k), dzdy_eval(k)) 999 if (change_dzdx) dzdx_eval(k) = 0.d0 1000 if (change_dzdy) dzdy_eval(k) = 0.d0 1001 1002 elseif (outmode .eq. NATURAL) then 1003 if (i .eq. 0) call near_interval(xx, x, nx, i) 1004 if (j .eq. 0) call near_interval(yy, y, ny, j) 1005 call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j), 1006 $ z_eval(k), 1007 $ dzdx_eval(k), dzdy_eval(k)) 1008 1009 endif 1010 1011 enddo 1012 1013 end 1014 1015 subroutine BiCubicInterpWithGradAndHes(x, y, C, nx, ny, 1016 $ x_eval, y_eval, z_eval, 1017 $ dzdx_eval, dzdy_eval, 1018 $ d2zdx2_eval, d2zdxy_eval, d2zdy2_eval, 1019 $ m, outmode) 1020* 1021* PURPOSE 1022* bicubic interpolation : 1023* the grid is defined by x(1..nx), y(1..ny) 1024* the known values are z(1..nx,1..ny), (z(i,j) being the value 1025* at point (x(i),y(j))) 1026* the interpolation is done on the points x_eval,y_eval(1..m) 1027* z_eval(k) is the result of the bicubic interpolation of 1028* (x_eval(k), y_eval(k)), [dzdx_eval(k), dzdy_eval(k)] is the gradient 1029* and [d2zdx2(k), d2zdxy(k), d2zdy2(k)] the Hessian 1030* 1031 implicit none 1032 integer nx, ny, m, outmode 1033 double precision x(nx), y(ny), C(4,4,nx-1,ny-1), 1034 $ x_eval(m), y_eval(m), z_eval(m), dzdx_eval(m), 1035 $ dzdy_eval(m), d2zdx2_eval(m), d2zdxy_eval(m), d2zdy2_eval(m) 1036 1037 double precision xx, yy 1038 integer k, i, j 1039 include 'constinterp.h.f' 1040 integer isanan 1041 double precision EvalBicubic 1042 external isanan, returnananfortran, EvalBicubic 1043 logical change_dzdx, change_dzdy 1044 1045 i = 0 1046 j = 0 1047 do k = 1, m 1048 xx = x_eval(k) 1049 call fast_int_search(xx, x, nx, i) 1050 yy = y_eval(k) 1051 call fast_int_search(yy, y, ny, j) 1052 1053 if (i .ne. 0 .and. j .ne. 0) then 1054 call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j), 1055 $ C(1,1,i,j), z_eval(k), 1056 $ dzdx_eval(k), dzdy_eval(k), 1057 $ d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k)) 1058 1059 elseif ( outmode .eq. BY_NAN .or. isanan(xx) .eq. 1 1060 $ .or. isanan(yy) .eq. 1) then 1061 CALL returnananfortran(z_eval(k)) 1062 dzdx_eval(k) = z_eval(k) 1063 dzdy_eval(k) = z_eval(k) 1064 d2zdx2_eval(k) = z_eval(k) 1065 d2zdxy_eval(k) = z_eval(k) 1066 d2zdy2_eval(k) = z_eval(k) 1067 1068 elseif (outmode .eq. BY_ZERO ) then 1069 z_eval(k) = 0.d0 1070 dzdx_eval(k) = 0.d0 1071 dzdy_eval(k) = 0.d0 1072 d2zdx2_eval(k) = 0.d0 1073 d2zdxy_eval(k) = 0.d0 1074 d2zdy2_eval(k) = 0.d0 1075 1076 elseif (outmode .eq. PERIODIC) then 1077 if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i) 1078 if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j) 1079 call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j), 1080 $ C(1,1,i,j), z_eval(k), 1081 $ dzdx_eval(k), dzdy_eval(k), 1082 $ d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k)) 1083 1084 elseif (outmode .eq. C0) then 1085 if (i .eq. 0) then 1086 call near_grid_point(xx, x, nx, i) 1087 change_dzdx = .true. 1088 else 1089 change_dzdx = .false. 1090 endif 1091 if (j .eq. 0) then 1092 call near_grid_point(yy, y, ny, j) 1093 change_dzdy = .true. 1094 else 1095 change_dzdy = .false. 1096 endif 1097 call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j), 1098 $ C(1,1,i,j), z_eval(k), 1099 $ dzdx_eval(k), dzdy_eval(k), 1100 $ d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k)) 1101 if (change_dzdx) then 1102 dzdx_eval(k) = 0.d0 1103 d2zdx2_eval(k) = 0.d0 1104 d2zdxy_eval(k) = 0.d0 1105 endif 1106 if (change_dzdy) then 1107 dzdy_eval(k) = 0.d0 1108 d2zdxy_eval(k) = 0.d0 1109 d2zdy2_eval(k) = 0.d0 1110 endif 1111 1112 elseif (outmode .eq. NATURAL) then 1113 if (i .eq. 0) call near_interval(xx, x, nx, i) 1114 if (j .eq. 0) call near_interval(yy, y, ny, j) 1115 call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j), 1116 $ C(1,1,i,j), z_eval(k), 1117 $ dzdx_eval(k), dzdy_eval(k), 1118 $ d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k)) 1119 endif 1120 1121 enddo 1122 1123 end 1124 1125 subroutine driverdb3val(xp, yp, zp, fp, np, tx, ty, tz, 1126 $ nx, ny, nz, kx, ky, kz, bcoef, work, 1127 $ xmin, xmax, ymin, ymax, zmin, zmax, 1128 $ outmode) 1129* 1130* PURPOSE 1131* driver on to db3val 1132* 1133 implicit none 1134 integer np, nx, ny, nz, kx, ky, kz, outmode 1135 double precision xp(np), yp(np), zp(np), fp(np), 1136 $ tx(*), ty(*), tz(*), bcoef(*), work(*), 1137 $ xmin, xmax, ymin, ymax, zmin, zmax 1138 1139 integer k 1140 logical flag_x, flag_y, flag_z 1141 double precision x, y, z 1142 include 'constinterp.h.f' 1143 integer isanan 1144 double precision db3val 1145 external isanan, returnananfortran, db3val 1146 1147 do k = 1, np 1148 x = xp(k) 1149 if ( xmin .le. x .and. x .le. xmax ) then 1150 flag_x = .true. 1151 else 1152 flag_x = .false. 1153 endif 1154 y = yp(k) 1155 if ( ymin .le. y .and. y .le. ymax ) then 1156 flag_y = .true. 1157 else 1158 flag_y = .false. 1159 endif 1160 z = zp(k) 1161 if ( zmin .le. z .and. z .le. zmax ) then 1162 flag_z = .true. 1163 else 1164 flag_z = .false. 1165 endif 1166 1167 if ( flag_x .and. flag_y .and. flag_z ) then 1168 fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, 1169 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1170 1171 elseif (outmode .eq. BY_NAN .or. isanan(x) .eq. 1 1172 $ .or. isanan(y) .eq. 1 1173 $ .or. isanan(z) .eq. 1) then 1174 CALL returnananfortran(fp(k)) 1175 1176 elseif (outmode .eq. BY_ZERO) then 1177 fp(k) = 0.d0 1178 1179 else 1180 if (outmode .eq. PERIODIC) then 1181 if (.not. flag_x) call proj_by_per(x, xmin, xmax) 1182 if (.not. flag_y) call proj_by_per(y, ymin, ymax) 1183 if (.not. flag_z) call proj_by_per(z, zmin, zmax) 1184 elseif (outmode .eq. C0) then 1185 if (.not. flag_x) call proj_on_grid(x, xmin, xmax) 1186 if (.not. flag_y) call proj_on_grid(y, ymin, ymax) 1187 if (.not. flag_z) call proj_on_grid(z, zmin, zmax) 1188 endif 1189 fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz, 1190 $ kx, ky, kz, bcoef, work) 1191 endif 1192 enddo 1193 1194 end 1195 1196 subroutine driverdb3valwithgrad(xp, yp, zp, fp, 1197 $ dfpdx, dfpdy, dfpdz, np, tx, ty, tz, 1198 $ nx, ny, nz, kx, ky, kz, bcoef, work, 1199 $ xmin, xmax, ymin, ymax, zmin, zmax, 1200 $ outmode) 1201* 1202* PURPOSE 1203* driver on to db3val with gradient computing 1204* 1205 implicit none 1206 integer np, nx, ny, nz, kx, ky, kz, outmode 1207 double precision xp(np), yp(np), zp(np), fp(np), 1208 $ dfpdx(np), dfpdy(np), dfpdz(np), 1209 $ tx(*), ty(*), tz(*), bcoef(*), work(*), 1210 $ xmin, xmax, ymin, ymax, zmin, zmax 1211 1212 integer k 1213 logical flag_x, flag_y, flag_z 1214 double precision x, y, z 1215 include 'constinterp.h.f' 1216 integer isanan 1217 double precision db3val 1218 external isanan, returnananfortran, db3val 1219 1220 do k = 1, np 1221 x = xp(k) 1222 if ( xmin .le. x .and. x .le. xmax ) then 1223 flag_x = .true. 1224 else 1225 flag_x = .false. 1226 endif 1227 y = yp(k) 1228 if ( ymin .le. y .and. y .le. ymax ) then 1229 flag_y = .true. 1230 else 1231 flag_y = .false. 1232 endif 1233 z = zp(k) 1234 if ( zmin .le. z .and. z .le. zmax ) then 1235 flag_z = .true. 1236 else 1237 flag_z = .false. 1238 endif 1239 1240 if ( flag_x .and. flag_y .and. flag_z ) then 1241 fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, 1242 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1243 dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz, 1244 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1245 dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz, 1246 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1247 dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz, 1248 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1249 1250 elseif (outmode .eq. BY_NAN .or. isanan(x) .eq. 1 1251 $ .or. isanan(y) .eq. 1 1252 $ .or. isanan(z) .eq. 1) then 1253 CALL returnananfortran(fp(k)) 1254 dfpdx(k) = fp(k) 1255 dfpdy(k) = fp(k) 1256 dfpdz(k) = fp(k) 1257 1258 elseif (outmode .eq. BY_ZERO) then 1259 fp(k) = 0.d0 1260 dfpdx(k) = 0.d0 1261 dfpdy(k) = 0.d0 1262 dfpdz(k) = 0.d0 1263 1264 elseif (outmode .eq. PERIODIC) then 1265 if (.not. flag_x) call proj_by_per(x, xmin, xmax) 1266 if (.not. flag_y) call proj_by_per(y, ymin, ymax) 1267 if (.not. flag_z) call proj_by_per(z, zmin, zmax) 1268 fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz, 1269 $ kx, ky, kz, bcoef, work) 1270 dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz, 1271 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1272 dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz, 1273 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1274 dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz, 1275 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1276 1277 elseif (outmode .eq. C0) then 1278 if (.not. flag_x) call proj_on_grid(x, xmin, xmax) 1279 if (.not. flag_y) call proj_on_grid(y, ymin, ymax) 1280 if (.not. flag_z) call proj_on_grid(z, zmin, zmax) 1281 fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz, 1282 $ kx, ky, kz, bcoef, work) 1283 if ( flag_x ) then 1284 dfpdx(k) = 0.d0 1285 else 1286 dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz, 1287 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1288 endif 1289 if ( flag_y ) then 1290 dfpdy(k) = 0.d0 1291 else 1292 dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz, 1293 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1294 endif 1295 if ( flag_z ) then 1296 dfpdz(k) = 0.d0 1297 else 1298 dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz, 1299 $ nx, ny, nz, kx, ky, kz, bcoef, work) 1300 endif 1301 1302 endif 1303 enddo 1304 1305 end 1306