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:	2722
Committed:	Thu Apr 20 18:24:24 2006 UTC (18 years, 2 months 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	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	2722	!! @version $Id: interpolation.F90,v 1.7 2006-04-20 18:24:24 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			cs%dx_i = 1.0d0 / (x(2) - x(1))
228			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			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	2722	this%n = 0
265	gezelter	2708
266			end subroutine deleteSpline
267
268	gezelter	2717	subroutine lookupNonuniformSpline(cs, xval, yval)
269	gezelter	2708
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	2722	real( kind = DP ) :: dx
276	gezelter	2708	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	2722	j = cs%n - 1
282	gezelter	2708
283	gezelter	2722	do i = 0, cs%n - 2
284	gezelter	2708
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	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
296	gezelter	2708
297			return
298	gezelter	2717	end subroutine lookupNonuniformSpline
299	gezelter	2708
300	gezelter	2717	subroutine lookupUniformSpline(cs, xval, yval)
301	gezelter	2708
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	2722	real( kind = DP ) :: dx
308	gezelter	2708	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	2722
313			j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
314
315	gezelter	2708	dx = xval - cs%x(j)
316	gezelter	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
317	gezelter	2708
318			return
319	gezelter	2717	end subroutine lookupUniformSpline
320	gezelter	2711
321	gezelter	2717	subroutine lookupUniformSpline1d(cs, xval, yval, dydx)
322
323			implicit none
324	gezelter	2711
325			type (cubicSpline), intent(in) :: cs
326	gezelter	2717	real( kind = DP ), intent(in) :: xval
327			real( kind = DP ), intent(out) :: yval, dydx
328	gezelter	2722	real( kind = DP ) :: dx
329	gezelter	2717	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	2722	j = MAX(1, MIN(cs%n-1, idint((xval-cs%x(1)) * cs%dx_i) + 1))
336
337	gezelter	2717	dx = xval - cs%x(j)
338	gezelter	2722	yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
339	gezelter	2717
340	gezelter	2722	dydx = cs%b(j) + dx(2.0d0 cs%c(j) + 3.0d0 * dx * cs%d(j))
341	gezelter	2717
342			return
343			end subroutine lookupUniformSpline1d
344
345			subroutine lookupSpline(cs, xval, yval)
346
347			type (cubicSpline), intent(in) :: cs
348	gezelter	2711	real( kind = DP ), intent(inout) :: xval
349			real( kind = DP ), intent(inout) :: yval
350
351			if (cs%isUniform) then
352	gezelter	2717	call lookupUniformSpline(cs, xval, yval)
353	gezelter	2711	else
354	gezelter	2717	call lookupNonuniformSpline(cs, xval, yval)
355	gezelter	2711	endif
356
357			return
358	gezelter	2717	end subroutine lookupSpline
359	gezelter	2708
360	gezelter	2717	end module interpolation