1! 2! Copyright (C) 2002-2005 FPMD-CPV groups 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!-----------------------------------------------------------------------------! 9! This module is basad on a similar module from CP2K 10!-----------------------------------------------------------------------------! 11 12 MODULE splines 13 14! routines for handling splines 15! allocate_spline: allocates x and y vectors for splines 16! init_spline: generate table for spline (allocate spl%y2) 17! spline: return value of spline for given abscissa (optional:also y1) 18! spline_1: return value of 1. derivative of spline for given abscissa 19! spline_int: return value of integral on given interval of spline 20! kill_spline: destructor ( spl%x,y und/oder spl%y2) 21 22! NB: splines are always "natural splines", i.e. values of first 23! derivative at the end-points cannot be specified 24!-----------------------------------------------------------------------------! 25 USE kinds, ONLY : DP 26 27 IMPLICIT NONE 28 29 PRIVATE 30 PUBLIC :: spline_data, allocate_spline, init_spline, spline, spline_1, & 31 spline_int, kill_spline, splineh, splinedx, splintdx, nullify_spline 32 33 TYPE spline_data 34 REAL (DP), POINTER :: x(:) ! array containing x values 35 REAL (DP), POINTER :: y(:) ! array containing y values 36 ! y(i) is the function value corresponding 37 ! to x(i) in the interpolation table 38 REAL (DP), POINTER :: y2(:) ! second derivative of interpolating function 39 INTEGER :: n ! number of element in the interpolation table 40 INTEGER :: pos 41 REAL (DP) :: h, invh, h26, h16 42 REAL (DP) :: xmin, xmax ! ... added by Carlo Cavazzoni 43 END TYPE spline_data 44 45!-----------------------------------------------------------------------------! 46 47 CONTAINS 48 49!-----------------------------------------------------------------------------! 50 51 SUBROUTINE nullify_spline( spl ) 52 TYPE (spline_data), INTENT (INOUT) :: spl 53 NULLIFY( spl%x ) 54 NULLIFY( spl%y ) 55 NULLIFY( spl%y2 ) 56 spl%n = 0 57 spl%pos = 0 58 spl%h = 0.0d0 59 spl%invh = 0.0d0 60 spl%h26 = 0.0d0 61 spl%h16 = 0.0d0 62 spl%xmin = 0.0d0 63 spl%xmax = 0.0d0 64 RETURN 65 END SUBROUTINE nullify_spline 66 67 SUBROUTINE allocate_spline( spl, nn, xmin, xmax ) 68 69 IMPLICIT NONE 70 71 TYPE (spline_data), INTENT (INOUT) :: spl 72 INTEGER, INTENT (IN) :: nn 73 REAL(DP), INTENT (IN), OPTIONAL :: xmin, xmax 74 75 INTEGER err 76 77 IF( PRESENT( xmin ) .AND. .NOT. PRESENT( xmax ) ) & 78 CALL errore(' allocate_spline ', ' wrong number of arguments ', 1 ) 79 80 spl%n = nn 81 82 IF ( associated(spl%x) ) THEN 83 DEALLOCATE (spl%x,STAT=err) 84 IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%x',1) 85 NULLIFY (spl%x) 86 END IF 87 88 ! note that spl%x is not allocated if we use a regular x grid 89 90 IF( PRESENT( xmin ) .AND. PRESENT( xmax ) ) THEN 91 IF( xmin >= xmax ) & 92 CALL errore(' allocate_spline ', ' wrong interval ', 1) 93 spl%xmin = xmin 94 spl%xmax = xmax 95 spl%h = ( xmax - xmin ) / DBLE( nn - 1 ) 96 spl%invh = 1.0d0 / spl%h 97 ELSE 98 spl%xmin = 0 99 spl%xmax = 0 100 ALLOCATE (spl%x(1:nn),STAT=err) 101 IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%x',1) 102 END IF 103 104 IF (associated(spl%y)) THEN 105 DEALLOCATE (spl%y,STAT=err) 106 IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%y',1) 107 NULLIFY (spl%y) 108 END IF 109 110 ALLOCATE (spl%y(1:nn),STAT=err) 111 IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%y',1) 112 113 IF (associated(spl%y2)) THEN 114 DEALLOCATE (spl%y2,STAT=err) 115 IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%y2',1) 116 NULLIFY (spl%y2) 117 END IF 118 119 ALLOCATE (spl%y2(1:nn),STAT=err) 120 IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%y2',1) 121 122 RETURN 123 END SUBROUTINE allocate_spline 124 125 126!----------------------------------------------------------------------- 127 128 SUBROUTINE init_spline( spl, endpt, y1a, y1b ) 129 130! endpt: 's': regular spacing 131! 'l': left; 'r': right; 'b': both = specify 1-deriv for each endpoints 132 133 IMPLICIT NONE 134 TYPE (spline_data), INTENT (INOUT) :: spl 135 CHARACTER (len=*), INTENT (IN), OPTIONAL :: endpt 136 REAL (DP), INTENT (IN), OPTIONAL :: y1a, y1b ! end point derivative 137 INTEGER :: err, i, k, n 138 REAL (DP) :: p, qn, sig, un, y1l, y1r, dyp, dym, dxp, dxm, dxpm 139 REAL (DP), POINTER :: ww(:) 140 CHARACTER (len=8) :: ep 141 LOGICAL :: reg, lep, rep 142 143 ! shortcat for regular mesh without table of x values 144 145 IF( .NOT. ASSOCIATED( spl%x ) ) THEN 146 CALL splinedx( spl%xmin, spl%xmax, spl%y(:), spl%n, 0.0d0, 0.0d0, spl%y2(:) ) 147 RETURN 148 END IF 149 150 ! Find out if y first derivative is given at endpoints 151 152 IF ( .NOT. present(endpt) ) THEN 153 ep = ' ' 154 ELSE 155 ep = endpt 156 END IF 157 reg = ( scan(ep,'sS') > 0 ) 158 lep = ( scan(ep,'lL') > 0 ) .OR. ( scan(ep,'bB') > 0 ) 159 rep = ( scan(ep,'rR') > 0 ) .OR. ( scan(ep,'bB') > 0 ) 160 161 ! check input parameter consistency 162 163 IF ( ( lep .OR. rep ) .AND. .NOT. present(y1a) ) & 164 CALL errore( 'init_spline', 'first deriv. at end-point missing', 1 ) 165 IF ( lep .AND. rep .AND. .NOT. present(y1b) ) & 166 CALL errore( 'init_spline', 'first deriv. at end-point missing', 1 ) 167 168 ! define endpoints derivative 169 170 IF ( lep ) y1l = y1a 171 IF ( rep .AND. .NOT. lep ) y1r = y1a 172 IF ( rep .AND. lep ) y1r = y1b 173 174 spl%pos = 1 175 ALLOCATE ( ww( 1 : spl%n ), STAT = err ) 176 IF (err/=0) CALL errore('init_spline','could not allocate ww',1) 177 178 n = spl % n 179 180 IF ( lep ) THEN 181 spl%y2(1) = -0.5d0 182 dxp = spl%x(2) - spl%x(1) 183 dyp = spl%y(2) - spl%y(1) 184 ww(1) = ( 3.0d0 / dxp ) * ( dyp / dxp - y1l ) 185 ELSE 186 spl%y2(1) = 0 187 ww(1) = 0.d0 188 END IF 189 190 DO i = 2, n - 1 191 192 dxp = spl%x(i+1) - spl%x(i) 193 dxm = spl%x(i) - spl%x(i-1) 194 dxpm = spl%x(i+1) - spl%x(i-1) 195 196 sig = dxm / dxpm 197 p = sig * spl%y2(i-1) + 2.0d0 198 spl%y2(i) = ( sig - 1.0d0 ) / p 199 200 dyp = spl%y(i+1) - spl%y(i) 201 dym = spl%y(i) - spl%y(i-1) 202 203 ww(i) = ( 6.0d0 * ( dyp / dxp - dym / dxm ) / dxpm - sig * ww(i-1) ) / p 204 205 END DO 206 207 IF ( rep ) THEN 208 qn = 0.5d0 209 dxm = spl%x(n) - spl%x(n-1) 210 dym = spl%y(n) - spl%y(n-1) 211 un = ( 3.0d0 / dxm ) * ( y1r - dym / dxm ) 212 ELSE 213 qn = 0 214 un = 0 215 END IF 216 217 spl % y2(n) = ( un - qn * ww(n-1) ) / ( qn * spl%y2(n-1) + 1.0d0 ) 218 219 DO k = n - 1, 1, -1 220 spl % y2(k) = spl%y2(k) * spl%y2(k+1) + ww(k) 221 END DO 222 223 DEALLOCATE ( ww, STAT = err ) 224 IF (err/=0) CALL errore('init_spline','could not deallocate ww',1) 225 226 IF ( reg ) THEN 227 spl%h = ( spl%x(n) - spl%x(1) ) / ( n - 1.0d0 ) 228 spl%h16 = spl%h / 6.0d0 229 spl%h26 = spl%h**2 / 6.0d0 230 spl%invh = 1 / spl%h 231 ELSE 232 spl%h = 0.0d0 233 spl%invh = 0.0d0 234 END IF 235 236 RETURN 237 END SUBROUTINE init_spline 238 239!----------------------------------------------------------------------- 240 241 FUNCTION interv( spl, xx ) 242 243 IMPLICIT NONE 244 245 TYPE (spline_data), INTENT (IN) :: spl 246 REAL (DP), INTENT (IN) :: xx 247 INTEGER :: interv 248 INTEGER :: khi, klo, i, p, n, k 249 250 ! if we have a regular mesh use a quick position search 251 252 IF ( spl%h /= 0 ) THEN 253 i = ( xx - spl%x(1) ) * spl%invh + 1 254 IF ( i < 1 .OR. i > spl%n ) & 255 CALL errore('interv', 'illegal x-value passed to interv',1) 256 interv = i 257 RETURN 258 END IF 259 260 p = spl%pos 261 IF ( p >= spl%n .OR. p <= 1 ) p = spl%n / 2 262 i = 0 263 n = spl%n 264 265 ! check if interval is close to previous interval 266 267 IF ( xx < spl%x(p+1) ) THEN 268 IF ( xx >= spl%x(p) ) THEN 269 i = spl%pos 270 ELSE IF ( p > 1 .AND. xx >= spl%x(p-1) ) THEN 271 i = p - 1 272 ELSE 273 klo = 1 274 khi = p + 1 275 END IF 276 ELSE IF ( (p + 2) <= n .AND. xx < spl%x(p+2) ) THEN 277 i = p + 1 278 ELSE 279 klo = p + 1 280 khi = n 281 END IF 282 283 ! perform binary search 284 285 IF ( i == 0 ) THEN 286 IF ( xx < spl%x(1) .OR. xx > spl%x(n) ) & 287 CALL errore('interv', 'xx value out of spline-range',1) 288 DO WHILE ( (khi - klo) > 1 ) 289 k = ( khi + klo ) / 2 290 IF ( spl%x(k) > xx ) THEN 291 khi = k 292 ELSE 293 klo = k 294 END IF 295 END DO 296 i = klo 297 END IF 298 299 interv = i 300 RETURN 301 END FUNCTION interv 302 303 304!----------------------------------------------------------------------- 305 FUNCTION spline( spl, xx, y1 ) 306 307 IMPLICIT NONE 308 309 TYPE (spline_data), INTENT (INOUT) :: spl 310 REAL (DP), INTENT (IN) :: xx 311 REAL (DP), INTENT (OUT), OPTIONAL :: y1 312 REAL (DP) :: spline 313 314 INTEGER :: khi, klo 315 REAL (DP) :: a, b, h, invh, ylo, yhi, y2lo, y2hi, a3ma, b3mb 316 317 ! shortcat for regular mesh without table of x values 318 319 IF( .NOT. ASSOCIATED( spl%x ) ) THEN 320 IF( PRESENT( y1 ) ) & 321 CALL errore(' spline ', ' y1 without x table not implemented ', 1 ) 322 CALL splintdx( spl%xmin, spl%xmax, spl%y, spl%y2, spl%n, xx, a ) 323 spline = a 324 RETURN 325 END IF 326 327 spl%pos = interv( spl, xx ) 328 klo = spl%pos 329 khi = spl%pos + 1 330 331 IF ( spl%h /= 0 ) THEN 332 h = spl%h 333 invh = spl%invh 334 ELSE 335 h = spl%x( khi ) - spl%x( klo ) 336 invh = spl%invh 337 IF ( h == 0.0d0 ) & 338 CALL errore('spline','bad spl%x input',1) 339 END IF 340 341 a = ( spl%x( khi ) - xx ) * invh 342 b = 1 - a 343 a3ma = a**3 - a 344 b3mb = b**3 - b 345 ylo = spl%y( klo ) 346 yhi = spl%y( khi ) 347 y2lo = spl%y2( klo ) 348 y2hi = spl%y2( khi ) 349 spline = a * ylo + b * yhi + ( a3ma * y2lo + b3mb * y2hi ) * ( h**2 ) / 6.0d0 350 351 IF ( present( y1 ) ) then 352 y1 = ( yhi - ylo ) * invh + & 353 ( ( 1.0d0 - 3 * a**2 ) * y2lo + ( 3 * b**2 - 1.0d0 ) * y2hi ) * h / 6.0d0 354 END IF 355 356 RETURN 357 END FUNCTION spline 358 359!----------------------------------------------------------------------- 360 361 FUNCTION splineh(spl,xx,y1) 362 IMPLICIT NONE 363 TYPE (spline_data), INTENT (IN) :: spl 364 REAL (DP), INTENT (IN) :: xx 365 REAL (DP), INTENT (OUT) :: y1 366 REAL (DP) :: splineh 367 368 INTEGER :: khi, klo, i 369 REAL (DP) :: a, b, h, invh, t, ylo, yhi, y2lo, y2hi, d, d0 370 371! fast spline for pair potentials without checks 372 h = spl%h 373 invh = spl%invh 374 d=xx-spl%x(1); i=INT(d*spl%invh); d0=DBLE(i)*h; i=i+1 375 i = (xx-spl%x(1))*invh + 1 376 377 a = (spl%x(i+1)-xx)*invh 378 b = 1 - a 379 t = -a*b 380! b=(d-d0)*invh; a=1-b; t=-a*b 381 ylo = spl%y(i) 382 yhi = spl%y(i+1) 383 y2lo = spl%y2(i) 384 y2hi = spl%y2(i+1) 385 splineh = a*ylo + b*yhi + ((a+1)*y2lo+(b+1)*y2hi)*t*spl%h26 386 y1 = (yhi-ylo)*invh + ((1.d0-3*a*a)*y2lo+(3*b*b-1.d0)*y2hi)*spl%h16 387 388 END FUNCTION splineh 389!----------------------------------------------------------------------- 390 FUNCTION spline_1(spl,xx) 391 IMPLICIT NONE 392 TYPE (spline_data), INTENT (INOUT) :: spl 393 REAL (DP), INTENT (IN) :: xx 394 REAL (DP) :: spline_1 395 396 INTEGER :: khi, klo 397 REAL (DP) :: a, b, h 398 399 spl%pos = interv(spl,xx) 400 klo = spl%pos 401 khi = spl%pos + 1 402 403 h = spl%x(khi) - spl%x(klo) 404 IF (h==0.d0) CALL errore('spline','bad spl%x input',1) 405 a = (spl%x(khi)-xx)/h 406 b = 1 - a 407 spline_1 = (spl%y(khi)-spl%y(klo))/h + ((1.d0-3*a**2)*spl%y2(klo)+(3*b** & 408 2-1.d0)*spl%y2(khi))*h/6.d0 409 410 RETURN 411 END FUNCTION spline_1 412 413 414 415!----------------------------------------------------------------------- 416 FUNCTION stamm(spl,p,x) 417 IMPLICIT NONE 418 TYPE (spline_data), INTENT (IN) :: spl 419 INTEGER, INTENT (IN) :: p 420 REAL (DP), INTENT (IN) :: x 421 REAL (DP) :: stamm 422 REAL (DP) :: a, b, aa, bb, h 423 424 h = spl%x(p+1) - spl%x(p) 425 b = (x-spl%x(p))/h 426 a = 1 - b 427 aa = a**2 428 bb = b**2 429 stamm = 0.5d0*h*(bb*spl%y(p+1)-aa*spl%y(p)) + h**3/12.d0*(aa*(1-0.5d0*aa)* & 430 spl%y2(p)-bb*(1-0.5d0*bb)*spl%y2(p+1)) 431 432 RETURN 433 END FUNCTION stamm 434 435 436 437!----------------------------------------------------------------------- 438 FUNCTION spline_int(spl,x0,x1) 439 IMPLICIT NONE 440 TYPE (spline_data), INTENT (INOUT) :: spl 441 REAL (DP), INTENT (IN) :: x0, x1 442 REAL (DP) :: spline_int 443 444 INTEGER :: j, pa, pb 445 REAL (DP) :: h, vorz, xa, xb, i1, i2 446 447 vorz = 1 448 xa = min(x0,x1) 449 xb = max(x0,x1) 450 IF (x0>x1) vorz = -1 451 IF (xa<spl%x(1) .OR. xb>spl%x(spl%n)) CALL errore('spline_int', & 452 'illegal integration range',1) 453 454 pa = interv(spl,xa) 455 pb = interv(spl,xb) 456 457 IF (pa==pb) THEN 458 spline_int = vorz*(stamm(spl,pa,xb)-stamm(spl,pa,xa)) 459 RETURN 460 END IF 461 462 i1 = 0 463 i2 = 0 464 DO j = pa + 1, pb - 1 465 h = spl%x(j+1) - spl%x(j) 466 i1 = i1 + h*(spl%y(j)+spl%y(j+1)) 467 i2 = i2 + h**3*(spl%y2(j)+spl%y2(j+1)) 468 END DO 469 h = spl%x(pa+1) - spl%x(pa) 470 i1 = i1 + h*spl%y(pa+1) 471 i2 = i2 + h**3*spl%y2(pa+1) 472 h = spl%x(pb+1) - spl%x(pb) 473 i1 = i1 + h*spl%y(pb) 474 i2 = i2 + h**3*spl%y2(pb) 475 476 spline_int = vorz*(i1/2.-i2/24.d0+stamm(spl,pb,xb)-stamm(spl,pa,xa)) 477 478 RETURN 479 END FUNCTION spline_int 480 481!----------------------------------------------------------------------- 482 SUBROUTINE kill_spline(spl,what) 483! deallocate splines 484! what=='a' or not present: deallocate all (spl%x, spl%y, spl%y2) 485! what=='d': deallocate only data (spl%x, spl%y) 486! what=='2': deallocate only table of 2. derivatives (spl%y2) 487 IMPLICIT NONE 488 TYPE (spline_data), INTENT (INOUT) :: spl 489 CHARACTER, INTENT (IN), OPTIONAL :: what 490 CHARACTER :: w 491 INTEGER :: err 492 493 w = 'a' 494 IF (present(what)) w = what 495 SELECT CASE (w) 496 CASE ('d','D') 497 IF (associated(spl%x)) THEN 498 DEALLOCATE (spl%x,STAT=err) 499 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%x',1) 500 NULLIFY (spl%x) 501 END IF 502 IF (associated(spl%y)) THEN 503 DEALLOCATE (spl%y,STAT=err) 504 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y',1) 505 NULLIFY (spl%y) 506 END IF 507 CASE ('2') 508 IF (associated(spl%y2)) THEN 509 DEALLOCATE (spl%y2,STAT=err) 510 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y2',1) 511 NULLIFY (spl%y2) 512 END IF 513 CASE ('a','A') 514 IF (associated(spl%x)) THEN 515 DEALLOCATE (spl%x,STAT=err) 516 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%x',1) 517 NULLIFY (spl%x) 518 END IF 519 IF (associated(spl%y)) THEN 520 DEALLOCATE (spl%y,STAT=err) 521 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y',1) 522 NULLIFY (spl%y) 523 END IF 524 IF (associated(spl%y2)) THEN 525 DEALLOCATE (spl%y2,STAT=err) 526 IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y2',1) 527 NULLIFY (spl%y2) 528 END IF 529 END SELECT 530 531 RETURN 532 END SUBROUTINE kill_spline 533 534!=-----------------------------------------------------------------------=! 535! Subroutines: splinedx, splintdx 536! added for compatibility with SISSA code 537! Carlo Cavazzoni 15-03-2000 538!=-----------------------------------------------------------------------=! 539 540 SUBROUTINE splinedx(xmin,xmax,y,n,yp1,ypn,y2) 541 USE kinds 542 IMPLICIT NONE 543 INTEGER, INTENT(IN) :: n 544 REAL(DP), INTENT(IN) :: yp1,ypn,xmin,xmax,y(:) 545 REAL(DP), INTENT(OUT) :: y2(:) 546 INTEGER :: i, k 547 REAL(DP) :: p, qn, sig, un, dx 548 REAL(DP) :: u(n) 549 550 dx = (xmax-xmin)/DBLE(n-1) 551 if ( yp1 .gt. 0.99d30 ) then 552 y2(1)=0.d0 553 u(1)=0.0d0 554 else 555 y2(1)=-0.5d0 556 u(1)=(3.d0/dx) * ( (y(2)-y(1)) / dx - yp1 ) 557 endif 558 do i=2,n-1 559 sig=0.5d0 560 p=sig*y2(i-1)+2.d0 561 y2(i)=(sig-1.d0)/p 562 u(i) = (6.0d0 * ( (y(i+1)-y(i))/ dx - (y(i)-y(i-1))/ dx ) & 563 / (2.0d0*dx) - sig * u(i-1) ) / p 564 end do 565 if ( ypn .gt. 0.99d30 ) then 566 qn=0.d0 567 un=0.d0 568 else 569 qn=0.5d0 570 un= ( 3.d0 / dx ) * ( ypn - (y(n)-y(n-1)) / dx ) 571 endif 572 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0d0) 573 do k=n-1,1,-1 574 y2(k)=y2(k)*y2(k+1)+u(k) 575 end do 576 return 577 END SUBROUTINE splinedx 578 579 580 SUBROUTINE splintdx(xmin,xmax,ya,y2a,n,x,y) 581 USE kinds 582 IMPLICIT NONE 583 INTEGER, INTENT(IN) :: n 584 REAL(DP), INTENT(IN) :: x,xmin,xmax,ya(:),y2a(:) 585 REAL(DP), INTENT(OUT) :: y 586 INTEGER :: khi,klo 587 REAL(DP) :: a,b,h,dx,xhi,xlo 588 dx = (xmax-xmin)/DBLE(n-1) 589 klo = INT(x/dx) 590 khi = klo+1 591 IF(klo.LT.1) THEN 592 CALL errore(' splintdx ',' klo less than one ',klo) 593 END IF 594 IF(khi.GT.n) THEN 595 CALL errore(' splintdx ',' khi grether than N ',khi) 596 END IF 597 xlo = xmin + DBLE(klo-1) * dx 598 xhi = xmin + DBLE(khi-1) * dx 599 600 a = (xhi-x)/dx 601 b = (x-xlo)/dx 602 y = a*ya(klo) + b*ya(khi) + & 603 ( (a*a*a-a)*y2a(klo) + (b*b*b-b)*y2a(khi) ) * (dx*dx)/6.0d0 604 RETURN 605 END SUBROUTINE splintdx 606 607!----------------------------------------------------------------------- 608 609 SUBROUTINE nr_spline( x, y, n, yp1, ypn, y2 ) 610 INTEGER :: n 611 REAL(DP) :: yp1, ypn, x(n), y(n), y2(n) 612 INTEGER :: i, k 613 REAL(DP) :: p, qn, sig, un 614 REAL(DP) :: u( n ) 615 if ( yp1 .gt. 0.99d30 ) then 616 y2(1)=0.d0 617 u(1)=0.d0 618 else 619 y2(1)=-0.5d0 620 u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) 621 endif 622 do i=2,n-1 623 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) 624 p=sig*y2(i-1)+2.d0 625 y2(i)=(sig-1.d0)/p 626 u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / & 627 (x(i)-x(i-1))) / (x(i+1)-x(i-1))-sig*u(i-1))/p 628 end do 629 if ( ypn .gt. 0.99d30 ) then 630 qn=0.d0 631 un=0.d0 632 else 633 qn=0.5d0 634 un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) 635 endif 636 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0) 637 do k=n-1,1,-1 638 y2(k)=y2(k)*y2(k+1)+u(k) 639 end do 640 return 641 END SUBROUTINE nr_spline 642 643 644 SUBROUTINE nr_splint( xa, ya, y2a, n, x, y ) 645 INTEGER :: n 646 REAL(DP) :: x,y,xa(n),y2a(n),ya(n) 647 INTEGER :: k,khi,klo 648 REAL(DP) :: a,b,h 649 klo=1 650 khi=n 6511 if (khi-klo.gt.1) then 652 k=(khi+klo)/2 653 if(xa(k).gt.x)then 654 khi=k 655 else 656 klo=k 657 endif 658 goto 1 659 endif 660 h=xa(khi)-xa(klo) 661 if (h.eq.0.) & 662 CALL errore(' nr_splint ', 'bad xa input in splint' , 1 ) 663 a=(xa(khi)-x)/h 664 b=(x-xa(klo))/h 665 y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))* & 666 (h**2)/6.d0 667 return 668 END SUBROUTINE nr_splint 669 670 671 SUBROUTINE nr_splie2( x1a, x2a, ya, m, n, y2a ) 672 INTEGER :: m, n 673 REAL(DP) :: x1a(m), x2a(n), y2a(m,n), ya(m,n) 674 INTEGER :: j,k 675 REAL(DP) :: y2tmp(n), ytmp(n) 676 do j = 1, m 677 do k = 1, n 678 ytmp(k) = ya(j,k) 679 end do 680 call nr_spline( x2a, ytmp, n, 1.d30, 1.d30, y2tmp ) 681 do k = 1, n 682 y2a(j,k) = y2tmp(k) 683 end do 684 end do 685 return 686 END SUBROUTINE nr_splie2 687 688 689 SUBROUTINE nr_splin2( x1a, x2a, ya, y2a, m, n, x1, x2, y ) 690 INTEGER :: m, n 691 REAL(DP) :: x1, x2, y, x1a(m), x2a(n), y2a(m,n), ya(m,n) 692 INTEGER :: j, k 693 REAL(DP) :: y2tmp( MAX(n,m) ), ytmp( n ), yytmp( MAX(n,m) ) 694 do j = 1, m 695 do k = 1, n 696 ytmp(k) = ya(j,k) 697 y2tmp(k) = y2a(j,k) 698 end do 699 call nr_splint( x2a, ytmp, y2tmp, n, x2, yytmp(j) ) 700 end do 701 call nr_spline( x1a, yytmp, m, 1.d30, 1.d30, y2tmp ) 702 call nr_splint( x1a, yytmp, y2tmp, m, x1, y ) 703 return 704 END SUBROUTINE nr_splin2 705 706!----------------------------------------------------------------------- 707 708 END MODULE splines 709