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 interplelation 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.2 2006-04-14 20:04:31 gezelter Exp $


module  INTERPOLATION
  use definitions
  use status
  implicit none
  PRIVATE

  character(len = statusMsgSize) :: errMSG

  type, public :: cubicSpline
     private
     integer :: np = 0
     real(kind=dp) :: dx
     real(kind=dp) :: dx_i
     real (kind=dp), pointer,dimension(:)   :: x => null()
     real (kind=dp), pointer,dimension(:,:) :: c => null()
  end type cubicSpline

  interface newSpline
     module procedure newSplineWithoutDerivs
     module procedure newSplineWithDerivs
  end interface

  public :: deleteSpline

contains


  subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary)

    !************************************************************************
    !
    ! newSplineWithoutDerivs 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 TAU are assumed to be
    !    strictly increasing.
    !
    !    Input, real y(I), contains the function value at x(I) for 
    !      I = 1, N.
    !
    !    yp1 contains the slope at x(1) and ypn contains
    !    the slope at x(N).
    !
    !    On output, the intermediate slopes at x(I) have been
    !    stored in cs%C(2,I), for I = 2 to N-1.

    implicit none

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

    alloc_error = 0

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

    ! make sure the sizes match

    if (size(x) .ne. size(y)) then
       call handleError("interpolation::newSplineWithoutDerivs", &
            "Array size mismatch")
    end if

    np = size(x)
    cs%np = np

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

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

    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
         (boundary.eq.'b').or.(boundary.eq.'B')) then
       cs%c(2,1) = yp1
    else
       cs%c(2,1) = 0.0_DP
    endif
    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
         (boundary.eq.'b').or.(boundary.eq.'B')) then
       cs%c(2,1) = ypn
    else
       cs%c(2,1) = 0.0_DP
    endif

    !
    !  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 = dx
    cs%dx_i = 1.0_DP / dx
    return
  end subroutine newSplineWithoutDerivs

  subroutine newSplineWithDerivs(cs, x, y, yp)

    !************************************************************************
    !
    ! newSplineWithDerivs 

    implicit none

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

    alloc_error = 0

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

    ! make sure the sizes match

    if ((size(x) .ne. size(y)).or.(size(x) .ne. size(yp))) then
       call handleError("interpolation::newSplineWithDerivs", &
            "Array size mismatch")
    end if
    
    np = size(x)
    cs%np = np

    allocate(cs%x(np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSplineWithDerivs", &
            "Error in allocating storage for x")
    endif
    
    allocate(cs%c(4,np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSplineWithDerivs", &
            "Error in allocating storage for c")
    endif
    
    do i = 1, np
       cs%x(i) = x(i)
       cs%c(1,i) = y(i)       
       cs%c(2,i) = yp(i)
    enddo
    !
    !  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 = dx
    cs%dx_i = 1.0_DP / dx

    return
  end subroutine newSplineWithDerivs

  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 lookup_nonuniform_spline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookup_nonuniform_spline 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 lookup_nonuniform_spline

  subroutine lookup_uniform_spline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookup_uniform_spline 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 = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))

    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 lookup_uniform_spline
  
end module INTERPOLATION
Revision:	932
Committed:	Fri Apr 14 20:04:31 2006 UTC (19 years ago) by gezelter
Original Path:	*trunk/src/utils/interpolation.F90*
File size:	12399 byte(s)
Log Message:	bug fixes for interpolation module
#	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			!! PURPOSE: Generic Spline interplelation routines. These routines assume that we are on a uniform grid for
47			!! precomputation of spline parameters.
48			!!
49			!! @author Charles F. Vardeman II
50	gezelter	932	!! @version $Id: interpolation.F90,v 1.2 2006-04-14 20:04:31 gezelter Exp $
51	gezelter	931
52
53			module INTERPOLATION
54			use definitions
55			use status
56			implicit none
57			PRIVATE
58
59			character(len = statusMsgSize) :: errMSG
60
61			type, public :: cubicSpline
62			private
63			integer :: np = 0
64			real(kind=dp) :: dx
65			real(kind=dp) :: dx_i
66			real (kind=dp), pointer,dimension(:) :: x => null()
67	gezelter	932	real (kind=dp), pointer,dimension(:,:) :: c => null()
68	gezelter	931	end type cubicSpline
69
70			interface newSpline
71			module procedure newSplineWithoutDerivs
72			module procedure newSplineWithDerivs
73			end interface
74
75			public :: deleteSpline
76
77			contains
78
79
80			subroutine newSplineWithoutDerivs(cs, x, y, yp1, ypn, boundary)
81
82			!************************************************************************
83			!
84			! newSplineWithoutDerivs solves for slopes defining a cubic spline.
85			!
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			! the data points. The entries of TAU are assumed to be
102			! strictly increasing.
103			!
104			! Input, real y(I), contains the function value at x(I) for
105			! I = 1, N.
106			!
107			! yp1 contains the slope at x(1) and ypn contains
108			! the slope at x(N).
109			!
110			! On output, the intermediate slopes at x(I) have been
111			! stored in cs%C(2,I), for I = 2 to N-1.
112
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			character(len=*), intent(in) :: boundary
119			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			call handleWarning("interpolation::newSplineWithoutDerivs", &
126			"Type was already created")
127			call deleteSpline(cs)
128			end if
129
130			! make sure the sizes match
131
132			if (size(x) .ne. size(y)) then
133			call handleError("interpolation::newSplineWithoutDerivs", &
134			"Array size mismatch")
135			end if
136
137			np = size(x)
138			cs%np = np
139
140			allocate(cs%x(np), stat=alloc_error)
141			if(alloc_error .ne. 0) then
142			call handleError("interpolation::newSplineWithoutDerivs", &
143			"Error in allocating storage for x")
144			endif
145
146			allocate(cs%c(4,np), stat=alloc_error)
147			if(alloc_error .ne. 0) then
148			call handleError("interpolation::newSplineWithoutDerivs", &
149			"Error in allocating storage for c")
150			endif
151
152			do i = 1, np
153			cs%x(i) = x(i)
154			cs%c(1,i) = y(i)
155			enddo
156
157			if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
158			(boundary.eq.'b').or.(boundary.eq.'B')) then
159			cs%c(2,1) = yp1
160			else
161			cs%c(2,1) = 0.0_DP
162			endif
163			if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
164			(boundary.eq.'b').or.(boundary.eq.'B')) then
165			cs%c(2,1) = ypn
166			else
167			cs%c(2,1) = 0.0_DP
168			endif
169
170			!
171			! Set up the right hand side of the linear system.
172			!
173			do i = 2, cs%np - 1
174			cs%c(2,i) = 3.0_DP * ( &
175			(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
176			(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
177			end do
178			!
179			! Set the diagonal coefficients.
180			!
181			cs%c(4,1) = 1.0_DP
182			do i = 2, cs%np - 1
183			cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
184			end do
185	gezelter	932	cs%c(4,cs%np) = 1.0_DP
186	gezelter	931	!
187			! Set the off-diagonal coefficients.
188			!
189			cs%c(3,1) = 0.0_DP
190			do i = 2, cs%np
191			cs%c(3,i) = x(i) - x(i-1)
192			end do
193			!
194			! Forward elimination.
195			!
196			do i = 2, cs%np - 1
197			g = -cs%c(3,i+1) / cs%c(4,i-1)
198			cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
199			cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
200			end do
201			!
202			! Back substitution for the interior slopes.
203			!
204			do i = cs%np - 1, 2, -1
205			cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
206			end do
207			!
208			! Now compute the quadratic and cubic coefficients used in the
209			! piecewise polynomial representation.
210			!
211			do i = 1, cs%np - 1
212			dx = x(i+1) - x(i)
213			divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
214			divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
215			cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
216			cs%c(4,i) = divdif3 / ( dx * dx )
217			end do
218
219	gezelter	932	cs%c(3,cs%np) = 0.0_DP
220			cs%c(4,cs%np) = 0.0_DP
221	gezelter	931
222			cs%dx = dx
223	gezelter	932	cs%dx_i = 1.0_DP / dx
224	gezelter	931	return
225			end subroutine newSplineWithoutDerivs
226
227			subroutine newSplineWithDerivs(cs, x, y, yp)
228
229			!************************************************************************
230			!
231			! newSplineWithDerivs
232
233			implicit none
234
235			type (cubicSpline), intent(inout) :: cs
236			real( kind = DP ), intent(in) :: x(:), y(:), yp(:)
237			real( kind = DP ) :: g, divdif1, divdif3, dx
238			integer :: i, alloc_error, np
239
240			alloc_error = 0
241
242			if (cs%np .ne. 0) then
243			call handleWarning("interpolation::newSplineWithDerivs", &
244			"Type was already created")
245			call deleteSpline(cs)
246			end if
247
248			! make sure the sizes match
249
250			if ((size(x) .ne. size(y)).or.(size(x) .ne. size(yp))) then
251			call handleError("interpolation::newSplineWithDerivs", &
252			"Array size mismatch")
253			end if
254
255			np = size(x)
256			cs%np = np
257
258			allocate(cs%x(np), stat=alloc_error)
259			if(alloc_error .ne. 0) then
260			call handleError("interpolation::newSplineWithDerivs", &
261			"Error in allocating storage for x")
262			endif
263
264			allocate(cs%c(4,np), stat=alloc_error)
265			if(alloc_error .ne. 0) then
266			call handleError("interpolation::newSplineWithDerivs", &
267			"Error in allocating storage for c")
268			endif
269
270			do i = 1, np
271			cs%x(i) = x(i)
272			cs%c(1,i) = y(i)
273			cs%c(2,i) = yp(i)
274			enddo
275			!
276			! Set the diagonal coefficients.
277			!
278			cs%c(4,1) = 1.0_DP
279			do i = 2, cs%np - 1
280			cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
281			end do
282	gezelter	932	cs%c(4,cs%np) = 1.0_DP
283	gezelter	931	!
284			! Set the off-diagonal coefficients.
285			!
286			cs%c(3,1) = 0.0_DP
287			do i = 2, cs%np
288			cs%c(3,i) = x(i) - x(i-1)
289			end do
290			!
291			! Forward elimination.
292			!
293			do i = 2, cs%np - 1
294			g = -cs%c(3,i+1) / cs%c(4,i-1)
295			cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
296			cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
297			end do
298			!
299			! Back substitution for the interior slopes.
300			!
301			do i = cs%np - 1, 2, -1
302			cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
303			end do
304			!
305			! Now compute the quadratic and cubic coefficients used in the
306			! piecewise polynomial representation.
307			!
308			do i = 1, cs%np - 1
309			dx = x(i+1) - x(i)
310			divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
311			divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
312			cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
313			cs%c(4,i) = divdif3 / ( dx * dx )
314			end do
315
316	gezelter	932	cs%c(3,cs%np) = 0.0_DP
317			cs%c(4,cs%np) = 0.0_DP
318	gezelter	931
319			cs%dx = dx
320	gezelter	932	cs%dx_i = 1.0_DP / dx
321	gezelter	931
322			return
323	gezelter	932	end subroutine newSplineWithDerivs
324	gezelter	931
325			subroutine deleteSpline(this)
326
327			type(cubicSpline) :: this
328
329			if(associated(this%x)) then
330			deallocate(this%x)
331			this%x => null()
332			end if
333			if(associated(this%c)) then
334			deallocate(this%c)
335			this%c => null()
336			end if
337
338			this%np = 0
339
340			end subroutine deleteSpline
341
342			subroutine lookup_nonuniform_spline(cs, xval, yval)
343
344			!*************************************************************************
345			!
346			! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
347			!
348			! Discussion:
349			!
350			! newSpline must be called first, to set up the
351			! spline data from the raw function and derivative data.
352			!
353			! Modified:
354			!
355			! 06 April 1999
356			!
357			! Reference:
358			!
359			! Conte and de Boor,
360			! Algorithm PCUBIC,
361			! Elementary Numerical Analysis,
362			! 1973, page 234.
363			!
364			! Parameters:
365			!
366			implicit none
367
368			type (cubicSpline), intent(in) :: cs
369			real( kind = DP ), intent(in) :: xval
370			real( kind = DP ), intent(out) :: yval
371	gezelter	932	real( kind = DP ) :: dx
372	gezelter	931	integer :: i, j
373			!
374			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
375			! or is nearest to xval.
376			!
377			j = cs%np - 1
378
379			do i = 1, cs%np - 2
380
381			if ( xval < cs%x(i+1) ) then
382			j = i
383			exit
384			end if
385
386			end do
387			!
388			! Evaluate the cubic polynomial.
389			!
390			dx = xval - cs%x(j)
391
392			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
393
394			return
395			end subroutine lookup_nonuniform_spline
396
397			subroutine lookup_uniform_spline(cs, xval, yval)
398
399			!*************************************************************************
400			!
401			! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
402			!
403			! Discussion:
404			!
405			! newSpline must be called first, to set up the
406			! spline data from the raw function and derivative data.
407			!
408			! Modified:
409			!
410			! 06 April 1999
411			!
412			! Reference:
413			!
414			! Conte and de Boor,
415			! Algorithm PCUBIC,
416			! Elementary Numerical Analysis,
417			! 1973, page 234.
418			!
419			! Parameters:
420			!
421			implicit none
422
423			type (cubicSpline), intent(in) :: cs
424			real( kind = DP ), intent(in) :: xval
425			real( kind = DP ), intent(out) :: yval
426	gezelter	932	real( kind = DP ) :: dx
427	gezelter	931	integer :: i, j
428			!
429			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
430			! or is nearest to xval.
431
432	gezelter	932	j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
433	gezelter	931
434			dx = xval - cs%x(j)
435
436			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
437
438			return
439			end subroutine lookup_uniform_spline
440
441			end module INTERPOLATION