Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/bvalue.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

bvalue.f90
View Options

!>
!> @file bvalue.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>
!>
function bvalue ( t, bcoef, n, k, x, jderiv )
!*************************************************************************
!
!! BVALUE evaluates a derivative of a spline from its B-spline representation.
!
! Discussion:
!
! The spline is taken to be continuous from the right.
!
! The nontrivial knot interval (T(I),T(I+1)) containing X is
! located with the aid of INTERV. The K B-spline coefficients
! of F relevant for this interval are then obtained from BCOEF,
! or are taken to be zero if not explicitly available, and are
! then differenced JDERIV times to obtain the B-spline
! coefficients of (D**JDERIV)F relevant for that interval.
!
! Precisely, with J = JDERIV, we have from X.(12) of the text that:
!
! (D**J)F = sum ( BCOEF(.,J)*B(.,K-J,T) )
!
! where
! / BCOEF(.), , J == 0
! /
! BCOEF(.,J) = / BCOEF(.,J-1) - BCOEF(.-1,J-1)
! / -----------------------------, 0 < J
! / (T(.+K-J) - T(.))/(K-J)
!
! Then, we use repeatedly the fact that
!
! sum ( A(.)*B(.,M,T)(X) ) = sum ( A(.,X)*B(.,M-1,T)(X) )
!
! with
! (X - T(.))*A(.) + (T(.+M-1) - X)*A(.-1)
! A(.,X) = ---------------------------------------
! (X - T(.)) + (T(.+M-1) - X)
!
! to write (D**J)F(X) eventually as a linear combination of
! B-splines of order 1, and the coefficient for B(I,1,T)(X)
! must then be the desired number (D**J)F(X).
! See x.(17)-(19) of text.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) T(N+K), the knot sequence. T is assumed
! to be nondecreasing.
!
! Input, real ( kind = 8 ) BCOEF(N), B-spline coefficient sequence.
!
! Input, integer N, the length of BCOEF.
!
! Input, integer K, the order of the spline.
!
! Input, real ( kind = 8 ) X, the point at which to evaluate.
!
! Input, integer JDERIV, the order of the derivative to
! be evaluated. JDERIV is assumed to be zero or positive.
!
! Output, real ( kind = 8 ) BVALUE, the value of the (JDERIV)-th
! derivative of the spline at X.
!
implicit none
integer k
integer n
real ( kind = 8 ) aj(k)
real ( kind = 8 ) bcoef(n)
real ( kind = 8 ) bvalue
real ( kind = 8 ) dl(k)
real ( kind = 8 ) dr(k)
integer i
integer ilo
integer j
integer jc
integer jcmax
integer jcmin
integer jderiv
integer jj
integer mflag
real ( kind = 8 ) t(n+k)
real ( kind = 8 ) x
bvalue = 0.0D+00
if ( k <= jderiv ) then
return
end if
!
! Find I so that 1 <= i < n+k and t(i) < t(i+1) and t(i) <= x < t(i+1).
!
! If no such i can be found, X lies
! outside the support of the spline F and bvalue=0.
! (the asymmetry in this choice of i makes F rightcontinuous)
!
call interv ( t, n+k, x, i, mflag )
if ( mflag /= 0 ) then
return
end if
!
! If K=1 (and jderiv = 0), bvalue = bcoef(i).
!
if ( k <= 1 ) then
bvalue = bcoef(i)
return
end if
!
! Store the K b-spline coefficients relevant for the knot interval
! (T(i),T(i+1)) in aj(1),...,aj(k) and compute dl(j)=x-t(i+1-j),
! dr(j)=T(i+j)-x, j=1,...,k-1 . set any of the aj not obtainable
! from input to zero. Set any T's not obtainable equal to T(1) or
! to T(n+k) appropriately.
!
jcmin = 1
if ( k <= i ) then
do j = 1, k-1
dl(j) = x-t(i+1-j)
end do
else
jcmin = 1-(i-k)
do j = 1, i
dl(j) = x-t(i+1-j)
end do
do j = i, k-1
aj(k-j) = 0.0D+00
dl(j) = dl(i)
end do
end if
jcmax = k
if ( i <= n ) then
go to 90
end if
jcmax = k + n - i
do j = 1, k+n-i
dr(j) = t(i+j)-x
end do
do j = k+n-i, k-1
aj(j+1) = 0.0D+00
dr(j) = dr(k+n-i)
end do
go to 110
90 continue
do j = 1, k-1
dr(j) = t(i+j)-x
end do
110 continue
do jc = jcmin, jcmax
aj(jc) = bcoef(i-k+jc)
end do
!
! Difference the coefficients JDERIV times.
!
do j = 1, jderiv
ilo = k-j
do jj = 1, k-j
aj(jj) = ((aj(jj+1)-aj(jj))/(dl(ilo)+dr(jj))) * real ( k - j, kind = 8 )
ilo = ilo-1
end do
end do
!
! Compute value at X in (t(i),t(i+1)) of jderiv-th derivative,
! given its relevant b-spline coefficients in aj(1),...,aj(k-jderiv).
!
do j = jderiv+1, k-1
ilo = k-j
do jj = 1, k-j
aj(jj) = ( aj(jj+1) * dl(ilo) + aj(jj) * dr(jj) ) &
/ ( dl(ilo) + dr(jj) )
ilo = ilo-1
end do
end do
bvalue = aj(1)
return
end
c4science · Help