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.9 2006-06-06 17:43:28 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:	2802
Committed:	Tue Jun 6 17:43:28 2006 UTC (19 years, 4 months ago) by gezelter
File size:	10755 byte(s)
Log Message:	testing GB, removing CM drift in Langevin Dynamics, fixing a memory leak, adding a visitor
#	User	Rev	Content
1	gezelter	2708	!!
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	2722	!! PURPOSE: Generic Spline interpolation routines.
47	gezelter	2708	!!
48			!! @author Charles F. Vardeman II
49	gezelter	2802	!! @version $Id: interpolation.F90,v 1.9 2006-06-06 17:43:28 gezelter Exp $
50	gezelter	2708
51
52	gezelter	2717	module interpolation
53	gezelter	2708	use definitions
54			use status
55			implicit none
56			PRIVATE
57
58			type, public :: cubicSpline
59	gezelter	2711	logical :: isUniform = .false.
60	gezelter	2722	integer :: n = 0
61	gezelter	2708	real(kind=dp) :: dx_i
62			real (kind=dp), pointer,dimension(:) :: x => null()
63	gezelter	2722	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	2708	end type cubicSpline
68
69	gezelter	2711	public :: newSpline
70	gezelter	2708	public :: deleteSpline
71	gezelter	2717	public :: lookupSpline
72			public :: lookupUniformSpline
73			public :: lookupNonuniformSpline
74			public :: lookupUniformSpline1d
75	gezelter	2711
76	gezelter	2708	contains
77	gezelter	2711
78	gezelter	2708
79	gezelter	2722	subroutine newSpline(cs, x, y, isUniform)
80	gezelter	2711
81	gezelter	2708	implicit none
82
83			type (cubicSpline), intent(inout) :: cs
84			real( kind = DP ), intent(in) :: x(:), y(:)
85	gezelter	2722	real( kind = DP ) :: fp1, fpn, p
86			REAL( KIND = DP), DIMENSION(size(x)-1) :: diff_y, H
87
88	gezelter	2711	logical, intent(in) :: isUniform
89	gezelter	2722	integer :: i, alloc_error, n, k
90	gezelter	2708
91			alloc_error = 0
92
93	gezelter	2722	if (cs%n .ne. 0) then
94	gezelter	2711	call handleWarning("interpolation::newSpline", &
95			"cubicSpline struct was already created")
96	gezelter	2708	call deleteSpline(cs)
97			end if
98
99			! make sure the sizes match
100
101	gezelter	2722	n = size(x)
102
103			if ( size(y) .ne. size(x) ) then
104	gezelter	2711	call handleError("interpolation::newSpline", &
105	gezelter	2708	"Array size mismatch")
106			end if
107	gezelter	2711
108	gezelter	2722	cs%n = n
109	gezelter	2711	cs%isUniform = isUniform
110	gezelter	2708
111	gezelter	2722	allocate(cs%x(n), stat=alloc_error)
112	gezelter	2708	if(alloc_error .ne. 0) then
113	gezelter	2711	call handleError("interpolation::newSpline", &
114	gezelter	2708	"Error in allocating storage for x")
115			endif
116
117	gezelter	2722	allocate(cs%y(n), stat=alloc_error)
118	gezelter	2708	if(alloc_error .ne. 0) then
119	gezelter	2711	call handleError("interpolation::newSpline", &
120	gezelter	2722	"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	2708	"Error in allocating storage for c")
133			endif
134	gezelter	2722
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	2708	cs%x(i) = x(i)
153	gezelter	2722	cs%y(i) = y(i)
154			end do
155	gezelter	2708
156	gezelter	2722	! Calculate coefficients for the tridiagonal system: store
157			! sub-diagonal in B, diagonal in D, difference quotient in C.
158	gezelter	2708
159	gezelter	2722	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	2710
188	gezelter	2722	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	2711
200	gezelter	2722	! 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	2708	end do
206	gezelter	2722	cs%c(1) = 3.0_dp * (cs%c(1) - fp1)
207	gezelter	2711
208	gezelter	2722	! 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	2708	end do
215	gezelter	2722	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	2708	end do
219
220	gezelter	2722	! Calculate the coefficients defining the spline.
221	gezelter	2708
222	gezelter	2722	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	2711
226	gezelter	2722	if (isUniform) then
227	gezelter	2756	cs%dx_i = 1.0_dp / (x(2) - x(1))
228	gezelter	2722	endif
229
230	gezelter	2708	return
231	gezelter	2722
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	2712	end subroutine newSpline
250	gezelter	2722
251	gezelter	2708	subroutine deleteSpline(this)
252
253			type(cubicSpline) :: this
254
255	gezelter	2802	if(associated(this%d)) then
256			deallocate(this%d)
257			this%d => null()
258	gezelter	2708	end if
259			if(associated(this%c)) then
260			deallocate(this%c)
261			this%c => null()
262			end if
263	gezelter	2802	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	2708
276	gezelter	2722	this%n = 0
277	gezelter	2708
278			end subroutine deleteSpline
279
280	gezelter	2717	subroutine lookupNonuniformSpline(cs, xval, yval)
281	gezelter	2708
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	2722	real( kind = DP ) :: dx
288	gezelter	2708	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	2722	j = cs%n - 1
294	gezelter	2708
295	gezelter	2722	do i = 0, cs%n - 2
296	gezelter	2708
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	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
308	gezelter	2708
309			return
310	gezelter	2717	end subroutine lookupNonuniformSpline
311	gezelter	2708
312	gezelter	2717	subroutine lookupUniformSpline(cs, xval, yval)
313	gezelter	2708
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	2722	real( kind = DP ) :: dx
320	gezelter	2708	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	2722
325	gezelter	2756	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
326	gezelter	2722
327	gezelter	2708	dx = xval - cs%x(j)
328	gezelter	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
329	gezelter	2708
330			return
331	gezelter	2717	end subroutine lookupUniformSpline
332	gezelter	2711
333	gezelter	2717	subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
334
335			implicit none
336	gezelter	2711
337			type (cubicSpline), intent(in) :: cs
338	gezelter	2717	real( kind = DP ), intent(in) :: xval
339			real( kind = DP ), intent(out) :: yval, dydx
340	gezelter	2722	real( kind = DP ) :: dx
341	gezelter	2717	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	2756	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
348	gezelter	2722
349	gezelter	2717	dx = xval - cs%x(j)
350	gezelter	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
351	gezelter	2717
352	gezelter	2756	dydx = cs%b(j) + dx(2.0_dp cs%c(j) + 3.0_dp * dx * cs%d(j))
353	gezelter	2717
354			return
355			end subroutine lookupUniformSpline1d
356
357			subroutine lookupSpline(cs, xval, yval)
358
359			type (cubicSpline), intent(in) :: cs
360	gezelter	2711	real( kind = DP ), intent(inout) :: xval
361			real( kind = DP ), intent(inout) :: yval
362
363			if (cs%isUniform) then
364	gezelter	2717	call lookupUniformSpline(cs, xval, yval)
365	gezelter	2711	else
366	gezelter	2717	call lookupNonuniformSpline(cs, xval, yval)
367	gezelter	2711	endif
368
369			return
370	gezelter	2717	end subroutine lookupSpline
371	gezelter	2708
372	gezelter	2717	end module interpolation