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