Homec4science
Diffusion SPClibs (master) Browse

spclibs/pppack/setupq.f90
93a1e9712d04master

Last Change
Raw File

setupq.f90
View Options

!>
!> @file setupq.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 setupq ( x, dx, y, npoint, v, qty )
!*************************************************************************
!
!! SETUPQ is to be called in SMOOTH.
!
! Discussion:
!
! put delx=x(.+1)-x(.) into v(.,4),
! put the three bands of q-transp*d into v(.,1-3), and
! put the three bands of (d*q)-transp*(d*q) at and above the diagonal
! into v(.,5-7) .
!
! here, q is the tridiagonal matrix of order (npoint-2,npoint)
! with general row 1/delx(i) , -1/delx(i)-1/delx(i+1) , 1/delx(i+1)
! and d is the diagonal matrix with general row dx(i) .
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
implicit none
integer npoint
real ( kind = 8 ) diff
real ( kind = 8 ) dx(npoint)
integer i
real ( kind = 8 ) prev
real ( kind = 8 ) qty(npoint)
real ( kind = 8 ) v(npoint,7)
real ( kind = 8 ) x(npoint)
real ( kind = 8 ) y(npoint)
v(1,4) = x(2)-x(1)
do i = 2, npoint-1
v(i,4) = x(i+1) - x(i)
v(i,1) = dx(i-1) / v(i-1,4)
v(i,2) = -dx(i) / v(i,4) - dx(i) / v(i-1,4)
v(i,3) = dx(i+1) / v(i,4)
end do
v(npoint,1) = 0.0D+00
do i = 2, npoint-1
v(i,5) = v(i,1)**2 + v(i,2)**2 + v(i,3)**2
end do
do i = 3, npoint-1
v(i-1,6) = v(i-1,2)*v(i,1)+v(i-1,3)*v(i,2)
end do
v(npoint-1,6) = 0.0D+00
do i = 4, npoint-1
v(i-2,7) = v(i-2,3) * v(i,1)
end do
v(npoint-2,7) = 0.0D+00
v(npoint-1,7) = 0.0D+00
!
! Construct q-transp. * y in QTY.
!
prev = (y(2)-y(1)) / v(1,4)
do i = 2, npoint-1
diff = (y(i+1)-y(i)) / v(i,4)
qty(i) = diff-prev
prev = diff
end do
return
end
c4science · Help