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. These routines 
!!           assume that we are on a uniform grid for precomputation of 
!!           spline parameters.
!!
!! @author Charles F. Vardeman II 
!! @version $Id: interpolation.F90,v 1.6 2006-04-17 21:49:12 gezelter Exp $


module interpolation
  use definitions
  use status
  implicit none
  PRIVATE

  character(len = statusMsgSize) :: errMSG

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

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

  subroutine newSpline(cs, x, y, yp1, ypn, isUniform)
    
    !************************************************************************
    !
    ! newSpline solves for slopes defining a cubic spline.
    !
    !  Discussion:
    !
    !    A tridiagonal linear system for the unknown slopes S(I) of
    !    F at x(I), I=1,..., N, is generated and then solved by Gauss
    !    elimination, with S(I) ending up in cs%C(2,I), for all I.
    !
    !  Reference:
    !
    !    Carl DeBoor,
    !    A Practical Guide to Splines,
    !    Springer Verlag.
    !
    !  Parameters:
    !
    !    Input, real x(N), the abscissas or X values of
    !    the data points.  The entries of x are assumed to be
    !    strictly increasing.
    !
    !    Input, real y(I), contains the function value at x(I) for 
    !      I = 1, N.
    !
    !    Input, real yp1 contains the slope at x(1) 
    !    Input, real ypn contains the slope at x(N)
    !
    !    On output, the slopes at x(I) have been stored in 
    !               cs%C(2,I), for I = 1 to N.

    implicit none

    type (cubicSpline), intent(inout) :: cs
    real( kind = DP ), intent(in) :: x(:), y(:)
    real( kind = DP ), intent(in) :: yp1, ypn
    logical, intent(in) :: isUniform
    real( kind = DP ) :: g, divdif1, divdif3, dx
    integer :: i, alloc_error, np

    alloc_error = 0

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

    ! make sure the sizes match

    np = size(x)

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

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

    allocate(cs%c(4,np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "Error in allocating storage for c")
    endif
       
    do i = 1, np
       cs%x(i) = x(i)
       cs%c(1,i) = y(i)       
    enddo

    ! Set the first derivative of the function to the second coefficient of 
    ! each of the endpoints

    cs%c(2,1) = yp1
    cs%c(2,np) = ypn
    
    !
    !  Set up the right hand side of the linear system.
    !

    do i = 2, cs%np - 1
       cs%c(2,i) = 3.0_DP * ( &
            (x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
            (x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
    end do

    !
    !  Set the diagonal coefficients.
    !
    cs%c(4,1) = 1.0_DP
    do i = 2, cs%np - 1 
       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
    end do
    cs%c(4,cs%np) = 1.0_DP
    !
    !  Set the off-diagonal coefficients.
    !
    cs%c(3,1) = 0.0_DP
    do i = 2, cs%np
       cs%c(3,i) = x(i) - x(i-1)
    end do
    !
    !  Forward elimination.
    !
    do i = 2, cs%np - 1
       g = -cs%c(3,i+1) / cs%c(4,i-1)
       cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
       cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
    end do
    !
    !  Back substitution for the interior slopes.
    !
    do i = cs%np - 1, 2, -1
       cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
    end do
    !
    !  Now compute the quadratic and cubic coefficients used in the 
    !  piecewise polynomial representation.
    !
    do i = 1, cs%np - 1
       dx = x(i+1) - x(i)
       divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
       divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
       cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
       cs%c(4,i) = divdif3 / ( dx * dx )
    end do

    cs%c(3,cs%np) = 0.0_DP
    cs%c(4,cs%np) = 0.0_DP

    cs%dx_i = 1.0_DP / dx

    return
  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%np = 0
    
  end subroutine deleteSpline

  subroutine lookupNonuniformSpline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookupNonuniformSpline evaluates a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    newSpline must be called first, to set up the
    !    spline data from the raw function and derivative data.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Reference:
    !
    !    Conte and de Boor,
    !    Algorithm PCUBIC,
    !    Elementary Numerical Analysis, 
    !    1973, page 234.
    !
    !  Parameters:
    !
    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%np - 1

    do i = 1, cs%np - 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%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
    
    return
  end subroutine lookupNonuniformSpline

  subroutine lookupUniformSpline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookupUniformSpline evaluates a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    newSpline must be called first, to set up the
    !    spline data from the raw function and derivative data.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Reference:
    !
    !    Conte and de Boor,
    !    Algorithm PCUBIC,
    !    Elementary Numerical Analysis, 
    !    1973, page 234.
    !
    !  Parameters:
    !
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) :: a, b, c, d, 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%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))

    dx = xval - cs%x(j)

    a = cs%c(1,j)
    b = cs%c(2,j)
    c = cs%c(3,j)
    d = cs%c(4,j)

    yval = c + dx * d
    yval = b + dx * yval  
    yval = a + dx * yval
    
    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 ) :: a, b, c, d, 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%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))

    dx = xval - cs%x(j)

    a = cs%c(1,j)
    b = cs%c(2,j)
    c = cs%c(3,j)
    d = cs%c(4,j)

    yval = c + dx * d
    yval = b + dx * yval  
    yval = a + dx * yval

    dydx = 2.0d0 * c + 3.0d0 * d * dx
    dydx = b + dx * dydx
       
    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:	2717
Committed:	Mon Apr 17 21:49:12 2006 UTC (19 years, 6 months ago) by gezelter
File size:	10805 byte(s)
Log Message:	Many performance improvements
#	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	2711	!! PURPOSE: Generic Spline interpolation routines. These routines
47			!! assume that we are on a uniform grid for precomputation of
48			!! spline parameters.
49	gezelter	2708	!!
50			!! @author Charles F. Vardeman II
51	gezelter	2717	!! @version $Id: interpolation.F90,v 1.6 2006-04-17 21:49:12 gezelter Exp $
52	gezelter	2708
53
54	gezelter	2717	module interpolation
55	gezelter	2708	use definitions
56			use status
57			implicit none
58			PRIVATE
59
60			character(len = statusMsgSize) :: errMSG
61
62			type, public :: cubicSpline
63	gezelter	2711	logical :: isUniform = .false.
64	gezelter	2708	integer :: np = 0
65			real(kind=dp) :: dx_i
66			real (kind=dp), pointer,dimension(:) :: x => null()
67	gezelter	2709	real (kind=dp), pointer,dimension(:,:) :: c => null()
68	gezelter	2708	end type cubicSpline
69
70	gezelter	2711	public :: newSpline
71	gezelter	2708	public :: deleteSpline
72	gezelter	2717	public :: lookupSpline
73			public :: lookupUniformSpline
74			public :: lookupNonuniformSpline
75			public :: lookupUniformSpline1d
76	gezelter	2711
77	gezelter	2708	contains
78	gezelter	2711
79	gezelter	2708
80	gezelter	2711	subroutine newSpline(cs, x, y, yp1, ypn, isUniform)
81
82	gezelter	2708	!************************************************************************
83			!
84	gezelter	2711	! newSpline solves for slopes defining a cubic spline.
85	gezelter	2708	!
86			! Discussion:
87			!
88			! A tridiagonal linear system for the unknown slopes S(I) of
89			! F at x(I), I=1,..., N, is generated and then solved by Gauss
90			! elimination, with S(I) ending up in cs%C(2,I), for all I.
91			!
92			! Reference:
93			!
94			! Carl DeBoor,
95			! A Practical Guide to Splines,
96			! Springer Verlag.
97			!
98			! Parameters:
99			!
100			! Input, real x(N), the abscissas or X values of
101	chrisfen	2710	! the data points. The entries of x are assumed to be
102	gezelter	2708	! strictly increasing.
103			!
104			! Input, real y(I), contains the function value at x(I) for
105			! I = 1, N.
106			!
107	gezelter	2711	! Input, real yp1 contains the slope at x(1)
108			! Input, real ypn contains the slope at x(N)
109	gezelter	2708	!
110	gezelter	2711	! On output, the slopes at x(I) have been stored in
111			! cs%C(2,I), for I = 1 to N.
112	gezelter	2708
113			implicit none
114
115			type (cubicSpline), intent(inout) :: cs
116			real( kind = DP ), intent(in) :: x(:), y(:)
117			real( kind = DP ), intent(in) :: yp1, ypn
118	gezelter	2711	logical, intent(in) :: isUniform
119	gezelter	2708	real( kind = DP ) :: g, divdif1, divdif3, dx
120			integer :: i, alloc_error, np
121
122			alloc_error = 0
123
124			if (cs%np .ne. 0) then
125	gezelter	2711	call handleWarning("interpolation::newSpline", &
126			"cubicSpline struct was already created")
127	gezelter	2708	call deleteSpline(cs)
128			end if
129
130			! make sure the sizes match
131
132	gezelter	2711	np = size(x)
133
134			if ( size(y) .ne. np ) then
135			call handleError("interpolation::newSpline", &
136	gezelter	2708	"Array size mismatch")
137			end if
138	gezelter	2711
139	gezelter	2708	cs%np = np
140	gezelter	2711	cs%isUniform = isUniform
141	gezelter	2708
142			allocate(cs%x(np), stat=alloc_error)
143			if(alloc_error .ne. 0) then
144	gezelter	2711	call handleError("interpolation::newSpline", &
145	gezelter	2708	"Error in allocating storage for x")
146			endif
147
148			allocate(cs%c(4,np), stat=alloc_error)
149			if(alloc_error .ne. 0) then
150	gezelter	2711	call handleError("interpolation::newSpline", &
151	gezelter	2708	"Error in allocating storage for c")
152			endif
153
154			do i = 1, np
155			cs%x(i) = x(i)
156			cs%c(1,i) = y(i)
157			enddo
158
159	chrisfen	2710	! Set the first derivative of the function to the second coefficient of
160			! each of the endpoints
161	gezelter	2708
162	chrisfen	2710	cs%c(2,1) = yp1
163			cs%c(2,np) = ypn
164
165	gezelter	2708	!
166			! Set up the right hand side of the linear system.
167			!
168	gezelter	2711
169	gezelter	2708	do i = 2, cs%np - 1
170			cs%c(2,i) = 3.0_DP * ( &
171			(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
172			(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
173			end do
174	gezelter	2711
175	gezelter	2708	!
176			! Set the diagonal coefficients.
177			!
178			cs%c(4,1) = 1.0_DP
179			do i = 2, cs%np - 1
180			cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
181			end do
182	gezelter	2709	cs%c(4,cs%np) = 1.0_DP
183	gezelter	2708	!
184			! Set the off-diagonal coefficients.
185			!
186			cs%c(3,1) = 0.0_DP
187			do i = 2, cs%np
188			cs%c(3,i) = x(i) - x(i-1)
189			end do
190			!
191			! Forward elimination.
192			!
193			do i = 2, cs%np - 1
194			g = -cs%c(3,i+1) / cs%c(4,i-1)
195			cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
196			cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
197			end do
198			!
199			! Back substitution for the interior slopes.
200			!
201			do i = cs%np - 1, 2, -1
202			cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
203			end do
204			!
205			! Now compute the quadratic and cubic coefficients used in the
206			! piecewise polynomial representation.
207			!
208			do i = 1, cs%np - 1
209			dx = x(i+1) - x(i)
210			divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
211			divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
212			cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
213			cs%c(4,i) = divdif3 / ( dx * dx )
214			end do
215
216	gezelter	2709	cs%c(3,cs%np) = 0.0_DP
217			cs%c(4,cs%np) = 0.0_DP
218	gezelter	2708
219	gezelter	2709	cs%dx_i = 1.0_DP / dx
220	gezelter	2711
221	gezelter	2708	return
222	gezelter	2712	end subroutine newSpline
223	gezelter	2708
224			subroutine deleteSpline(this)
225
226			type(cubicSpline) :: this
227
228			if(associated(this%x)) then
229			deallocate(this%x)
230			this%x => null()
231			end if
232			if(associated(this%c)) then
233			deallocate(this%c)
234			this%c => null()
235			end if
236
237			this%np = 0
238
239			end subroutine deleteSpline
240
241	gezelter	2717	subroutine lookupNonuniformSpline(cs, xval, yval)
242	gezelter	2708
243			!*************************************************************************
244			!
245	gezelter	2717	! lookupNonuniformSpline evaluates a piecewise cubic Hermite interpolant.
246	gezelter	2708	!
247			! Discussion:
248			!
249			! newSpline must be called first, to set up the
250			! spline data from the raw function and derivative data.
251			!
252			! Modified:
253			!
254			! 06 April 1999
255			!
256			! Reference:
257			!
258			! Conte and de Boor,
259			! Algorithm PCUBIC,
260			! Elementary Numerical Analysis,
261			! 1973, page 234.
262			!
263			! Parameters:
264			!
265			implicit none
266
267			type (cubicSpline), intent(in) :: cs
268			real( kind = DP ), intent(in) :: xval
269			real( kind = DP ), intent(out) :: yval
270	gezelter	2709	real( kind = DP ) :: dx
271	gezelter	2708	integer :: i, j
272			!
273			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
274			! or is nearest to xval.
275			!
276			j = cs%np - 1
277
278			do i = 1, cs%np - 2
279
280			if ( xval < cs%x(i+1) ) then
281			j = i
282			exit
283			end if
284
285			end do
286			!
287			! Evaluate the cubic polynomial.
288			!
289			dx = xval - cs%x(j)
290
291			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
292
293			return
294	gezelter	2717	end subroutine lookupNonuniformSpline
295	gezelter	2708
296	gezelter	2717	subroutine lookupUniformSpline(cs, xval, yval)
297	gezelter	2708
298			!*************************************************************************
299			!
300	gezelter	2717	! lookupUniformSpline evaluates a piecewise cubic Hermite interpolant.
301	gezelter	2708	!
302			! Discussion:
303			!
304			! newSpline must be called first, to set up the
305			! spline data from the raw function and derivative data.
306			!
307			! Modified:
308			!
309			! 06 April 1999
310			!
311			! Reference:
312			!
313			! Conte and de Boor,
314			! Algorithm PCUBIC,
315			! Elementary Numerical Analysis,
316			! 1973, page 234.
317			!
318			! Parameters:
319			!
320			implicit none
321
322			type (cubicSpline), intent(in) :: cs
323			real( kind = DP ), intent(in) :: xval
324			real( kind = DP ), intent(out) :: yval
325	gezelter	2717	real( kind = DP ) :: a, b, c, d, dx
326	gezelter	2708	integer :: i, j
327			!
328			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
329			! or is nearest to xval.
330
331	gezelter	2709	j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
332	gezelter	2708
333			dx = xval - cs%x(j)
334
335	gezelter	2717	a = cs%c(1,j)
336			b = cs%c(2,j)
337			c = cs%c(3,j)
338			d = cs%c(4,j)
339
340			yval = c + dx * d
341			yval = b + dx * yval
342			yval = a + dx * yval
343	gezelter	2708
344			return
345	gezelter	2717	end subroutine lookupUniformSpline
346	gezelter	2711
347	gezelter	2717	subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
348
349			implicit none
350	gezelter	2711
351			type (cubicSpline), intent(in) :: cs
352	gezelter	2717	real( kind = DP ), intent(in) :: xval
353			real( kind = DP ), intent(out) :: yval, dydx
354			real( kind = DP ) :: a, b, c, d, dx
355			integer :: i, j
356
357			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
358			! or is nearest to xval.
359
360			j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
361
362			dx = xval - cs%x(j)
363
364			a = cs%c(1,j)
365			b = cs%c(2,j)
366			c = cs%c(3,j)
367			d = cs%c(4,j)
368
369			yval = c + dx * d
370			yval = b + dx * yval
371			yval = a + dx * yval
372
373			dydx = 2.0d0 * c + 3.0d0 * d * dx
374			dydx = b + dx * dydx
375
376			return
377			end subroutine lookupUniformSpline1d
378
379			subroutine lookupSpline(cs, xval, yval)
380
381			type (cubicSpline), intent(in) :: cs
382	gezelter	2711	real( kind = DP ), intent(inout) :: xval
383			real( kind = DP ), intent(inout) :: yval
384
385			if (cs%isUniform) then
386	gezelter	2717	call lookupUniformSpline(cs, xval, yval)
387	gezelter	2711	else
388	gezelter	2717	call lookupNonuniformSpline(cs, xval, yval)
389	gezelter	2711	endif
390
391			return
392	gezelter	2717	end subroutine lookupSpline
393	gezelter	2708
394	gezelter	2717	end module interpolation