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.1 2006-04-14 19:57:04 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(4,:) :: c => null()
  end type cubicSpline

  interface splineLookup
     module procedure multiSplint
     module procedure splintd
     module procedure splintd1
     module procedure splintd2
  end interface

  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,n) = 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,np) = 0.0_DP
    cs%c(4,np) = 0.0_DP

    cs%dx = dx
    cs%dxi = 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,n) = 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,np) = 0.0_DP
    cs%c(4,np) = 0.0_DP

    cs%dx = dx
    cs%dxi = 1.0_DP / dx

    return
  end subroutine newSplineWithoutDerivs

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