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. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. 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.
!!
!!
!!  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.7 2006-04-20 18:24:24 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.0d0 / (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%x)) then
       deallocate(this%x)
       this%x => null()
    end if
    if(associated(this%c)) then
       deallocate(this%c)
       this%c => 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, idint((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, idint((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.0d0 * cs%c(j) + 3.0d0 * 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:	939
Committed:	Thu Apr 20 18:24:24 2006 UTC (19 years ago) by gezelter
File size:	10477 byte(s)
Log Message:	Complete rewrite of spline code and everything that uses it.
#	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			!! 1. Acknowledgement of the program authors must be made in any
10			!! publication of scientific results based in part on use of the
11			!! program. An acceptable form of acknowledgement is citation of
12			!! the article in which the program was described (Matthew
13			!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			!! Parallel Simulation Engine for Molecular Dynamics,"
16			!! J. Comput. Chem. 26, pp. 252-271 (2005))
17			!!
18			!! 2. Redistributions of source code must retain the above copyright
19			!! notice, this list of conditions and the following disclaimer.
20			!!
21			!! 3. Redistributions in binary form must reproduce the above copyright
22			!! notice, this list of conditions and the following disclaimer in the
23			!! documentation and/or other materials provided with the
24			!! distribution.
25			!!
26			!! This software is provided "AS IS," without a warranty of any
27			!! kind. All express or implied conditions, representations and
28			!! warranties, including any implied warranty of merchantability,
29			!! fitness for a particular purpose or non-infringement, are hereby
30			!! excluded. The University of Notre Dame and its licensors shall not
31			!! be liable for any damages suffered by licensee as a result of
32			!! using, modifying or distributing the software or its
33			!! derivatives. In no event will the University of Notre Dame or its
34			!! licensors be liable for any lost revenue, profit or data, or for
35			!! direct, indirect, special, consequential, incidental or punitive
36			!! damages, however caused and regardless of the theory of liability,
37			!! arising out of the use of or inability to use software, even if the
38			!! University of Notre Dame has been advised of the possibility of
39			!! such damages.
40			!!
41			!!
42			!! 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	939	!! @version $Id: interpolation.F90,v 1.7 2006-04-20 18:24:24 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			cs%dx_i = 1.0d0 / (x(2) - x(1))
228			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			if(associated(this%x)) then
256			deallocate(this%x)
257			this%x => null()
258			end if
259			if(associated(this%c)) then
260			deallocate(this%c)
261			this%c => null()
262			end if
263
264	gezelter	939	this%n = 0
265	gezelter	931
266			end subroutine deleteSpline
267
268	gezelter	938	subroutine lookupNonuniformSpline(cs, xval, yval)
269	gezelter	931
270			implicit none
271
272			type (cubicSpline), intent(in) :: cs
273			real( kind = DP ), intent(in) :: xval
274			real( kind = DP ), intent(out) :: yval
275	gezelter	939	real( kind = DP ) :: dx
276	gezelter	931	integer :: i, j
277			!
278			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
279			! or is nearest to xval.
280			!
281	gezelter	939	j = cs%n - 1
282	gezelter	931
283	gezelter	939	do i = 0, cs%n - 2
284	gezelter	931
285			if ( xval < cs%x(i+1) ) then
286			j = i
287			exit
288			end if
289
290			end do
291			!
292			! Evaluate the cubic polynomial.
293			!
294			dx = xval - cs%x(j)
295	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
296	gezelter	931
297			return
298	gezelter	938	end subroutine lookupNonuniformSpline
299	gezelter	931
300	gezelter	938	subroutine lookupUniformSpline(cs, xval, yval)
301	gezelter	931
302			implicit none
303
304			type (cubicSpline), intent(in) :: cs
305			real( kind = DP ), intent(in) :: xval
306			real( kind = DP ), intent(out) :: yval
307	gezelter	939	real( kind = DP ) :: dx
308	gezelter	931	integer :: i, j
309			!
310			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
311			! or is nearest to xval.
312	gezelter	939
313			j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
314
315	gezelter	931	dx = xval - cs%x(j)
316	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
317	gezelter	931
318			return
319	gezelter	938	end subroutine lookupUniformSpline
320	gezelter	934
321	gezelter	938	subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
322
323			implicit none
324	gezelter	934
325			type (cubicSpline), intent(in) :: cs
326	gezelter	938	real( kind = DP ), intent(in) :: xval
327			real( kind = DP ), intent(out) :: yval, dydx
328	gezelter	939	real( kind = DP ) :: dx
329	gezelter	938	integer :: i, j
330
331			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
332			! or is nearest to xval.
333
334
335	gezelter	939	j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
336
337	gezelter	938	dx = xval - cs%x(j)
338	gezelter	939	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
339	gezelter	938
340	gezelter	939	dydx = cs%b(j) + dx(2.0d0 cs%c(j) + 3.0d0 * dx * cs%d(j))
341	gezelter	938
342			return
343			end subroutine lookupUniformSpline1d
344
345			subroutine lookupSpline(cs, xval, yval)
346
347			type (cubicSpline), intent(in) :: cs
348	gezelter	934	real( kind = DP ), intent(inout) :: xval
349			real( kind = DP ), intent(inout) :: yval
350
351			if (cs%isUniform) then
352	gezelter	938	call lookupUniformSpline(cs, xval, yval)
353	gezelter	934	else
354	gezelter	938	call lookupNonuniformSpline(cs, xval, yval)
355	gezelter	934	endif
356
357			return
358	gezelter	938	end subroutine lookupSpline
359	gezelter	931
360	gezelter	938	end module interpolation