Branch data Line data Source code
1 : : !!****m* ABINIT/m_splines
2 : : !! NAME
3 : : !! m_splines
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains routines for spline interpolation.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2010-2026 ABINIT group (YP, BAmadon)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_splines
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : :
28 : : use m_fstrings, only : sjoin, itoa, ftoa
29 : : !use m_time, only : timab
30 : :
31 : : implicit none
32 : :
33 : : public :: splfit
34 : : public :: spline
35 : : public :: spline_bicubic
36 : : public :: spline_c
37 : : public :: spline_complex
38 : : public :: spline_integrate
39 : : public :: splint
40 : : public :: splint_complex
41 : :
42 : : !FIXME deprecated
43 : : public :: intrpl
44 : :
45 : : ! *************************************************************************
46 : :
47 : : contains
48 : : !!***
49 : :
50 : : !----------------------------------------------------------------------
51 : :
52 : : !!****f* m_splines/splfit
53 : : !! NAME
54 : : !! splfit
55 : : !!
56 : : !! FUNCTION
57 : : !! Evaluate cubic spline fit to get function values on input set of ORDERED, UNFORMLY SPACED points.
58 : : !! Optionally gives derivatives (first and second) at those points too.
59 : : !! If point lies outside the range of arg, assign the extremal
60 : : !! point values to these points, and zero derivative.
61 : : !!
62 : : !! INPUTS
63 : : !! arg(numarg)=equally spaced arguments (spacing de) for data to which spline was fit.
64 : : !! fun(numarg,2)=function values to which spline was fit and spline
65 : : !! fit to second derivatives (from Numerical Recipes spline).
66 : : !! ider= see above
67 : : !! newarg(numnew)=new values of arguments at which function is desired.
68 : : !! numarg=number of arguments at which spline was fit.
69 : : !! numnew=number of arguments at which function values are desired.
70 : : !!
71 : : !! OUTPUT
72 : : !! derfun(numnew)=(optional) values of first or second derivative of function.
73 : : !! This is only computed for ider=1 or 2; otherwise derfun not used.
74 : : !! newfun(numnew)=values of function at newarg(numnew).
75 : : !! This is only computed for ider=0 or 1.
76 : : !!
77 : : !! NOTES
78 : : !! if ider=0, compute only the function (contained in fun)
79 : : !! if ider=1, compute the function (contained in fun) and its first derivative (in derfun)
80 : : !! if ider=2, compute only the second derivative of the function (in derfun)
81 : : !!
82 : : !! SOURCE
83 : :
84 : 17435364 : subroutine splfit(arg, derfun, fun, ider, newarg, newfun, numarg, numnew)
85 : :
86 : : integer, intent(in) :: ider, numarg, numnew
87 : : real(dp), intent(in) :: arg(numarg), fun(numarg,2), newarg(numnew)
88 : : real(dp), intent(out) :: derfun(numnew)
89 : : real(dp), intent(inout) :: newfun(numnew)
90 : :
91 : : !Local variables---------------------------------------
92 : : integer :: i,jspl
93 : : real(dp) :: argmin,de,d,aa,bb,cc,dd,de2_dby_six,de_dby_six
94 : : !real(dp) :: tsec(2)
95 : :
96 : : ! *************************************************************************
97 : :
98 : : ! argmin is smallest x value in spline fit; de is uniform spacing of spline argument
99 : 17435364 : argmin = arg(1)
100 : 17435364 : de = (arg(numarg) - argmin) / dble(numarg-1)
101 : 17435364 : de2_dby_six = de**2 / six
102 : 17435364 : de_dby_six = de / six
103 : :
104 [ - + ]: 17435364 : if (de < tol12) then
105 : 0 : ABI_ERROR(sjoin('spacing should be strictly positive, while de is: ', ftoa(de)))
106 : : endif
107 : :
108 : 17435364 : jspl = -1
109 : :
110 : : ! Do one loop for no grads, other for grads
111 : 17435364 : select case (ider)
112 : : case (0)
113 : :
114 : : ! Spline index loop for no grads:
115 [ + + ]: 2748730710 : do i=1,numnew
116 [ + + ]: 2748730710 : if (newarg(i) >= arg(numarg)) then
117 : : ! function values are being requested outside range of data.',a1,'
118 : : ! Function and slope will be set to values at upper end of data.
119 : :
120 : 42013021 : newfun(i) = fun(numarg,1)
121 : :
122 [ + + ]: 2695638205 : else if (newarg(i) <= arg(1)) then
123 : 60669662 : newfun(i) = fun(1,1)
124 : :
125 : : else
126 : 2634968543 : jspl = 1 + int((newarg(i) - argmin)/de)
127 : 2634968543 : d = newarg(i) - arg(jspl)
128 : 2634968543 : bb = d / de
129 : 2634968543 : aa = one - bb
130 : 2634968543 : cc = aa*(aa**2 -one) * de2_dby_six
131 : 2634968543 : dd = bb*(bb**2 -one) * de2_dby_six
132 : 2634968543 : newfun(i)= aa * fun(jspl,1) + bb*fun(jspl+1,1) + cc*fun(jspl,2) + dd*fun(jspl+1,2)
133 : : end if
134 : : enddo
135 : :
136 : : case (1)
137 : :
138 : : ! Spline index loop includes grads:
139 [ + + ]: 427378429 : do i=1,numnew
140 : :
141 [ + + ]: 427378429 : if (newarg(i) >= arg(numarg)) then
142 : 29125942 : newfun(i) = fun(numarg,1)
143 : 29125942 : derfun(i) = zero
144 : :
145 [ + + ]: 393053047 : else if (newarg(i) <= arg(1)) then
146 : 32536578 : newfun(i) = fun(1,1)
147 : 32536578 : derfun(i) = zero
148 : :
149 : : else
150 : : ! cubic spline interpolation:
151 : 360516469 : jspl = 1 + int((newarg(i) - arg(1)) / de)
152 : 360516469 : d = newarg(i) - arg(jspl)
153 : 360516469 : bb = d / de
154 : 360516469 : aa = one - bb
155 : 360516469 : cc = aa*(aa**2 - one) * de2_dby_six
156 : 360516469 : dd = bb*(bb**2 - one) * de2_dby_six
157 : 360516469 : newfun(i) = aa*fun(jspl,1) + bb*fun(jspl+1,1) + cc*fun(jspl,2) + dd*fun(jspl+1,2)
158 : : ! spline fit to first derivative:
159 : : ! note correction of Numerical Recipes sign error
160 : : derfun(i) = (fun(jspl+1,1)-fun(jspl,1)) / de + &
161 : 360516469 : (-(3.d0*aa**2 -one) * fun(jspl,2) + (3.d0*bb**2 -one) * fun(jspl+1,2)) * de_dby_six
162 : :
163 : : end if
164 : : enddo
165 : :
166 : : case (2)
167 : :
168 [ + + ]: 26829454 : do i=1,numnew
169 : :
170 [ - + ]: 26829454 : if (newarg(i) >= arg(numarg)) then
171 : 0 : derfun(i) = zero
172 : :
173 [ + + ]: 25673014 : else if (newarg(i) <= arg(1)) then
174 : 98 : derfun(i) = zero
175 : :
176 : : else
177 : : ! cubic spline interpolation:
178 : 25672916 : jspl = 1 + int((newarg(i) - argmin) / de)
179 : 25672916 : d = newarg(i) - arg(jspl)
180 : 25672916 : bb = d / de
181 : 25672916 : aa = one - bb
182 : : ! second derivative of spline (piecewise linear function)
183 : 25672916 : derfun(i) = aa*fun(jspl,2) + bb*fun(jspl+1,2)
184 : :
185 : : end if
186 : : enddo
187 : :
188 : : case default
189 [ + + + - ]: 17435364 : ABI_ERROR(sjoin("Invalid ider:", itoa(ider)))
190 : : end select
191 : :
192 : 17435364 : end subroutine splfit
193 : : !!***
194 : :
195 : : !----------------------------------------------------------------------
196 : :
197 : : !!****f* m_splines/spline
198 : : !! NAME
199 : : !! spline
200 : : !!
201 : : !! FUNCTION
202 : : !! SPLINE (originally SPLINE_CUBIC_SET) computes the second derivatives
203 : : !! of a cubic spline.
204 : : !!
205 : : !! INPUTS
206 : : !! Input, integer N, the number of data points; N must be at least 2.
207 : : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
208 : : !! spline will actually be linear.
209 : : !!
210 : : !! Input, double precision T(N), the knot values, that is, the points where data
211 : : !! is specified. The knot values should be distinct, and increasing.
212 : : !!
213 : : !! Input, double precision Y(N), the data values to be interpolated.
214 : : !!
215 : : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
216 : : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
217 : : !!
218 : : !! OUTPUT
219 : : !! Output, double precision YPP(N), the second derivatives of the cubic spline.
220 : : !! Work space, double precision DIAG(N) - should be removed ...
221 : : !!
222 : : !! SOURCE
223 : :
224 : 4449344 : subroutine spline( t, y, n, ybcbeg, ybcend, ypp )
225 : :
226 : : !*******************************************************************************
227 : : !
228 : : ! Discussion:
229 : : !
230 : : ! For data interpolation, the user must call SPLINE_CUBIC_SET to
231 : : ! determine the second derivative data, passing in the data to be
232 : : ! interpolated, and the desired boundary conditions.
233 : : !
234 : : ! The data to be interpolated, plus the SPLINE_CUBIC_SET output,
235 : : ! defines the spline. The user may then call SPLINE_CUBIC_VAL to
236 : : ! evaluate the spline at any point.
237 : : !
238 : : ! The cubic spline is a piecewise cubic polynomial. The intervals
239 : : ! are determined by the "knots" or abscissas of the data to be
240 : : ! interpolated. The cubic spline has continous first and second
241 : : ! derivatives over the entire interval of interpolation.
242 : : !
243 : : ! For any point T in the interval T(IVAL), T(IVAL+1), the form of
244 : : ! the spline is
245 : : !
246 : : ! SPL(T) = A(IVAL)
247 : : ! + B(IVAL) * ( T - T(IVAL) )
248 : : ! + C(IVAL) * ( T - T(IVAL) )**2
249 : : ! + D(IVAL) * ( T - T(IVAL) )**3
250 : : !
251 : : ! If we assume that we know the values Y(*) and YPP(*), which represent
252 : : ! the values and second derivatives of the spline at each knot, then
253 : : ! the coefficients can be computed as:
254 : : !
255 : : ! A(IVAL) = Y(IVAL)
256 : : ! B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
257 : : ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
258 : : ! C(IVAL) = YPP(IVAL) / 2
259 : : ! D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
260 : : !
261 : : ! Since the first derivative of the spline is
262 : : !
263 : : ! SPL'(T) = B(IVAL)
264 : : ! + 2 * C(IVAL) * ( T - T(IVAL) )
265 : : ! + 3 * D(IVAL) * ( T - T(IVAL) )**2,
266 : : !
267 : : ! the requirement that the first derivative be continuous at interior
268 : : ! knot I results in a total of N-2 equations, of the form:
269 : : !
270 : : ! B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
271 : : ! + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
272 : : !
273 : : ! or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
274 : : !
275 : : ! ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
276 : : ! - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
277 : : ! + YPP(IVAL-1) * H(IVAL-1)
278 : : ! + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
279 : : ! =
280 : : ! ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
281 : : ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
282 : : !
283 : : ! or
284 : : !
285 : : ! YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
286 : : ! + YPP(IVAL) * H(IVAL)
287 : : ! =
288 : : ! 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
289 : : ! - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
290 : : !
291 : : ! Boundary conditions must be applied at the first and last knots.
292 : : ! The resulting tridiagonal system can be solved for the YPP values.
293 : : !
294 : : ! Modified:
295 : : !
296 : : ! 07 February 1999
297 : : ! 28 November 2004 XGonze : double precision
298 : : ! make arguments similar to the Numeric Recipes routine
299 : : ! also use algorithmics similar to the Numeric Recipes routine
300 : : !
301 : : ! Author:
302 : : !
303 : : ! John Burkardt
304 : : ! (XGonze got it from http://www.psc.edu/~burkardt/src/spline/spline.html)
305 : : !
306 : : ! Parameters:
307 : : !
308 : : ! Input, integer N, the number of data points; N must be at least 2.
309 : : ! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
310 : : ! spline will actually be linear.
311 : : !
312 : : ! Input, double precision T(N), the knot values, that is, the points where data
313 : : ! is specified. The knot values should be distinct, and increasing.
314 : : !
315 : : ! Input, double precision Y(N), the data values to be interpolated.
316 : : !
317 : : ! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
318 : : ! conditions if IBCBEG or IBCEND is equal to 1 or 2.
319 : : !
320 : : ! Output, double precision YPP(N), the second derivatives of the cubic spline.
321 : : !
322 : : ! Work space, double precision DIAG(N) - should be removed ...
323 : : !
324 : : !
325 : : ! XG041127 : In the initial implementation, one had the control on
326 : : ! IBCBEG and IBCEND. Now, they are determined by the values
327 : : ! of YBCBEG, YBCEND. Option 2 has been disabled.
328 : : !
329 : : ! Input, integer IBCBEG, left boundary condition flag:
330 : : !
331 : : ! 0: the spline should be a quadratic over the first interval;
332 : : ! 1: the first derivative at the left endpoint should be YBCBEG;
333 : : ! 2: the second derivative at the left endpoint should be YBCBEG.
334 : : !
335 : : ! Input, integer IBCEND, right boundary condition flag:
336 : : !
337 : : ! 0: the spline should be a quadratic over the last interval;
338 : : ! 1: the first derivative at the right endpoint should be YBCEND;
339 : : ! 2: the second derivative at the right endpoint should be YBCEND.
340 : :
341 : : integer, intent(in) :: n
342 : : real(dp), intent(in) :: t(n)
343 : : real(dp), intent(in) :: y(n)
344 : : real(dp), intent(in) :: ybcbeg
345 : : real(dp), intent(in) :: ybcend
346 : :
347 : : real(dp), intent(out) :: ypp(n)
348 : :
349 : : integer :: ibcbeg
350 : : integer :: ibcend
351 : : integer :: i,k
352 : : real(dp) :: ratio,pinv
353 : 4449344 : real(dp), allocatable :: tmp(:)
354 : : !
355 : : ! Check.
356 : : !
357 [ - + ]: 4449344 : if ( n <= 1 ) then
358 : 0 : write(std_out,* ) ' '
359 : 0 : write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
360 : 0 : write(std_out,* ) ' The number of knots must be at least 2.'
361 : 0 : write(std_out,* ) ' The input value of N = ', n
362 : 0 : ABI_ERROR("Fatal error")
363 : : end if
364 : :
365 [ - + + - : 13348032 : ABI_MALLOC(tmp,(n))
+ - + - -
+ - + ]
366 : :
367 [ + + ]: 117484941 : do i = 1, n-1
368 [ - + ]: 117484941 : if ( t(i) >= t(i+1) ) then
369 : 0 : write(std_out,* ) ' '
370 : 0 : write(std_out,* ) 'SPLINE_CUBIC_SET - Fatal error!'
371 : 0 : write(std_out,* ) ' The knots must be strictly increasing, but'
372 : 0 : write(std_out,* ) ' T(', i,') = ', t(i)
373 : 0 : write(std_out,* ) ' T(',i+1,') = ', t(i+1)
374 : 0 : ABI_ERROR("Fatal error")
375 : : end if
376 : : end do
377 : : !
378 : : ! XG041127
379 : 4449344 : ibcbeg=1 ; ibcend=1
380 [ + + ]: 4449344 : if(ybcbeg>1.0d+30)ibcbeg=0
381 [ + + ]: 4449344 : if(ybcend>1.0d+30)ibcend=0
382 : : !
383 : : ! Set the first and last equations.
384 : : !
385 [ + + ]: 4449344 : if ( ibcbeg == 0 ) then
386 : 38916 : ypp(1) = 0.d0
387 : 38916 : tmp(1) = 0.d0
388 : : else if ( ibcbeg == 1 ) then
389 : 4410428 : ypp(1) = -0.5d0
390 : 4410428 : tmp(1) = (3.d0/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
391 : : end if
392 [ + + ]: 4449344 : if ( ibcend == 0 ) then
393 : 38916 : ypp(n) = 0.d0
394 : 38916 : tmp(n) = 0.d0
395 : : else if ( ibcend == 1 ) then
396 : 4410428 : ypp(n) = 0.5d0
397 : 4410428 : tmp(n) = (3.d0/(t(n)-t(n-1)))*(ybcend-(y(n)-y(n-1))/(t(n)-t(n-1)))
398 : : end if
399 : :
400 : : !
401 : : ! Set the intermediate equations.
402 : : !
403 [ + + ]: 113035597 : do i=2,n-1
404 : 108586253 : ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
405 : 108586253 : pinv = 1.0d0/(ratio*ypp(i-1) + 2.0d0)
406 : 108586253 : ypp(i) = (ratio-1.0d0)*pinv
407 : : tmp(i)=(6.0d0*((y(i+1)-y(i))/(t(i+1)-t(i))-(y(i)-y(i-1)) &
408 : 108586253 : & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
409 [ + + ]: 113035597 : if (abs(tmp(i))<1.d5*tiny(0.d0)) tmp(i)=0.d0 !MT20050927
410 : : enddo
411 : :
412 : : ! Solve the equations
413 : 4449344 : ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0d0)
414 [ + + ]: 117484941 : do k=n-1,1,-1
415 : 117484941 : ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
416 : : enddo
417 : :
418 : 4449344 : ABI_FREE(tmp)
419 : 4449344 : end subroutine spline
420 : : !!***
421 : :
422 : : !----------------------------------------------------------------------
423 : :
424 : : !!****f* m_splines/spline_bicubic
425 : : !! NAME
426 : : !! spline_bicubic
427 : : !!
428 : : !! FUNCTION
429 : : !! Generates coefficients for bicubic spline interpolation.
430 : : !!
431 : : !! INPUTS
432 : : !! n1 = length of first dimension
433 : : !! n2 = length of second dimension
434 : : !! x1 = positions on first dimension
435 : : !! x2 = positions on second dimension
436 : : !! y = function values on the (x1,x2) grid
437 : : !! der1_x1 = first derivative of y wrt x1
438 : : !! der1_x2 = first derivative of y wrt x2
439 : : !! der2_x1x2 = second-order cross-derivative of y wrt x1x2
440 : : !!
441 : : !! OUTPUT
442 : : !! spl_c = spline coefficients
443 : : !!
444 : : !! NOTES
445 : : !! Adapted from Numerical Recipes and libbci.
446 : : !!
447 : : !! SOURCE
448 : :
449 : 0 : subroutine spline_bicubic(n1,n2,x1,x2,y,der1_x1,der1_x2,der2_x1x2,spl_c)
450 : :
451 : : integer,intent(in) :: n1,n2
452 : : real(dp),intent(in) :: x1(n1),x2(n2),y(n1,n2)
453 : : real(dp),intent(in) :: der1_x1(n1,n2),der1_x2(n1,n2),der2_x1x2(n1,n2)
454 : : real(dp),intent(out):: spl_c(4,4,n1,n2)
455 : :
456 : : integer :: i1,i2
457 : : real(dp) :: dx1,dx2,wt(16,16),z(16)
458 : :
459 : : data wt /1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4, &
460 : : & 8*0,3,0,-9,6,-2,0,6,-4,10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6, &
461 : : & 2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3, &
462 : : & 2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,0,1,-2,1,5*0, &
463 : : & -3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,10*0,-3,3,2*0,2,-2,2*0, &
464 : : & -1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2, &
465 : : & -1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/
466 : :
467 : : ! Set coefficients for i1<n1 and i2<n2
468 [ # # ]: 0 : do i2 = 1,n2-1
469 [ # # ]: 0 : do i1 = 1,n1-1
470 : 0 : dx1 = x1(i1+1) - x1(i1)
471 : 0 : dx2 = x2(i2+1) - x2(i2)
472 : 0 : z(1) = y(i1,i2)
473 : 0 : z(2) = y(i1+1,i2)
474 : 0 : z(3) = y(i1+1,i2+1)
475 : 0 : z(4) = y(i1,i2+1)
476 : 0 : z(5) = der1_x1(i1,i2) * dx1
477 : 0 : z(6) = der1_x1(i1+1,i2) * dx1
478 : 0 : z(7) = der1_x1(i1+1,i2+1) * dx1
479 : 0 : z(8) = der1_x1(i1,i2+1) * dx1
480 : 0 : z(9) = der1_x2(i1,i2) * dx2
481 : 0 : z(10) = der1_x2(i1+1,i2) * dx2
482 : 0 : z(11) = der1_x2(i1+1,i2+1) * dx2
483 : 0 : z(12) = der1_x2(i1,i2+1) * dx2
484 : 0 : z(13) = der2_x1x2(i1,i2) * dx1 * dx2
485 : 0 : z(14) = der2_x1x2(i1+1,i2) * dx1 * dx2
486 : 0 : z(15) = der2_x1x2(i1+1,i2+1) * dx1 * dx2
487 : 0 : z(16) = der2_x1x2(i1,i2+1) * dx1 * dx2
488 [ # # # # ]: 0 : z = matmul(wt,z)
489 : 0 : spl_c(:,:,i1,i2) = reshape(z,(/4,4/),order=(/2,1/))
490 : : end do
491 : : end do
492 : :
493 : : ! Set coefficients for i1=n1 and i2=n2 (valid only at the border)
494 [ # # # # : 0 : spl_c(:,:,n1,:) = 0
# # ]
495 [ # # # # : 0 : spl_c(:,:,:,n2) = 0
# # ]
496 [ # # ]: 0 : spl_c(1,1,n1,:) = y(n1,:)
497 [ # # ]: 0 : spl_c(1,1,:,n2) = y(:,n2)
498 : :
499 : 0 : end subroutine spline_bicubic
500 : : !!***
501 : :
502 : : !----------------------------------------------------------------------
503 : :
504 : : !!****f* m_splines/spline_r
505 : : !! NAME
506 : : !! spline_r
507 : : !!
508 : : !! FUNCTION
509 : : !! Computes the spline of a real function.
510 : : !! If point lies outside the range of original grids, assign the extremal
511 : : !! point values to either head or tail.
512 : : !!
513 : : !! INPUTS
514 : : !! nomega_lo = number of point in the non regular grid (e.g. !logarithmic)
515 : : !! nomega_li = number of point in the regular grid on which the spline is computed
516 : : !! omega_lo = value of freq on the 1st grid
517 : : !! omega_li = value of freq on the 2nd grid
518 : : !! tospline_lo = function on the 1st grid
519 : : !!
520 : : !! OUTPUT
521 : : !! splined_lo = spline (on the 2nd grid)
522 : : !!
523 : : !! SOURCE
524 : :
525 : 0 : subroutine spline_r( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo, extrapolate)
526 : :
527 : : !Arguments --------------------------------------------
528 : : !scalars
529 : : integer, intent(in) :: nomega_lo, nomega_li
530 : : real(dp), intent(in) :: omega_lo(nomega_lo)
531 : : real(dp), intent(in) :: omega_li(nomega_li)
532 : : real(dp), intent(in) :: tospline_lo(nomega_lo)
533 : : real(dp), intent(out) :: splined_li(nomega_li)
534 : : logical, intent(in), optional :: extrapolate
535 : :
536 : : !Local variables---------------------------------------
537 : : !scalars
538 : : integer :: begin, end
539 : : real(dp) :: ybcbeg, ybcend
540 : 0 : real(dp) :: ysplin2_lo(nomega_lo)
541 : :
542 : 0 : ybcbeg=zero
543 : 0 : ybcend=zero
544 : :
545 : 0 : begin = 1
546 : 0 : end = nomega_li
547 : :
548 : 0 : call spline(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo)
549 [ # # ]: 0 : if (present(extrapolate)) then
550 [ # # ]: 0 : if (extrapolate) then
551 [ # # ]: 0 : do begin = 1, nomega_li
552 [ # # ]: 0 : if (omega_li(begin) >= omega_lo(1)) exit
553 : : end do
554 [ # # ]: 0 : do end = nomega_li, 1, -1
555 [ # # ]: 0 : if (omega_li(end) <= omega_lo(nomega_lo)) exit
556 : : end do
557 [ # # ]: 0 : ABI_CHECK(begin <= end, 'spline_c: omega_li not properly ordered')
558 : : end if
559 : : end if
560 [ # # ]: 0 : ABI_CHECK(begin <= end, 'spline_r: omega_li not properly ordered')
561 : 0 : call splint(nomega_lo, omega_lo, tospline_lo, ysplin2_lo, end-begin+1, omega_li(begin:end), splined_li(begin:end))
562 [ # # ]: 0 : if (present(extrapolate)) then
563 [ # # ]: 0 : if (extrapolate) then
564 [ # # ]: 0 : splined_li(1:begin-1) = tospline_lo(1)
565 [ # # ]: 0 : splined_li(end+1:nomega_li) = tospline_lo(nomega_lo)
566 : : end if
567 : : end if
568 : :
569 : 0 : end subroutine spline_r
570 : : !!***
571 : :
572 : : !----------------------------------------------------------------------
573 : :
574 : : !!****f* m_splines/spline_c
575 : : !! NAME
576 : : !! spline_c
577 : : !!
578 : : !! FUNCTION
579 : : !! Computes the spline of a complex function.
580 : : !!
581 : : !! INPUTS
582 : : !! nomega_lo = number of point in the non regular grid (e.g. !logarithmic)
583 : : !! nomega_li = number of point in the regular grid on which the spline is computed
584 : : !! omega_lo = value of freq on the 1st grid
585 : : !! omega_li = value of freq on the 2nd grid
586 : : !! tospline_lo = function on the 1st grid
587 : : !!
588 : : !! OUTPUT
589 : : !! splined_lo = spline (on the 2nd grid)
590 : : !!
591 : : !! SOURCE
592 : :
593 : 11570 : subroutine spline_c( nomega_lo, nomega_li, omega_lo, omega_li, splined_li, tospline_lo, extrapolate)
594 : :
595 : : !Arguments --------------------------------------------
596 : : !scalars
597 : : integer, intent(in) :: nomega_lo, nomega_li
598 : : real(dp), intent(in) :: omega_lo(nomega_lo)
599 : : real(dp), intent(in) :: omega_li(nomega_li)
600 : : complex(dp), intent(in) :: tospline_lo(nomega_lo)
601 : : complex(dp), intent(out) :: splined_li(nomega_li)
602 : : logical, intent(in), optional :: extrapolate
603 : :
604 : : !Local variables---------------------------------------
605 : : !scalars
606 : : integer :: begin, end
607 : : complex(dp) :: ybcbeg, ybcend
608 : 23140 : complex(dp) :: ysplin2_lo(nomega_lo)
609 : :
610 : 11570 : ybcbeg=czero
611 : 11570 : ybcend=czero
612 : :
613 : 11570 : begin = 1
614 : 11570 : end = nomega_li
615 : :
616 : 11570 : call spline_complex(omega_lo, tospline_lo, nomega_lo, ybcbeg, ybcend, ysplin2_lo)
617 [ + + ]: 11570 : if (present(extrapolate)) then
618 [ + - ]: 11520 : if (extrapolate) then
619 [ - + ]: 23040 : do begin = 1, nomega_li
620 [ + + ]: 23040 : if (omega_li(begin) >= omega_lo(1)) exit
621 : : end do
622 [ + - ]: 23040 : do end = nomega_li, 1, -1
623 [ + + ]: 23040 : if (omega_li(end) <= omega_lo(nomega_lo)) exit
624 : : end do
625 [ - + ]: 11520 : ABI_CHECK(begin <= end, 'spline_c: omega_li not properly ordered')
626 : : end if
627 : : end if
628 : 11570 : call splint_complex( nomega_lo, omega_lo, tospline_lo,ysplin2_lo, end-begin+1, omega_li(begin:end), splined_li(begin:end))
629 [ + + ]: 11570 : if (present(extrapolate)) then
630 [ + - ]: 11520 : if (extrapolate) then
631 [ + + ]: 23040 : splined_li(1:begin-1) = tospline_lo(1)
632 [ + + ]: 23040 : splined_li(end+1:nomega_li) = tospline_lo(nomega_lo)
633 : : end if
634 : : end if
635 : :
636 : 11570 : end subroutine spline_c
637 : : !!***
638 : :
639 : : !----------------------------------------------------------------------
640 : :
641 : : !!****f* m_splines/spline_complex
642 : : !! NAME
643 : : !! spline_complex
644 : : !!
645 : : !! FUNCTION
646 : : !! spline_complex interfaces the usual spline routine in a case of a
647 : : !! complex function
648 : : !!
649 : : !! INPUTS
650 : : !! Input, integer N, the number of data points; N must be at least 2.
651 : : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
652 : : !! spline will actually be linear.
653 : : !!
654 : : !! Input, double precision T(N), the knot values, that is, the points where data
655 : : !! is specified. The knot values should be distinct, and increasing.
656 : : !!
657 : : !! Input, complex Y(N), the data values to be interpolated.
658 : : !!
659 : : !! Input, complex YBCBEG, YBCEND, the values to be used in the boundary
660 : : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
661 : : !!
662 : : !! OUTPUT
663 : : !! Output, complex YPP(N), the second derivatives of the cubic spline.
664 : : !!
665 : : !! SOURCE
666 : :
667 : 29994 : subroutine spline_complex( t, y, n, ybcbeg, ybcend, ypp )
668 : :
669 : : integer, intent(in) :: n
670 : : real(dp), intent(in) :: t(n)
671 : : complex(dp), intent(in) :: y(n)
672 : : complex(dp), intent(in) :: ybcbeg
673 : : complex(dp), intent(in) :: ybcend
674 : : complex(dp), intent(out) :: ypp(n)
675 : :
676 : 29994 : real(dp), allocatable :: y_r(:)
677 : : real(dp) :: ybcbeg_r
678 : : real(dp) :: ybcend_r
679 : 29994 : real(dp), allocatable :: ypp_r(:)
680 : 29994 : real(dp), allocatable :: y_i(:)
681 : : real(dp) :: ybcbeg_i
682 : : real(dp) :: ybcend_i
683 : 29994 : real(dp), allocatable :: ypp_i(:)
684 : :
685 [ - + + - : 89982 : ABI_MALLOC(y_r,(n))
+ - + - -
+ - + ]
686 [ - + + - : 59988 : ABI_MALLOC(ypp_r,(n))
- + - + ]
687 [ - + + - : 59988 : ABI_MALLOC(y_i,(n))
- + - + ]
688 [ - + + - : 59988 : ABI_MALLOC(ypp_i,(n))
- + - + ]
689 [ - + + + ]: 8145078 : y_r=real(y)
690 [ - + + + ]: 8145078 : y_i=aimag(y) !vz_d
691 : 29994 : ybcbeg_r=real(ybcbeg)
692 : 29994 : ybcbeg_i=aimag(ybcbeg) !vz_d
693 : 29994 : ybcend_r=real(ybcend)
694 : 29994 : ybcend_i=aimag(ybcend) !vz_d
695 : 29994 : call spline( t, y_r, n, ybcbeg_r, ybcend_r, ypp_r )
696 : 29994 : call spline( t, y_i, n, ybcbeg_i, ybcend_i, ypp_i )
697 [ + + ]: 8115084 : ypp=cmplx(ypp_r,ypp_i)
698 [ - + ]: 29994 : ABI_FREE(y_r)
699 [ - + ]: 29994 : ABI_FREE(ypp_r)
700 [ - + ]: 29994 : ABI_FREE(y_i)
701 [ - + ]: 29994 : ABI_FREE(ypp_i)
702 : :
703 : 29994 : end subroutine spline_complex
704 : : !!***
705 : :
706 : : !----------------------------------------------------------------------
707 : :
708 : : !!****f* m_splines/splint
709 : : !! NAME
710 : : !! splint
711 : : !!
712 : : !! FUNCTION
713 : : !! Compute spline interpolation. There is no hypothesis
714 : : !! about the spacing of the input grid points.
715 : : !!
716 : : !! INPUTS
717 : : !! nspline: number of grid points of input mesh
718 : : !! xspline(nspline): input mesh
719 : : !! yspline(nspline): function on input mesh
720 : : !! ysplin2(nspline): second derivative of yspline on input mesh
721 : : !! nfit: number of points of output mesh
722 : : !! xfit(nfit): output mesh
723 : : !!
724 : : !! OUTPUT
725 : : !! yfit(nfit): function on output mesh
726 : : !! [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline).
727 : : !! The input value is incremented by the number of such points.
728 : : !!
729 : : !! SOURCE
730 : :
731 : 50913889 : subroutine splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr)
732 : :
733 : : integer, intent(in) :: nfit, nspline
734 : : integer,optional,intent(out) :: ierr
735 : : real(dp), intent(in) :: xspline(nspline)
736 : : real(dp), intent(in) :: yspline(nspline)
737 : : real(dp), intent(in) :: ysplin2(nspline)
738 : : real(dp), intent(in) :: xfit(nfit)
739 : : real(dp), intent(out) :: yfit(nfit)
740 : :
741 : : !local
742 : : integer :: left,i,k,right,my_err
743 : : real(dp) :: delarg,invdelarg,aa,bb
744 : :
745 : : !source
746 : :
747 : 50913889 : my_err=0
748 : :
749 : 50913889 : left = 1
750 [ + + ]: 12823368573 : do i=1, nfit
751 : 12772454684 : yfit(i)=0.d0 ! Initialize for the unlikely event that rmax exceed r(mesh)
752 : : !
753 [ + + ]: 13051317340 : do k=left+1, nspline
754 [ + + ]: 13051317340 : if(xspline(k) >= xfit(i)) then
755 [ - + ]: 12757177982 : if(xspline(k-1) <= xfit(i)) then
756 : : right = k
757 : : left = k-1
758 : : else
759 [ # # ]: 0 : if (k-1.eq.1 .and. i.eq.1) then
760 : 0 : ABI_ERROR('xfit(1) < xspline(1)')
761 : : !my_err=my_err+1
762 : : !exit
763 : : else
764 : 0 : ABI_ERROR('xfit not properly ordered')
765 : : end if
766 : : end if
767 : 12757177982 : delarg= xspline(right) - xspline(left)
768 : 12757177982 : invdelarg= 1.0d0/delarg
769 : 12757177982 : aa= (xspline(right)-xfit(i))*invdelarg
770 : 12757177982 : bb= (xfit(i)-xspline(left))*invdelarg
771 : :
772 : : yfit(i) = aa*yspline(left) + bb*yspline(right) &
773 : : & +( (aa*aa*aa-aa)*ysplin2(left) + &
774 : 12757177982 : & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0d0
775 : 12757177982 : exit
776 : : end if
777 : : end do ! k
778 : : !
779 [ + + ]: 12823368573 : if (k==nspline+1) my_err=my_err+1 ! xfit not found
780 : : end do ! i
781 : :
782 [ + + ]: 50913889 : if (PRESENT(ierr)) ierr=my_err
783 : :
784 : 50913889 : end subroutine splint
785 : : !!***
786 : :
787 : : !----------------------------------------------------------------------
788 : :
789 : : !!****f* m_splines/splint_complex
790 : : !! NAME
791 : : !! splint_complex
792 : : !!
793 : : !! FUNCTION
794 : : !! Interface to the usual splint to compute *complex* spline interpolation. There is no hypothesis
795 : : !! about the spacing of the input grid points.
796 : : !!
797 : : !! INPUTS
798 : : !! nspline: number of grid points of input mesh
799 : : !! xspline(nspline): input mesh
800 : : !! yspline(nspline): complex function on input mesh
801 : : !! ysplin2(nspline): second derivative of yspline on input mesh
802 : : !! nfit: number of points of output mesh
803 : : !! xfit(nfit): output mesh
804 : : !!
805 : : !! OUTPUT
806 : : !! yfit(nfit): complex function on output mesh
807 : : !!
808 : : !! SOURCE
809 : :
810 : 29994 : subroutine splint_complex (nspline,xspline,yspline,ysplin2,nfit,xfit,yfit)
811 : :
812 : : integer, intent(in) :: nfit, nspline
813 : : real(dp), intent(in) :: xspline(nspline)
814 : : complex(dp), intent(in) :: yspline(nspline)
815 : : complex(dp), intent(in) :: ysplin2(nspline)
816 : : real(dp), intent(in) :: xfit(nfit)
817 : : complex(dp), intent(out) :: yfit(nfit)
818 : :
819 : 29994 : real(dp), allocatable :: ysplin2_r(:)
820 : 29994 : real(dp), allocatable :: ysplin2_i(:)
821 : 29994 : real(dp), allocatable :: yspline_r(:)
822 : 29994 : real(dp), allocatable :: yspline_i(:)
823 : 29994 : real(dp), allocatable :: yfit_r(:)
824 : 29994 : real(dp), allocatable :: yfit_i(:)
825 : :
826 [ - + + - : 89982 : ABI_MALLOC(yspline_r,(nspline))
+ - + - -
+ - + ]
827 [ - + + - : 59988 : ABI_MALLOC(yspline_i,(nspline))
- + - + ]
828 [ - + + - : 59988 : ABI_MALLOC(ysplin2_r,(nspline))
- + - + ]
829 [ - + + - : 59988 : ABI_MALLOC(ysplin2_i,(nspline))
- + - + ]
830 [ - + + - : 89982 : ABI_MALLOC(yfit_r,(nfit))
+ - + - -
+ - + ]
831 [ - + + - : 59988 : ABI_MALLOC(yfit_i,(nfit))
- + - + ]
832 : :
833 : : !local
834 : :
835 : : !source
836 [ - + + + ]: 8145078 : yspline_r=real(yspline)
837 [ - + + + ]: 8145078 : yspline_i=aimag(yspline) !vz_d
838 [ - + + + ]: 8145078 : ysplin2_r=real(ysplin2)
839 [ - + + + ]: 8145078 : ysplin2_i=aimag(ysplin2) !vz_d
840 : 29994 : call splint (nspline,xspline,yspline_r,ysplin2_r,nfit,xfit,yfit_r)
841 : 29994 : call splint (nspline,xspline,yspline_i,ysplin2_i,nfit,xfit,yfit_i)
842 [ + + ]: 6260473804 : yfit=cmplx(yfit_r,yfit_i)
843 [ - + ]: 29994 : ABI_FREE(yspline_r)
844 [ - + ]: 29994 : ABI_FREE(yspline_i)
845 [ - + ]: 29994 : ABI_FREE(ysplin2_r)
846 [ - + ]: 29994 : ABI_FREE(ysplin2_i)
847 [ - + ]: 29994 : ABI_FREE(yfit_r)
848 [ - + ]: 29994 : ABI_FREE(yfit_i)
849 : :
850 : 29994 : end subroutine splint_complex
851 : : !!***
852 : :
853 : : !!****f* m_splines/spline_integrate
854 : : !! NAME
855 : : !! spline_integrate
856 : : !!
857 : : !! FUNCTION
858 : : !! Calculates an integral using cubic spline interpolation.
859 : : !!
860 : : !! INPUTS
861 : : !! npts= number of grid points of input mesh
862 : : !! dx= step of input mesh
863 : : !! integrand= function on input mesh
864 : : !!
865 : : !! OUTPUT
866 : : !! integral= integral of the input function
867 : : !!
868 : : !! SOURCE
869 : :
870 : 0 : subroutine spline_integrate(integral,npts,dx,integrand)
871 : :
872 : : integer,intent(in) :: npts
873 : : real(dp),intent(out) :: integral
874 : : real(dp),intent(in) :: dx,integrand(npts)
875 : :
876 : : integer :: ix
877 : 0 : real(dp) :: ptmp,sf(npts),sf_der2(npts),sf_mesh(npts),utmp(npts)
878 : :
879 : : ! Prepare mesh
880 [ # # ]: 0 : forall (ix=1:npts) sf_mesh(ix) = (ix - 1) * dx
881 : :
882 : : ! Calculate second derivative of integrand (adapted from Numercial Recipes)
883 : 0 : sf_der2(1) = zero
884 : 0 : sf_der2(npts) = zero
885 : 0 : utmp(1) = zero
886 : :
887 [ # # ]: 0 : do ix=2,npts-1
888 : 0 : ptmp = half * sf_der2(ix-1) + two
889 : 0 : sf_der2(ix) = (half - one) / ptmp
890 : : utmp(ix) = (three * (integrand(ix+1) + integrand(ix-1) - &
891 : 0 : & two*integrand(ix)) / (dx**2) - half * utmp(ix-1)) / ptmp
892 : : end do
893 [ # # ]: 0 : do ix=npts-1,1,-1
894 : 0 : sf_der2(ix) = sf_der2(ix) * sf_der2(ix+1) + utmp(ix)
895 : : end do
896 : :
897 : : ! Actually calculate integral
898 [ # # ]: 0 : sf(:) = integrand(:) * dx
899 : : integral = (sf(1) + sf(npts)) / 2.0_dp - &
900 : : & (sf_der2(1) + sf_der2(npts)) / 24.0_dp + &
901 [ # # # # ]: 0 : & sum(sf(2:npts-1)) - sum(sf_der2(2:npts-1)) / 12.0_dp
902 : :
903 : 0 : end subroutine spline_integrate
904 : : !!***
905 : :
906 : : !!****f* m_splines/intrpl
907 : : !! NAME
908 : : !! intrpl
909 : : !!
910 : : !! FUNCTION
911 : : !!
912 : : !! DOUBLE PRECISION INTERPOLATION OF A SINGLE VALUED FUNCTION.
913 : : !! THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
914 : : !! GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE
915 : : !! AND FOR A GIVEN SET OF X VALUES(ABSCISSAE),THE VALUES OF
916 : : !! A SINGLE VALUED FUNCTION Y=Y(X).
917 : : !!
918 : : !! THE SUBROUTINE ALSO CALCULATES FIRST DERIVATIVES DV(X) AND
919 : : !! SECOND DERIVATIVE DV2(X)
920 : : !
921 : : !! THE INPUT PARAMETERS ARE;
922 : : !!
923 : : !! L=NUMBER OF DATA POINTS
924 : : !! (MUST BE TWO OR GREATER)
925 : : !! X=ARRAY OF DIMENSION L STORING THE X VALUES
926 : : !! OF INPUT DATA POINTS (IN ASCENDING ORDER)
927 : : !! Y=ARRAY OF DIMENSION L STORING THE Y VALUES OF INPUT DATA POINTS
928 : : !! N=NUMBER OF POINTS AT WHICH INTERPOLATION OF THE Y-VALUES
929 : : !! IS REQUIRED (MUST BE 1 OR GREATER)
930 : : !! U=ARRAY OF DIMENSION N STORING THE X VALUES
931 : : !! OF THE DESIRED POINTS
932 : : !!
933 : : !! THE OUTPUT PARAMETER IS V=ARRAY OF DIMENSION N WHERE THE
934 : : !! INTERPOLATED Y VALUES ARE TO BE DISPLAYED
935 : : !!
936 : : !! INPUTS
937 : : !! CUBIC SPLINE INTERPOLATION
938 : : !!
939 : : !! OUTPUT
940 : : !!
941 : : !! NOTES
942 : : !! This routine is deprecated and will be replaced by the other routines of this module.
943 : : !!
944 : : !! SOURCE
945 : :
946 : 0 : SUBROUTINE INTRPL(L,X,Y,N,U,V,dv,dv2,ideriv)
947 : :
948 : : IMPLICIT DOUBLE PRECISION (A-H,O-Z)
949 : : IMPLICIT INTEGER(I-N)
950 : : !
951 : : PARAMETER (NQQ=12000)
952 : :
953 : : COMMON/QQ/ QQ(4,NQQ)
954 : : DIMENSION X(L),Y(L),U(N),V(N),DV(NQQ),DV2(NQQ)
955 : : EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3)
956 : : REAL*8 M1,M2,M3,M4,M5
957 : : EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5),&
958 : : & (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
959 : : !
960 : : ! PRELIMINARY PROCESSING
961 : :
962 : 0 : L0=L
963 : 0 : LM1=L0-1
964 : 0 : LM2=LM1-1
965 : 0 : LP1=L0+1
966 : 0 : N0=N
967 [ # # ]: 0 : IF(N0.GT.NQQ) THEN
968 : 0 : NQQV=NQQ
969 : 0 : write(std_out,2089) NQQV,N0
970 : : ! CALL EXIT
971 : : END IF
972 [ # # ]: 0 : IF(LM2.LT.0) GO TO 90
973 [ # # ]: 0 : IF(N0.LE.0) GO TO 91
974 [ # # ]: 0 : DO 11 I=2,L0
975 : :
976 : : ! IF(X(I-1)-X(I))11,95,96
977 [ # # ]: 0 : IF(X(I-1)-X(I).EQ.0.0D0) GO TO 95
978 [ # # ]: 0 : IF(X(I-1)-X(I).GT.0.0D0) GO TO 96
979 : 0 : 11 CONTINUE
980 : : IPV=0
981 : : !
982 : : !*** MAIN LOOP
983 : : FINT=0.0D0
984 [ # # ]: 0 : DO 80 K=1,N0
985 : 0 : UK=U(K)
986 : : !
987 : : !*** ROUTINE TO LOCATE THE DESIRED POINT
988 [ # # ]: 0 : IF(UK.GE.X(L0)) GO TO 26
989 [ # # ]: 0 : IF(UK.LT.X(1)) GO TO 25
990 : 0 : IMN=2
991 : 0 : IMX=L0
992 : 0 : 21 I=(IMN+IMX)/2
993 [ # # ]: 0 : IF(UK.GE.X(I)) GO TO 23
994 : 0 : IMX=I
995 : 0 : GO TO 24
996 : 0 : 23 IMN=I+1
997 [ # # ]: 0 : 24 IF(IMX.GT.IMN) GO TO 21
998 : 0 : I=IMX
999 : 0 : GO TO 30
1000 : 0 : 25 I=1
1001 : 0 : GO TO 30
1002 : 0 : 26 I=LP1
1003 : 0 : GO TO 30
1004 : : !
1005 : : !*** CHECK IF I=IPV
1006 [ # # ]: 0 : 30 IF(I.EQ.IPV) GO TO 70
1007 : 0 : IPV=I
1008 : : !
1009 : : !*** ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND TO
1010 : : !*** ESTIMATE THEM IF NECESSARY
1011 : 0 : J=I
1012 [ # # ]: 0 : IF(J.EQ.1) J=2
1013 [ # # ]: 0 : IF(J.EQ.LP1) J=L0
1014 : 0 : X3=X(J-1)
1015 : 0 : Y3=Y(J-1)
1016 : 0 : X4=X(J)
1017 : 0 : Y4=Y(J)
1018 : 0 : A3=X4-X3
1019 : 0 : M3=(Y4-Y3)/A3
1020 [ # # ]: 0 : IF(LM2.EQ.0) GO TO 43
1021 [ # # ]: 0 : IF(J.EQ.2) GO TO 41
1022 : 0 : X2=X(J-2)
1023 : 0 : Y2=Y(J-2)
1024 : 0 : A2=X3-X2
1025 : 0 : M2=(Y3-Y2)/A2
1026 [ # # ]: 0 : IF(J.EQ.L0) GO TO 42
1027 : 0 : 41 X5=X(J+1)
1028 : 0 : Y5=Y(J+1)
1029 : 0 : A4=X5-X4
1030 : 0 : M4=(Y5-Y4)/A4
1031 [ # # ]: 0 : IF(J.EQ.2) M2=M3+M3-M4
1032 : 0 : GO TO 45
1033 : 0 : 42 M4=M3+M3-M2
1034 : 0 : GO TO 45
1035 : : 43 M2=M3
1036 [ # # ]: 0 : 45 IF(J.LE.3) GO TO 46
1037 : 0 : A1=X2-X(J-3)
1038 : 0 : M1=(Y2-Y(J-3))/A1
1039 : 0 : GO TO 47
1040 : 0 : 46 M1=M2+M2-M3
1041 [ # # ]: 0 : 47 IF(J.GE.LM1) GO TO 48
1042 : 0 : A5=X(J+2)-X5
1043 : 0 : M5=(Y(J+2)-Y5)/A5
1044 : 0 : GO TO 50
1045 : 0 : 48 M5=M4+M4-M3
1046 : : !
1047 : : !*** NUMERICAL DIFFERENTIATION
1048 [ # # ]: 0 : 50 IF(I.EQ.LP1) GO TO 52
1049 : 0 : W2=ABS(M4-M3)
1050 : 0 : W3=ABS(M2-M1)
1051 : 0 : SW=W2+W3
1052 [ # # ]: 0 : IF(SW.NE.0.0) GO TO 51
1053 : : W2=0.5D0
1054 : : W3=0.5D0
1055 : 0 : SW=1.0D0
1056 : 0 : 51 T3=(W2*M2+W3*M3)/SW
1057 [ # # ]: 0 : IF(I.EQ.1) GO TO 54
1058 : 0 : 52 W3=ABS(M5-M4)
1059 : 0 : W4=ABS(M3-M2)
1060 : 0 : SW=W3+W4
1061 [ # # ]: 0 : IF(SW.NE.0.0) GO TO 53
1062 : : W3=0.5D0
1063 : : W4=0.5D0
1064 : 0 : SW=1.0D0
1065 : 0 : 53 T4=(W3*M3+W4*M4)/SW
1066 [ # # ]: 0 : IF(I.NE.LP1) GO TO 60
1067 : : T3=T4
1068 : 0 : SA=A2+A3
1069 : 0 : T4=0.5D0*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA))
1070 : : X3=X4
1071 : : Y3=Y4
1072 : 0 : A3=A2
1073 : 0 : M3=M4
1074 : 0 : GO TO 60
1075 : 0 : 54 T4=T3
1076 : 0 : SA=A3+A4
1077 : 0 : T3=0.5D0*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA))
1078 : 0 : X3=X3-A4
1079 : 0 : Y3=Y3-M2*A4
1080 : 0 : A3=A4
1081 : 0 : M3=M2
1082 : : !
1083 : : !*** COMPUTATION OF THE POLYNOMIAL
1084 : 0 : 60 Q2=(2.0D0*(M3-T3)+M3-T4)/A3
1085 : 0 : Q3=(-M3-M3+T3+T4)/(A3*A3)
1086 : 0 : 70 DX=UK-P0
1087 : 0 : V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3))
1088 : :
1089 [ # # ]: 0 : IF(IDERIV.EQ.0) GO TO 80
1090 : 0 : DV(K)=Q1+DX*(2.0D0*Q2+DX*3.0D0*Q3)
1091 : 0 : DV2(k)=6.0D0*Q3*DX+2.d0*Q2
1092 : 0 : QQ(1,K)=Q0
1093 : 0 : QQ(2,K)=Q1
1094 : 0 : QQ(3,K)=Q2
1095 : 0 : QQ(4,K)=Q3
1096 : 0 : 80 CONTINUE
1097 : 0 : RETURN
1098 : : !
1099 : : !*** ERROR EXIT
1100 : 0 : 90 write(std_out,2090)
1101 : 0 : GO TO 99
1102 : 0 : 91 write(std_out,2091)
1103 : 0 : GO TO 99
1104 : 0 : 95 write(std_out,2095)
1105 : 0 : GO TO 97
1106 : 0 : 96 write(std_out,2096)
1107 : 0 : 97 write(std_out,2097)I,X(I)
1108 : 0 : 99 write(std_out,2099) L0,N0
1109 : 0 : RETURN
1110 : : !
1111 : : !*** FORMAT STATEMENTS
1112 : : 2089 FORMAT( 'WARNING ERROR IN INTRPL. MAX ALLOWED VALUE OF N0 IS',&
1113 : : & I3,' HERE N0 IS',I3)
1114 : : 2090 FORMAT(1X/' N = 1 OR LESS.'/)
1115 : : 2091 FORMAT(1X/' N = 0 OR LESS.'/)
1116 : : 2095 FORMAT(1X/' IDENTICAL X VALUES.'/)
1117 : : 2096 FORMAT(1X/' X VALUES OUT OF SEQUENCE.'/)
1118 : : 2097 FORMAT(4X,'I =',I7,10X,6X,'X(I) =',E12.3)
1119 : : 2099 FORMAT(4X,'L =',I7,10X,3X,'N =',I7/ &
1120 : : & ' ERROR DETECTED IN ROUTINE INTRPL')
1121 : : !
1122 : : END subroutine intrpl
1123 : : !!***
1124 : :
1125 : : !!****f* m_splines/spline2
1126 : : !! NAME
1127 : : !! spline2
1128 : : !!
1129 : : !! FUNCTION
1130 : : !! SPLINE2 computes the first derivatives of a cubic spline.
1131 : : !!
1132 : : !! INPUTS
1133 : : !! Input, integer N, the number of data points; N must be at least 2.
1134 : : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
1135 : : !! spline will actually be linear.
1136 : : !!
1137 : : !! Input, double precision X(N), the knot values, that is, the points where data
1138 : : !! is specified. The knot values should be distinct, and increasing.
1139 : : !!
1140 : : !! Input, double precision Y(N), the data values to be interpolated.
1141 : : !!
1142 : : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
1143 : : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
1144 : : !!
1145 : : !! Input, integer IBCBEG, IBCEND, the type of boundary conditions at the first / last point
1146 : : !! If 1, the value of the first derivative is constrained to YBCBEG / YBCEND
1147 : : !! If 2, the value of the second derivative is constrained to YBCBEG / YBCEND
1148 : : !! If 3, the third derivative is continuous at the second / second to last point (not-a-knot)
1149 : : !!
1150 : : !! OUTPUT
1151 : : !! Output, double precision YP(N), the first derivatives of the cubic spline.
1152 : : !! On [x(i),x(i+1)], the spline y(i)+yp(i)*(x-x(i))+c*(x-x(i))**2+d*(x-x(i))**3 can then
1153 : : !! be reconstructed with c = (3*s-2*yp(i)-yp(i+1))/(x(i+1)-x(i))
1154 : : !! d = (yp(i)+yp(i+1)-2*s)/(x(i+1)-x(i))**2
1155 : : !! s = (y(i+1)-y(i))/(x(i+1)-x(i))
1156 : : !!
1157 : : !! Solving the tridiagonal system on the first derivatives rather than the
1158 : : !! second derivatives as in the regular spline subroutine allows for the
1159 : : !! fantastic "not-a-knot" boundary condition.
1160 : : !!
1161 : : !! SOURCE
1162 : :
1163 : 0 : subroutine spline2(x,y,n,yp,ybcbeg,ybcend,ibcbeg,ibcend)
1164 : :
1165 : : integer,intent(in) :: n,ibcbeg,ibcend
1166 : : real(dp),intent(in) :: ybcbeg,ybcend
1167 : : real(dp),intent(in) :: x(n),y(n)
1168 : : real(dp),intent(inout) :: yp(n)
1169 : :
1170 : : integer :: i
1171 : : real(dp) :: pinv,ratio
1172 : 0 : real(dp),allocatable :: tmp(:)
1173 : :
1174 [ # # # # : 0 : ABI_MALLOC(tmp,(n))
# # # # #
# # # ]
1175 : :
1176 [ # # ]: 0 : if (ibcbeg==1) then
1177 : 0 : yp(1)=zero ; tmp(1)=ybcbeg
1178 [ # # ]: 0 : else if (ibcbeg==2) then
1179 : 0 : yp(1)=-half ; tmp(1)=(six*((y(2)-y(1))/(x(2)-x(1)))-ybcbeg*(x(2)-x(1)))/four
1180 [ # # ]: 0 : else if (ibcbeg==3) then
1181 : 0 : yp(1)=(x(1)-x(3))/(x(3)-x(2))
1182 : : tmp(1)=((y(2)-y(1))*(x(3)-x(2))*(two*x(3)+x(2)-three*x(1))/(x(2)-x(1))+(y(3)-y(2))* &
1183 : 0 : & (x(2)-x(1))**2/(x(3)-x(2)))/((x(3)-x(2))*(x(3)-x(1)))
1184 : : else
1185 : 0 : ABI_BUG("Option not recognized for ibcbeg !")
1186 : : end if
1187 : :
1188 [ # # ]: 0 : if (ibcend==1) then
1189 : 0 : yp(n)=zero ; tmp(n)=ybcend
1190 [ # # ]: 0 : else if (ibcend==2) then
1191 : 0 : yp(n)=half ; tmp(n)=(six*(y(n)-y(n-1))/(x(n)-x(n-1))+ybcend*(x(n)-x(n-1)))/four
1192 [ # # ]: 0 : else if (ibcend==3) then
1193 : 0 : yp(n)=(x(n)-x(n-2))/(x(n-1)-x(n-2))
1194 : : tmp(n)=((y(n-1)-y(n-2))*(x(n)-x(n-1))**2/(x(n-1)-x(n-2))+ &
1195 : : & (y(n)-y(n-1))*(x(n-1)-x(n-2))*(three*x(n)-x(n-1)-two*x(n-2)) &
1196 : 0 : & /(x(n)-x(n-1)))/((x(n-1)-x(n-2))*(x(n)-x(n-2)))
1197 : : else
1198 : 0 : ABI_BUG("Option not recognized for ibcend !")
1199 : : end if
1200 : :
1201 [ # # ]: 0 : do i=2,n-1
1202 : 0 : ratio = (x(i+1)-x(i))/(x(i+1)-x(i-1))
1203 : 0 : pinv = one/(two+ratio*yp(i-1))
1204 : 0 : yp(i) = (ratio-one)*pinv
1205 : : tmp(i) = (three*((y(i)-y(i-1))*(x(i+1)-x(i))/(x(i)-x(i-1))+ &
1206 : : & (y(i+1)-y(i))*(x(i)-x(i-1))/(x(i+1)-x(i)))/(x(i+1)-x(i-1))-&
1207 : 0 : & ratio*tmp(i-1))*pinv
1208 [ # # ]: 0 : if (abs(tmp(i))<1.d5*tiny(zero)) tmp(i) = zero !MT20050927
1209 : : end do
1210 : :
1211 : 0 : yp(n) = (tmp(n)-yp(n)*tmp(n-1))/(yp(n)*yp(n-1)+one)
1212 [ # # ]: 0 : do i=n-1,1,-1
1213 : 0 : yp(i)=yp(i)*yp(i+1)+tmp(i)
1214 : : end do
1215 : :
1216 : 0 : ABI_FREE(tmp)
1217 : :
1218 : 0 : end subroutine spline2
1219 : : !!***
1220 : :
1221 : : !!****f* m_splines/spline2_complex
1222 : : !! NAME
1223 : : !! spline2_complex
1224 : : !!
1225 : : !! FUNCTION
1226 : : !! SPLINE2_COMPLEX computes the first derivatives of a cubic spline for a complex function.
1227 : : !!
1228 : : !! INPUTS
1229 : : !! Input, integer N, the number of data points; N must be at least 2.
1230 : : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
1231 : : !! spline will actually be linear.
1232 : : !!
1233 : : !! Input, double precision X(N), the knot values, that is, the points where data
1234 : : !! is specified. The knot values should be distinct, and increasing.
1235 : : !!
1236 : : !! Input, double precision Y(N), the data values to be interpolated.
1237 : : !!
1238 : : !! Input, double precision YBCBEG, YBCEND, the values to be used in the boundary
1239 : : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
1240 : : !!
1241 : : !! Input, integer IBCBEG, IBCEND, the type of boundary conditions at the first / last point
1242 : : !! If 1, the value of the first derivative is constrained to YBCBEG / YBCEND
1243 : : !! If 2, the value of the second derivative is constrained to YBCBEG / YBCEND
1244 : : !! If 3, the third derivative is continuous at the second / second to last point (not-a-knot)
1245 : : !!
1246 : : !! OUTPUT
1247 : : !! Output, double precision YP(N), the first derivatives of the cubic spline.
1248 : : !! On [x(i),x(i+1)], the spline y(i)+yp(i)*(x-x(i))+c*(x-x(i))**2+d*(x-x(i))**3 can then
1249 : : !! be reconstructed with c = (3*s-2*yp(i)-yp(i+1))/(x(i+1)-x(i))
1250 : : !! d = (yp(i)+yp(i+1)-2*s)/(x(i+1)-x(i))**2
1251 : : !! s = (y(i+1)-y(i))/(x(i+1)-x(i))
1252 : : !!
1253 : : !! Solving the tridiagonal system on the first derivatives rather than the
1254 : : !! second derivatives as in the regular spline subroutine allows for the
1255 : : !! fantastic "not-a-knot" boundary condition.
1256 : : !!
1257 : : !! SOURCE
1258 : :
1259 : 0 : subroutine spline2_complex(x,y,n,yp,ybcbeg,ybcend,ibcbeg,ibcend)
1260 : :
1261 : : integer,intent(in) :: n,ibcbeg,ibcend
1262 : : complex(dp),intent(in) :: ybcbeg,ybcend
1263 : : real(dp),intent(in) :: x(n)
1264 : : complex(dp),intent(in) :: y(n)
1265 : : complex(dp),intent(inout) :: yp(n)
1266 : :
1267 : 0 : real(dp),allocatable :: y_tmp(:),yp_i(:),yp_r(:)
1268 : :
1269 [ # # # # : 0 : ABI_MALLOC(y_tmp,(n))
# # # # #
# # # ]
1270 [ # # # # : 0 : ABI_MALLOC(yp_i,(n))
# # # # ]
1271 [ # # # # : 0 : ABI_MALLOC(yp_r,(n))
# # # # ]
1272 : :
1273 [ # # ]: 0 : y_tmp(:) = dble(y(:))
1274 : 0 : call spline2(x(:),y_tmp(:),n,yp_r(:),dble(ybcbeg),dble(ybcend),ibcbeg,ibcend)
1275 [ # # ]: 0 : y_tmp(:) = aimag(y(:))
1276 : 0 : call spline2(x(:),y_tmp(:),n,yp_i(:),aimag(ybcbeg),aimag(ybcend),ibcbeg,ibcend)
1277 [ # # ]: 0 : yp(:) = cmplx(yp_r(:),yp_i(:),kind=dp)
1278 : :
1279 [ # # ]: 0 : ABI_FREE(y_tmp)
1280 [ # # ]: 0 : ABI_FREE(yp_i)
1281 [ # # ]: 0 : ABI_FREE(yp_r)
1282 : :
1283 : 0 : end subroutine spline2_complex
1284 : : !!***
1285 : :
1286 : : end module m_splines
1287 : : !!***
|
