Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/cspint.f90
03540012fedefeature/refactoring

Last Change
Raw File

cspint.f90
View Options

!>
!> @file cspint.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine cspint ( ftab, xtab, ntab, a, b, y, e, work, result, ind )
!*************************************************************************
!
!! CSPINT estimates an integral using a spline interpolant.
!
! Discussion:
!
! CSPINT estimates the integral from A to B of F(X) by
! computing the natural spline S(X) that interpolates to F
! and integrating that exactly.
!
! F is supplied to the routine in the form of tabulated data.
!
! Other output from the program includes the definite integral
! from X(1) to X(I) of the spline, and the coefficients
! necessary for the user to evaluate the spline outside of
! this routine.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) FTAB(NTAB), contains the tabulated values
! of the functions, FTAB(I)=F(XTAB(I)).
!
! Input, real ( kind = 8 ) XTAB(NTAB), contains the points at
! which the function was evaluated. The XTAB's must be
! distinct and in ascending order.
!
! Input, integer NTAB, the number of entries in FTAB
! and XTAB. NTAB must be at least 3.
!
! Input, real ( kind = 8 ) A, lower limit of integration.
!
! Input, real ( kind = 8 ) B, upper limit of integration.
!
! Output, real ( kind = 8 ) Y(3,NTAB), will contain the coefficients
! of the interpolating natural spline over each subinterval.
!
! For XTAB(I) < = X <= XTAB(I+1),
!
! S(X) = FTAB(I) + Y(1,I)*(X-XTAB(I)) + Y(2,I)*(X-XTAB(I))**2
! + Y(3,I)*(X-XTAB(I))**3
!
! Output, real ( kind = 8 ) E(NTAB), E(I)=the definite integral
! from XTAB(1) to XTAB(I) of S(X).
!
! Workspace, real ( kind = 8 ) WORK(NTAB).
!
! Output, real ( kind = 8 ) RESULT, the estimated value of the integral.
!
! Output, integer IND, error flag.
! IND=0 if NTAB < 3 or the XTAB's are not distinct and in
! ascending order.
! IND=1 otherwise.
!
implicit none
integer ntab
real ( kind = 8 ) a
real ( kind = 8 ) b
real ( kind = 8 ) e(ntab)
real ( kind = 8 ) ftab(ntab)
integer i
integer ind
integer j
real ( kind = 8 ) r
real ( kind = 8 ) result
real ( kind = 8 ) s
real ( kind = 8 ) term
real ( kind = 8 ) u
real ( kind = 8 ) work(ntab)
real ( kind = 8 ) xtab(ntab)
real ( kind = 8 ) y(3,ntab)
ind = 0
if ( ntab < 3 ) then
write(*,*)' '
write(*,*)'CSPINT - Fatal error!'
write(*,*)' NTAB must be at least 3,'
write(*,*)' but your value was NTAB = ',ntab
stop
end if
do i = 1, ntab-1
if ( xtab(i+1) <= xtab(i) ) then
write(*,*)' '
write(*,*)'CSPINT - Fatal error!'
write(*,*)' Interval ',i,' is illegal.'
write(*,*)' XTAB(I) =',xtab(i)
write(*,*)' XTAB(I+1)=',xtab(i+1)
stop
end if
end do
s = 0.0D+00
do i = 1, ntab-1
r = ( ftab(i+1) - ftab(i) ) / ( xtab(i+1) - xtab(i) )
y(2,i) = r - s
s = r
end do
result = 0.0D+00
s = 0.0D+00
r = 0.0D+00
y(2,1) = 0.0D+00
y(2,ntab) = 0.0D+00
do i = 2, ntab-1
y(2,i) = y(2,i) + r * y(2,i-1)
work(i) = 2.0D+00 * ( xtab(i-1) - xtab(i+1) ) - r * s
s = xtab(i+1) - xtab(i)
r = s / work(i)
end do
do j = 2, ntab-1
i = ntab+1-j
y(2,i) = ((xtab(i+1)-xtab(i))*y(2,i+1)-y(2,i))/work(i)
end do
do i = 1, ntab-1
s = xtab(i+1) - xtab(i)
r = y(2,i+1) - y(2,i)
y(3,i) = r / s
y(2,i) = 3.0D+00 * y(2,i)
y(1,i) = ( ftab(i+1) - ftab(i) ) / s - ( y(2,i) + r ) * s
end do
e(1) = 0.0D+00
do i = 1, ntab-1
s = xtab(i+1) - xtab(i)
term = ( ( ( y(3,i) * 0.25D+00 * s &
+ y(2,i) / 3.0D+00 ) * s &
+ y(1,i) * 0.5D+00 ) * s + ftab(i) ) * s
e(i+1) = e(i) + term
end do
!
! Determine where the endpoints A and B lie in the mesh of XTAB's.
!
r = a
u = 1.0D+00
do j = 1, 2
if ( r <= xtab(1) ) then
result = result-u*((r-xtab(1))*y(1,1)* 0.5D+00 + ftab(1))*(r-xtab(1))
else if ( xtab(ntab) <= r ) then
result = result-u*(e(ntab)+(r-xtab(ntab))*(ftab(ntab) + 0.5D+00 * &
(ftab(ntab-1)+(xtab(ntab)-xtab(ntab-1))*y(1,ntab-1))*(r- &
xtab(ntab))))
else
do i = 1, ntab-1
if ( r <= xtab(i+1) ) then
r = r - xtab(i)
result = result-u*(e(i)+(((y(3,i)*0.25D+00*r+y(2,i)/3.0D+00)*r &
+y(1,i) * 0.5D+00 )*r+ftab(i))*r)
go to 100
end if
end do
end if
100 continue
u = -1.0D+00
r = b
end do
ind = 1
return
end
c4science · Help