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:	2709
Committed:	Fri Apr 14 20:04:31 2006 UTC (18 years, 2 months ago) by gezelter
File size:	12399 byte(s)
Log Message:	bug fixes for interpolation module
#	Content
1	!!
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.2 2006-04-14 20:04:31 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(:,:) :: c => null()
68	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	cs%c(4,cs%np) = 1.0_DP
186	!
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	cs%c(3,cs%np) = 0.0_DP
220	cs%c(4,cs%np) = 0.0_DP
221
222	cs%dx = dx
223	cs%dx_i = 1.0_DP / dx
224	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	cs%c(4,cs%np) = 1.0_DP
283	!
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	cs%c(3,cs%np) = 0.0_DP
317	cs%c(4,cs%np) = 0.0_DP
318
319	cs%dx = dx
320	cs%dx_i = 1.0_DP / dx
321
322	return
323	end subroutine newSplineWithDerivs
324
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	real( kind = DP ) :: dx
372	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	real( kind = DP ) :: dx
427	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	j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
433
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