1 Subroutine spline_real (n, x, y, r) 2 implicit none 3#include "dimqm_constants.fh" 4! ============================================================================ 5! Given arrays x(1:n) and y(1:n) containing a tabulated function, 6! i.e. y_i = f(x_i), with x_1 < x_2 ... < x_n, and given values yp1 and ypn 7! for the first derivative of the interpolation function at point 1 and n, 8! respectively, this routine returns an array r of length N that contains the 9! second derivatives of the interpolating function at the tabulated points x. 10! The interpolation is "natural", meaning that the second derivative on the 11! endpoints is zero. 12! Adopted from Numerical Recipes by Seth M. Morton, 2011 13! ============================================================================ 14 15 integer, Intent(in) :: n 16 double precision, Intent(In) :: x(n) ! Domain of data 17 double precision, Intent(In) :: y(n) ! Range of rata (i.e. f(x)) 18 double precision, Intent(Out) :: r(n) ! The second derivatives at data points 19 20 double precision :: a(n) 21 double precision :: b(n) 22 double precision :: c(n) 23 24 integer :: ier 25 26! Set up the tridiagonal equations 27 c(1:n-1) = x(2:n) - x(1:n-1) 28 r(1:n-1) = SIX * ( ( y(2:n) - y(1:n-1) ) / c(1:n-1) ) 29 r(2:n-1) = r(2:n-1) - r(1:n-1) 30 a(2:n-1) = c(1:n-2) 31 b(2:n-1) = TWO * ( c(2:n-1) + a(2:n-1) ) 32 b(1) = ONE 33 b(n) = ONE 34 35! --------------------------------------------------- 36! The lower and upper boundaries are set to "natural" 37! --------------------------------------------------- 38 39 r(1) = ZERO 40 c(1) = ZERO 41 r(n) = ZERO 42 c(n) = ZERO 43 44! --------------------------------------------------------------- 45! Solve the tridiagonal system of equations with a LAPACK routine 46! --------------------------------------------------------------- 47 48 call dgtsv (n, 1, a(2:n), b(1:n), c(1:n-1), r(1:n), n, ier) 49 if (ier /= 0) then 50 write(*,*) 'Error in spline' 51 stop 52 end if 53 54 End Subroutine spline_real 55 56 Subroutine spline (n, x, yr, yi, rr, ri) 57 implicit none 58 59! ========================================================= 60! An extention of spline_real to spline complex functions 61! Two real arrays representing real and imaginary are given 62! ========================================================= 63 64 Integer, Intent(in) :: n 65 double precision, Intent(In) :: x(n) ! Domain of data 66 double precision, Intent(In) :: yr(n) ! Range of rata (i.e. f(x)) (real) 67 double precision, Intent(In) :: yi(n) ! Range of rata (i.e. f(x)) (imag) 68 double precision, Intent(Out) :: rr(n) ! The second derivatives at data points 69 double precision, Intent(Out) :: ri(n) ! The second derivatives at data points 70 71! ------------------------------------ 72! Spline real and imaginary separately 73! ------------------------------------ 74 75 call spline_real (n, x, yr, rr) 76 call spline_real (n, x, yi, ri) 77 78 End Subroutine spline 79 80 Function interpolate_real (n, xa, ya, y2a, x) Result(y) 81 implicit none 82#include "dimqm_constants.fh" 83 84! ============================================================================= 85! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulates a function 86! (with the xa_i 's in order) and given the array y2a(1:n) which is output from 87! spline (above), and given a value of x, this routine returns a cubic-spline 88! interpolation value y. 89! iter was added to prevent a new search for each iteration in a sequencial 90! series of calls. Use iter=1 for each call if they are unrelated from each 91! other. 92! Adopted from Numerical Recipes by Seth M. Morton, 2011 93! ============================================================================= 94 95 integer, Intent(In) :: n 96 double precision, Intent(In) :: xa(n) 97 double precision, Intent(In) :: ya(n) 98 double precision, Intent(In) :: y2a(n) 99 double precision, Intent(In) :: x 100 double precision :: y 101 102 integer :: klo 103 integer :: khi 104 double precision :: a 105 double precision :: b 106 double precision :: h 107 integer locate 108 109! ----------------------------- 110! Find low and high index range 111! ----------------------------- 112 113 klo = MAX(MIN(locate(n, xa, x), n-1), 1) 114 khi = klo + 1 115 116 if (xa(klo) == xa(khi)) then 117 write(*,*) 'bad xa input in interpolate' 118 stop 119 end if 120 121! ------------------- 122! Evaluate the spline 123! ------------------- 124 125 h = xa(khi) - xa(klo) 126 a = ( xa(khi) - x ) / h 127 b = ( x - xa(klo) ) / h 128 y = a * ya(klo) + b * ya(khi) 129 $ + ( ( a**3 - a ) * y2a(klo) + ( b**3 - b ) * y2a(khi) ) 130 $ * ( h**2 ) / SIX 131 132 End Function interpolate_real 133 134 Function interpolate_complex1 (n, xa, yar, yai, y2ar, y2ai, x) 135 $ Result(y) 136 implicit none 137#include "dimqm_constants.fh" 138 139! =============================================================== 140! Extention of interpolate_real to interpolate a complex function 141! The original y-values are given in two separate real arrays 142! representing real and imaginary. 143! =============================================================== 144 145 integer, Intent(In) :: n 146 double precision, Intent(In) :: xa(n) 147 double precision, Intent(In) :: yar(n) 148 double precision, Intent(In) :: yai(n) 149 double precision, Intent(In) :: y2ar(n) 150 double precision, Intent(In) :: y2ai(n) 151 double precision, Intent(In) :: x 152 double complex :: y 153 154 integer :: klo 155 integer :: khi 156 double precision :: a 157 double precision :: b 158 double precision :: h 159 double precision :: t 160 integer locate 161 162! ----------------------------- 163! Find low and high index range 164! -----------------------------` 165 166 klo = MAX(MIN(locate(n, xa, x), n-1), 1) 167 khi = klo + 1 168 169 if (xa(klo) == xa(khi)) then 170 write(*,*) 'bad xa input in interpolate' 171 stop 172 end if 173 174! ------------------- 175! Evaluate the spline 176! ------------------- 177 178 h = xa(khi) - xa(klo) 179 a = ( xa(khi) - x ) / h 180 b = ( x - xa(klo) ) / h 181 y = a * yar(klo) + b * yar(khi) 182 $ + ( ( a**3 - a ) * y2ar(klo) + ( b**3 - b ) * y2ar(khi) ) 183 $ * ( h**2 ) / SIX 184 t = a * yai(klo) + b * yai(khi) 185 $ + ( ( a**3 - a ) * y2ai(klo) + ( b**3 - b ) * y2ai(khi) ) 186 $ * ( h**2 ) / SIX 187 y = CMPLX(REAL(y), t) 188 189 End Function interpolate_complex1 190 191 Function interpolate_complex2 (n, xa, ya, y2a, x) Result(y) 192 implicit none 193#include "dimqm_constants.fh" 194 195! =============================================================== 196! Extention of interpolate_real to interpolate a complex function 197! The original y-values are given as complex. 198! =============================================================== 199 200 integer, Intent(In) :: n 201 double precision, Intent(In) :: xa(n) 202 double complex, Intent(In) :: ya(n) 203 double complex, Intent(In) :: y2a(n) 204 double precision, Intent(In) :: x 205 double complex :: y 206 207 integer :: klo 208 integer :: khi 209 double precision :: a 210 double precision :: b 211 double precision :: h 212 integer locate 213 214! ----------------------------- 215! Find low and high index range 216! ----------------------------- 217 218 klo = MAX(MIN(locate(n, xa, x), n-1), 1) 219 khi = klo + 1 220 221 if (xa(klo) == xa(khi)) then 222 write(*,*) 'bad xa input in interpolate' 223 stop 224 end if 225 226! ------------------- 227! Evaluate the spline 228! ------------------- 229 230 h = xa(khi) - xa(klo) 231 a = ( xa(khi) - x ) / h 232 b = ( x - xa(klo) ) / h 233 y = a * ya(klo) + b * ya(khi) 234 $ + ( ( a**3 - a ) * y2a(klo) + ( b**3 - b ) * y2a(khi) ) 235 $ * ( h**2 ) / SIX_C 236 237 End Function interpolate_complex2 238 239 Function locate (n, xx, x) Result(i) 240 implicit none 241#include "dimqm_constants.fh" 242 243! ========================================================================== 244! Given an array xx(1:n), and given a value x, returns a value j such that 245! x is between xx(j) and xx(j+1). xx must be monotonic, either increasing or 246! decreasing. j =0 or j =n is returned to indicate that x is out of range. 247! Adopted from Numerical Recipes by Seth M. Morton, 2011 248! ========================================================================== 249 250 integer, Intent(In) :: n 251 double precision, Intent(In) :: xx(n) 252 double precision, Intent(In) :: x 253 integer :: i 254 255 integer :: jl 256 integer :: jm 257 integer :: ju 258 Logical :: ascnd 259 260! ------------------------ 261! Quick return if possible 262! ------------------------ 263 264 if (x == xx(1)) then 265 i = 1 266 return 267 else if (x == xx(n)) then 268 i = n - 1 269 return 270 end if 271 272! --------------------------------- 273! Initiallize limits and order flag 274! --------------------------------- 275 276! True if ascending order of table, false otherwise. 277 ascnd = xx(n) >= xx(1) 278 jl = 0 279 ju = n + 1 280 281! ----------------------------------------- 282! Find the index for the value immeadiately 283! below or equal to the requested value 284! ----------------------------------------- 285 286! Repeat until this condition is satisfied. 287 do while (ju - jl > 1) 288! Compute a midpoint, and replace either the lower limit 289! or the upper limit, as appropriate. 290 jm = ( ju + jl ) / 2 291 if (ascnd .eqv. ( x >= xx(jm) )) then 292 jl = jm 293 else 294 ju = jm 295 end if 296 end do 297 298! Set the output 299 i = jl 300 301 End Function locate 302