!>
!> @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/>.
!>
!> @author
!> (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