src/utils/interpolation.F90

!!
!! Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!
!! SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
!! research, please cite the appropriate papers when you publish your
!! work.  Good starting points are:
!!                                                                      
!! [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
!! [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
!! [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
!! [4]  Vardeman & Gezelter, in progress (2009).
!!
!!
!!  interpolation.F90
!!
!!  Created by Charles F. Vardeman II on 03 Apr 2006.
!!
!!  PURPOSE: Generic Spline interpolation routines. 
!!
!! @author Charles F. Vardeman II 
!! @version $Id: interpolation.F90,v 1.10 2009-11-25 20:02:05 gezelter Exp $


module interpolation
  use definitions
  use status
  implicit none
  PRIVATE

  type, public :: cubicSpline
     logical :: isUniform = .false.
     integer :: n = 0
     real(kind=dp) :: dx_i
     real (kind=dp), pointer,dimension(:)   :: x => null()
     real (kind=dp), pointer,dimension(:)   :: y => null()
     real (kind=dp), pointer,dimension(:)   :: b => null()
     real (kind=dp), pointer,dimension(:)   :: c => null()
     real (kind=dp), pointer,dimension(:)   :: d => null()
  end type cubicSpline

  public :: newSpline
  public :: deleteSpline
  public :: lookupSpline
  public :: lookupUniformSpline
  public :: lookupNonuniformSpline
  public :: lookupUniformSpline1d
  
contains
  

  subroutine newSpline(cs, x, y, isUniform)
    
    implicit none

    type (cubicSpline), intent(inout) :: cs
    real( kind = DP ), intent(in) :: x(:), y(:)
    real( kind = DP ) :: fp1, fpn, p
    REAL( KIND = DP), DIMENSION(size(x)-1) :: diff_y, H

    logical, intent(in) :: isUniform
    integer :: i, alloc_error, n, k

    alloc_error = 0

    if (cs%n .ne. 0) then
       call handleWarning("interpolation::newSpline", &
            "cubicSpline struct was already created")
       call deleteSpline(cs)
    end if

    ! make sure the sizes match

    n = size(x)
    
    if ( size(y) .ne. size(x) ) then
       call handleError("interpolation::newSpline", &
            "Array size mismatch")
    end if
    
    cs%n = n
    cs%isUniform = isUniform

    allocate(cs%x(n), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for x")
    endif

    allocate(cs%y(n), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for y")
    endif

    allocate(cs%b(n), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for b")
    endif

    allocate(cs%c(n), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for c")
    endif

    allocate(cs%d(n), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for d")
    endif

    ! make sure we are monotinically increasing in x:

    h = diff(x)
    if (any(h <= 0)) then
       call handleError("interpolation::newSpline", &
            "Negative dx interval found")
    end if

    ! load x and y values into the cubicSpline structure:

    do i = 1, n
       cs%x(i) = x(i)
       cs%y(i) = y(i)
    end do

    ! Calculate coefficients for the tridiagonal system: store
    ! sub-diagonal in B, diagonal in D, difference quotient in C.

    cs%b(1:n-1) = h
    diff_y = diff(y)
    cs%c(1:n-1) = diff_y / h

    if (n == 2) then
       ! Assume the derivatives at both endpoints are zero
       ! another assumption could be made to have a linear interpolant
       ! between the two points.  In that case, the b coefficients
       ! below would be diff_y(1)/h(1) and the c and d coefficients would
       ! both be zero.
       cs%b(1) = 0.0_dp
       cs%c(1) = -3.0_dp * (diff_y(1)/h(1))**2
       cs%d(1) = -2.0_dp * (diff_y(1)/h(1))**3
       cs%b(2) = cs%b(1)
       cs%c(2) = 0.0_dp
       cs%d(2) = 0.0_dp
       cs%dx_i = 1.0_dp / h(1)
      return
    end if

    cs%d(1) = 2.0_dp * cs%b(1)
    do i = 2, n-1
      cs%d(i) = 2.0_dp * (cs%b(i) + cs%b(i-1))
    end do
    cs%d(n) = 2.0_dp * cs%b(n-1)

    ! Calculate estimates for the end slopes using polynomials
    ! that interpolate the data nearest the end.
    
    fp1 = cs%c(1) - cs%b(1)*(cs%c(2) - cs%c(1))/(cs%b(1) + cs%b(2))
    if (n > 3) then
      fp1 = fp1 + cs%b(1)*((cs%b(1) + cs%b(2))*(cs%c(3) - cs%c(2))/ &
           (cs%b(2) + cs%b(3)) - cs%c(2) + cs%c(1))/(x(4) - x(1))
    end if
          
    fpn = cs%c(n-1) + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2))/(cs%b(n-2) + cs%b(n-1))
    if (n > 3) then
      fpn = fpn + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2) - (cs%b(n-2) + cs%b(n-1))* &
           (cs%c(n-2) - cs%c(n-3))/(cs%b(n-2) + cs%b(n-3)))/(x(n) - x(n-3))
    end if

    ! Calculate the right hand side and store it in C.

    cs%c(n) = 3.0_dp * (fpn - cs%c(n-1))
    do i = n-1,2,-1
      cs%c(i) = 3.0_dp * (cs%c(i) - cs%c(i-1))
    end do
    cs%c(1) = 3.0_dp * (cs%c(1) - fp1)

    ! Solve the tridiagonal system.

    do k = 2, n
      p = cs%b(k-1) / cs%d(k-1)
      cs%d(k) = cs%d(k) - p*cs%b(k-1)
      cs%c(k) = cs%c(k) - p*cs%c(k-1)
    end do
    cs%c(n) = cs%c(n) / cs%d(n)
    do k = n-1, 1, -1
      cs%c(k) = (cs%c(k) - cs%b(k) * cs%c(k+1)) / cs%d(k)
    end do

    ! Calculate the coefficients defining the spline.

    cs%d(1:n-1) = diff(cs%c) / (3.0_dp * h)
    cs%b(1:n-1) = diff_y / h - h * (cs%c(1:n-1) + h * cs%d(1:n-1))
    cs%b(n) = cs%b(n-1) + h(n-1) * (2.0_dp*cs%c(n-1) + h(n-1)*3.0_dp*cs%d(n-1))

    if (isUniform) then
       cs%dx_i = 1.0_dp / (x(2) - x(1))
    endif

    return
    
  contains
    
    function diff(v)
      ! Auxiliary function to compute the forward difference
      ! of data stored in a vector v.
      
      implicit none
      real (kind = dp), dimension(:), intent(in) :: v
      real (kind = dp), dimension(size(v)-1) :: diff
      
      integer :: n
      
      n = size(v)
      diff = v(2:n) - v(1:n-1)
      return
    end function diff
    
  end subroutine newSpline
       
  subroutine deleteSpline(this)

    type(cubicSpline) :: this
    
    if(associated(this%d)) then
       deallocate(this%d)
       this%d => null()
    end if
    if(associated(this%c)) then
       deallocate(this%c)
       this%c => null()
    end if
    if(associated(this%b)) then
       deallocate(this%b)
       this%b => null()
    end if
    if(associated(this%y)) then
       deallocate(this%y)
       this%y => null()
    end if
    if(associated(this%x)) then
       deallocate(this%x)
       this%x => null()
    end if
    
    this%n = 0
    
  end subroutine deleteSpline

  subroutine lookupNonuniformSpline(cs, xval, yval)
    
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) ::  dx
    integer :: i, j
    !
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.
    !
    j = cs%n - 1

    do i = 0, cs%n - 2

       if ( xval < cs%x(i+1) ) then
          j = i
          exit
       end if

    end do
    !
    !  Evaluate the cubic polynomial.
    !
    dx = xval - cs%x(j)
    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
    
    return
  end subroutine lookupNonuniformSpline

  subroutine lookupUniformSpline(cs, xval, yval)
    
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) ::  dx
    integer :: i, j
    !
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.
    
    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
    
    dx = xval - cs%x(j)
    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
    
    return
  end subroutine lookupUniformSpline

  subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
    
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval, dydx
    real( kind = DP ) :: dx
    integer :: i, j
    
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.


    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
    
    dx = xval - cs%x(j)
    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))

    dydx = cs%b(j) + dx*(2.0_dp * cs%c(j) + 3.0_dp * dx * cs%d(j))
       
    return
  end subroutine lookupUniformSpline1d

  subroutine lookupSpline(cs, xval, yval)

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(inout) :: xval
    real( kind = DP ), intent(inout) :: yval
    
    if (cs%isUniform) then
       call lookupUniformSpline(cs, xval, yval)
    else
       call lookupNonuniformSpline(cs, xval, yval)
    endif

    return
  end subroutine lookupSpline
  
end module interpolation
Revision:	1390
Committed:	Wed Nov 25 20:02:06 2009 UTC (15 years, 5 months ago) by gezelter
Original Path:	*trunk/src/utils/interpolation.F90*
File size:	10773 byte(s)
Log Message:	Almost all of the changes necessary to create OpenMD out of our old project (OOPSE-4)
#	User	Rev	Content
1	gezelter	931	!!
2			!! Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9	gezelter	1390	!! 1. Redistributions of source code must retain the above copyright
10	gezelter	931	!! notice, this list of conditions and the following disclaimer.
11			!!
12	gezelter	1390	!! 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	931	!! notice, this list of conditions and the following disclaimer in the
14			!! documentation and/or other materials provided with the
15			!! distribution.
16			!!
17			!! This software is provided "AS IS," without a warranty of any
18			!! kind. All express or implied conditions, representations and
19			!! warranties, including any implied warranty of merchantability,
20			!! fitness for a particular purpose or non-infringement, are hereby
21			!! excluded. The University of Notre Dame and its licensors shall not
22			!! be liable for any damages suffered by licensee as a result of
23			!! using, modifying or distributing the software or its
24			!! derivatives. In no event will the University of Notre Dame or its
25			!! licensors be liable for any lost revenue, profit or data, or for
26			!! direct, indirect, special, consequential, incidental or punitive
27			!! damages, however caused and regardless of the theory of liability,
28			!! arising out of the use of or inability to use software, even if the
29			!! University of Notre Dame has been advised of the possibility of
30			!! such damages.
31			!!
32	gezelter	1390	!! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			!! research, please cite the appropriate papers when you publish your
34			!! work. Good starting points are:
35			!!
36			!! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			!! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			!! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			!! [4] Vardeman & Gezelter, in progress (2009).
40	gezelter	931	!!
41	gezelter	1390	!!
42	gezelter	931	!! interpolation.F90
43			!!
44			!! Created by Charles F. Vardeman II on 03 Apr 2006.
45			!!
46	gezelter	939	!! PURPOSE: Generic Spline interpolation routines.
47	gezelter	931	!!
48			!! @author Charles F. Vardeman II
49	gezelter	1390	!! @version $Id: interpolation.F90,v 1.10 2009-11-25 20:02:05 gezelter Exp $
50	gezelter	931
51
52	gezelter	938	module interpolation
53	gezelter	931	use definitions
54			use status
55			implicit none
56			PRIVATE
57
58			type, public :: cubicSpline
59	gezelter	934	logical :: isUniform = .false.
60	gezelter	939	integer :: n = 0
61	gezelter	931	real(kind=dp) :: dx_i
62			real (kind=dp), pointer,dimension(:) :: x => null()
63	gezelter	939	real (kind=dp), pointer,dimension(:) :: y => null()
64			real (kind=dp), pointer,dimension(:) :: b => null()
65			real (kind=dp), pointer,dimension(:) :: c => null()
66			real (kind=dp), pointer,dimension(:) :: d => null()
67	gezelter	931	end type cubicSpline
68
69	gezelter	934	public :: newSpline
70	gezelter	931	public :: deleteSpline
71	gezelter	938	public :: lookupSpline
72			public :: lookupUniformSpline
73			public :: lookupNonuniformSpline
74			public :: lookupUniformSpline1d
75	gezelter	934
76	gezelter	931	contains
77	gezelter	934
78	gezelter	931
79	gezelter	939	subroutine newSpline(cs, x, y, isUniform)
80	gezelter	934
81	gezelter	931	implicit none
82
83			type (cubicSpline), intent(inout) :: cs
84			real( kind = DP ), intent(in) :: x(:), y(:)
85	gezelter	939	real( kind = DP ) :: fp1, fpn, p
86			REAL( KIND = DP), DIMENSION(size(x)-1) :: diff_y, H
87
88	gezelter	934	logical, intent(in) :: isUniform
89	gezelter	939	integer :: i, alloc_error, n, k
90	gezelter	931
91			alloc_error = 0
92
93	gezelter	939	if (cs%n .ne. 0) then
94	gezelter	934	call handleWarning("interpolation::newSpline", &
95			"cubicSpline struct was already created")
96	gezelter	931	call deleteSpline(cs)
97			end if
98
99			! make sure the sizes match
100
101	gezelter	939	n = size(x)
102
103			if ( size(y) .ne. size(x) ) then
104	gezelter	934	call handleError("interpolation::newSpline", &
105	gezelter	931	"Array size mismatch")
106			end if
107	gezelter	934
108	gezelter	939	cs%n = n
109	gezelter	934	cs%isUniform = isUniform
110	gezelter	931
111	gezelter	939	allocate(cs%x(n), stat=alloc_error)
112	gezelter	931	if(alloc_error .ne. 0) then
113	gezelter	934	call handleError("interpolation::newSpline", &
114	gezelter	931	"Error in allocating storage for x")
115			endif
116
117	gezelter	939	allocate(cs%y(n), stat=alloc_error)
118	gezelter	931	if(alloc_error .ne. 0) then
119	gezelter	934	call handleError("interpolation::newSpline", &
120	gezelter	939	"Error in allocating storage for y")
121			endif
122
123			allocate(cs%b(n), stat=alloc_error)
124			if(alloc_error .ne. 0) then
125			call handleError("interpolation::newSpline", &
126			"Error in allocating storage for b")
127			endif
128
129			allocate(cs%c(n), stat=alloc_error)
130			if(alloc_error .ne. 0) then
131			call handleError("interpolation::newSpline", &
132	gezelter	931	"Error in allocating storage for c")
133			endif
134	gezelter	939
135			allocate(cs%d(n), stat=alloc_error)
136			if(alloc_error .ne. 0) then
137			call handleError("interpolation::newSpline", &
138			"Error in allocating storage for d")
139			endif
140
141			! make sure we are monotinically increasing in x:
142
143			h = diff(x)
144			if (any(h <= 0)) then
145			call handleError("interpolation::newSpline", &
146			"Negative dx interval found")
147			end if
148
149			! load x and y values into the cubicSpline structure:
150
151			do i = 1, n
152	gezelter	931	cs%x(i) = x(i)
153	gezelter	939	cs%y(i) = y(i)
154			end do
155	gezelter	931
156	gezelter	939	! Calculate coefficients for the tridiagonal system: store
157			! sub-diagonal in B, diagonal in D, difference quotient in C.
158	gezelter	931
159	gezelter	939	cs%b(1:n-1) = h
160			diff_y = diff(y)
161			cs%c(1:n-1) = diff_y / h
162
163			if (n == 2) then
164			! Assume the derivatives at both endpoints are zero
165			! another assumption could be made to have a linear interpolant
166			! between the two points. In that case, the b coefficients
167			! below would be diff_y(1)/h(1) and the c and d coefficients would
168			! both be zero.
169			cs%b(1) = 0.0_dp
170			cs%c(1) = -3.0_dp * (diff_y(1)/h(1))**2
171			cs%d(1) = -2.0_dp * (diff_y(1)/h(1))**3
172			cs%b(2) = cs%b(1)
173			cs%c(2) = 0.0_dp
174			cs%d(2) = 0.0_dp
175			cs%dx_i = 1.0_dp / h(1)
176			return
177			end if
178
179			cs%d(1) = 2.0_dp * cs%b(1)
180			do i = 2, n-1
181			cs%d(i) = 2.0_dp * (cs%b(i) + cs%b(i-1))
182			end do
183			cs%d(n) = 2.0_dp * cs%b(n-1)
184
185			! Calculate estimates for the end slopes using polynomials
186			! that interpolate the data nearest the end.
187	chrisfen	933
188	gezelter	939	fp1 = cs%c(1) - cs%b(1)*(cs%c(2) - cs%c(1))/(cs%b(1) + cs%b(2))
189			if (n > 3) then
190			fp1 = fp1 + cs%b(1)((cs%b(1) + cs%b(2))(cs%c(3) - cs%c(2))/ &
191			(cs%b(2) + cs%b(3)) - cs%c(2) + cs%c(1))/(x(4) - x(1))
192			end if
193
194			fpn = cs%c(n-1) + cs%b(n-1)*(cs%c(n-1) - cs%c(n-2))/(cs%b(n-2) + cs%b(n-1))
195			if (n > 3) then
196			fpn = fpn + cs%b(n-1)(cs%c(n-1) - cs%c(n-2) - (cs%b(n-2) + cs%b(n-1)) &
197			(cs%c(n-2) - cs%c(n-3))/(cs%b(n-2) + cs%b(n-3)))/(x(n) - x(n-3))
198			end if
199	gezelter	934
200	gezelter	939	! Calculate the right hand side and store it in C.
201
202			cs%c(n) = 3.0_dp * (fpn - cs%c(n-1))
203			do i = n-1,2,-1
204			cs%c(i) = 3.0_dp * (cs%c(i) - cs%c(i-1))
205	gezelter	931	end do
206	gezelter	939	cs%c(1) = 3.0_dp * (cs%c(1) - fp1)
207	gezelter	934
208	gezelter	939	! Solve the tridiagonal system.
209
210			do k = 2, n
211			p = cs%b(k-1) / cs%d(k-1)
212			cs%d(k) = cs%d(k) - p*cs%b(k-1)
213			cs%c(k) = cs%c(k) - p*cs%c(k-1)
214	gezelter	931	end do
215	gezelter	939	cs%c(n) = cs%c(n) / cs%d(n)
216			do k = n-1, 1, -1
217			cs%c(k) = (cs%c(k) - cs%b(k) * cs%c(k+1)) / cs%d(k)
218	gezelter	931	end do
219
220	gezelter	939	! Calculate the coefficients defining the spline.
221	gezelter	931
222	gezelter	939	cs%d(1:n-1) = diff(cs%c) / (3.0_dp * h)
223			cs%b(1:n-1) = diff_y / h - h * (cs%c(1:n-1) + h * cs%d(1:n-1))
224			cs%b(n) = cs%b(n-1) + h(n-1) * (2.0_dpcs%c(n-1) + h(n-1)3.0_dp*cs%d(n-1))
225	gezelter	934
226	gezelter	939	if (isUniform) then
227	gezelter	960	cs%dx_i = 1.0_dp / (x(2) - x(1))
228	gezelter	939	endif
229
230	gezelter	931	return
231	gezelter	939
232			contains
233
234			function diff(v)
235			! Auxiliary function to compute the forward difference
236			! of data stored in a vector v.
237
238			implicit none
239			real (kind = dp), dimension(:), intent(in) :: v
240			real (kind = dp), dimension(size(v)-1) :: diff
241
242			integer :: n
243
244			n = size(v)
245			diff = v(2:n) - v(1:n-1)
246			return
247			end function diff
248
249	gezelter	935	end subroutine newSpline
250	gezelter	939
251	gezelter	931	subroutine deleteSpline(this)
252
253			type(cubicSpline) :: this
254
255	gezelter	983	if(associated(this%d)) then
256			deallocate(this%d)
257			this%d => null()
258	gezelter	931	end if
259			if(associated(this%c)) then
260			deallocate(this%c)
261			this%c => null()
262			end if
263	gezelter	983	if(associated(this%b)) then
264			deallocate(this%b)
265			this%b => null()
266			end if
267			if(associated(this%y)) then
268			deallocate(this%y)
269			this%y => null()
270			end if
271			if(associated(this%x)) then
272			deallocate(this%x)
273			this%x => null()
274			end if
275	gezelter	931
276	gezelter	939	this%n = 0
277	gezelter	931
278			end subroutine deleteSpline
279
280	gezelter	938	subroutine lookupNonuniformSpline(cs, xval, yval)
281	gezelter	931
282			implicit none
283
284			type (cubicSpline), intent(in) :: cs
285			real( kind = DP ), intent(in) :: xval
286			real( kind = DP ), intent(out) :: yval
287	gezelter	939	real( kind = DP ) :: dx
288	gezelter	931	integer :: i, j
289			!
290			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
291			! or is nearest to xval.
292			!
293	gezelter	939	j = cs%n - 1
294	gezelter	931
295	gezelter	939	do i = 0, cs%n - 2
296	gezelter	931
297			if ( xval < cs%x(i+1) ) then
298			j = i
299			exit
300			end if
301
302			end do
303			!
304			! Evaluate the cubic polynomial.
305			!
306			dx = xval - cs%x(j)
307	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
308	gezelter	931
309			return
310	gezelter	938	end subroutine lookupNonuniformSpline
311	gezelter	931
312	gezelter	938	subroutine lookupUniformSpline(cs, xval, yval)
313	gezelter	931
314			implicit none
315
316			type (cubicSpline), intent(in) :: cs
317			real( kind = DP ), intent(in) :: xval
318			real( kind = DP ), intent(out) :: yval
319	gezelter	939	real( kind = DP ) :: dx
320	gezelter	931	integer :: i, j
321			!
322			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
323			! or is nearest to xval.
324	gezelter	939
325	gezelter	960	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
326	gezelter	939
327	gezelter	931	dx = xval - cs%x(j)
328	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
329	gezelter	931
330			return
331	gezelter	938	end subroutine lookupUniformSpline
332	gezelter	934
333	gezelter	938	subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
334
335			implicit none
336	gezelter	934
337			type (cubicSpline), intent(in) :: cs
338	gezelter	938	real( kind = DP ), intent(in) :: xval
339			real( kind = DP ), intent(out) :: yval, dydx
340	gezelter	939	real( kind = DP ) :: dx
341	gezelter	938	integer :: i, j
342
343			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
344			! or is nearest to xval.
345
346
347	gezelter	960	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
348	gezelter	939
349	gezelter	938	dx = xval - cs%x(j)
350	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
351	gezelter	938
352	gezelter	960	dydx = cs%b(j) + dx(2.0_dp cs%c(j) + 3.0_dp * dx * cs%d(j))
353	gezelter	938
354			return
355			end subroutine lookupUniformSpline1d
356
357			subroutine lookupSpline(cs, xval, yval)
358
359			type (cubicSpline), intent(in) :: cs
360	gezelter	934	real( kind = DP ), intent(inout) :: xval
361			real( kind = DP ), intent(inout) :: yval
362
363			if (cs%isUniform) then
364	gezelter	938	call lookupUniformSpline(cs, xval, yval)
365	gezelter	934	else
366	gezelter	938	call lookupNonuniformSpline(cs, xval, yval)
367	gezelter	934	endif
368
369			return
370	gezelter	938	end subroutine lookupSpline
371	gezelter	931
372	gezelter	938	end module interpolation