1
2! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
3! This file is distributed under the terms of the GNU Lesser General Public
4! License. See the file COPYING for license details.
5
6!BOP
7! !ROUTINE: rfinterp
8! !INTERFACE:
9subroutine rfinterp(ni,xi,wci,fi,no,xo,fo)
10! !INPUT/OUTPUT PARAMETERS:
11!   ni  : number of input points (in,integer)
12!   xi  : input abscissa array (in,real(ni))
13!   wci : input spline coefficient weights (in,real(12,ni))
14!   fi  : input data array (in,real(ni))
15!   no  : number of output points (in,integer)
16!   xo  : output abscissa array (in,real(no))
17!   fo  : output interpolated function (out,real(no))
18! !DESCRIPTION:
19!   Given a function defined on a set of input points, this routine uses a
20!   clamped cubic spline to interpolate the function on a different set of
21!   points. See routine {\tt spline}.
22!
23! !REVISION HISTORY:
24!   Created January 2005 (JKD)
25!   Arguments changed, April 2016 (JKD)
26!EOP
27!BOC
28implicit none
29! arguments
30integer, intent(in) :: ni
31real(8), intent(in) :: xi(ni),wci(12,ni),fi(ni)
32integer, intent(in) :: no
33real(8), intent(in) :: xo(no)
34real(8), intent(out) :: fo(no)
35! local variables
36integer i,j,k,l
37real(8) x,dx
38! automatic arrays
39real(8) cf(3,ni)
40if (ni.le.0) then
41  write(*,*)
42  write(*,'("Error(rfinterp): invalid number of input points : ",I8)') ni
43  write(*,*)
44  stop
45end if
46if (no.le.0) then
47  write(*,*)
48  write(*,'("Error(rfinterp): invalid number of output points : ",I8)') no
49  write(*,*)
50  stop
51end if
52if (ni.eq.1) then
53  fo(:)=fi(1)
54  return
55end if
56! compute the spline coefficients
57call splinew(ni,wci,fi,cf)
58! evaluate spline at output points
59i=1
60do l=1,no
61  x=xo(l)
62  if (i.ge.ni) i=1
63  if (x.lt.xi(i)) goto 10
64  if (x.le.xi(i+1)) goto 30
65! binary search
6610 continue
67  i=1
68  j=ni
6920 continue
70  k=(i+j)/2
71  if (x.lt.xi(k)) then
72    j=k
73  else
74    i=k
75  end if
76  if (j.gt.i+1) goto 20
7730 continue
78  dx=x-xi(i)
79  fo(l)=fi(i)+dx*(cf(1,i)+dx*(cf(2,i)+dx*cf(3,i)))
80end do
81end subroutine
82!EOC
83
84