Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/cubslo.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

cubslo.f90
View Options

!>
!> @file cubslo.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/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine cubslo ( tau, c, n )
!*******************************************************************************
!
!! CUBSLO solves for slopes defining a cubic spline.
!
! Discussion:
!
! A tridiagonal linear system for the unknown slopes S(I) of
! F at TAU(I), I=1,..., N, is generated and then solved by Gauss
! elimination, with S(I) ending up in C(2,I), for all I.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) TAU(N), the abscissas or X values of
! the data points. The entries of TAU are assumed to be
! strictly increasing.
!
! Input, integer N, the number of data points. N is
! assumed to be at least 2.
!
! Input/output, real ( kind = 8 ) C(4,N).
! On input, C(1,I) contains the function value at TAU(I),
! for I = 1 to N.
! C(2,1) contains the slope at TAU(1) and C(2,N) contains
! the slope at TAU(N).
! On output, the intermediate slopes at TAU(I) have been
! stored in C(2,I), for I = 2 to N-1.
!
implicit none
integer n
real ( kind = 8 ) c(4,n)
real ( kind = 8 ) g
integer i
integer ibcbeg
integer ibcend
real ( kind = 8 ) tau(n)
!
! Set up the right hand side of the linear system.
! C(2,1) and C(2,N) are presumably already set.
!
do i = 2, n-1
c(2,i) = 3.0D+00 * ( &
( tau(i) - tau(i-1) ) * ( c(1,i+1) - c(1,i) ) / ( tau(i+1) - tau(i) ) + &
( tau(i+1) - tau(i) ) * ( c(1,i) - c(1,i-1) ) / ( tau(i) - tau(i-1) ) )
end do
!
! Set the diagonal coefficients.
!
c(4,1) = 1.0D+00
do i = 2, n-1
c(4,i) = 2.0D+00 * ( tau(i+1) - tau(i-1) )
end do
c(4,n) = 1.0D+00
!
! Set the off-diagonal coefficients.
!
c(3,1) = 0.0D+00
do i = 2, n
c(3,i) = tau(i) - tau(i-1)
end do
!
! Forward elimination.
!
do i = 2, n-1
g = -c(3,i+1) / c(4,i-1)
c(4,i) = c(4,i) + g * c(3,i-1)
c(2,i) = c(2,i) + g * c(2,i-1)
end do
!
! Back substitution for the interior slopes.
!
do i = n-1, 2, -1
c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i)
end do
return
end
c4science · Help