!>
!> @file tautsp.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 tautsp ( tau, gtau, ntau, gamma, s, break, coef, l, k, iflag )

!*************************************************************************
!
!! TAUTSP constructs a cubic spline interpolant to given data.
!
!  Discussion:
!
!    If 0 < GAMMA, additional knots are introduced where needed to
!    make the interpolant more flexible locally.  This avoids extraneous
!    inflection points typical of cubic spline interpolation at knots to
!    rapidly changing data.
!
!  Method:
!
!    On the I-th interval, (TAU(I), TAU(I+1)), the interpolant is of the
!    form:
!
!    (*)  f(u(x))=a+b*u + c*h(u,z) + d*h(1-u,1-z) ,
!
!    with
!
!      U = U(X) = ( X - TAU(I) ) / DTAU(I).
!
!    Here,
!       z=z(i) = addg(i+1)/(addg(i)+addg(i+1))
!  (= .5, in case the denominator vanishes). with
!       addg(j)=abs(ddg(j)), ddg(j) = dg(j+1)-dg(j),
!       dg(j)=divdif ( j) = (gtau(j+1)-gtau(j))/dtau(j)
!  and
!       h(u,z)=alpha*u**3+(1-alpha)*(max(((u-zeta)/(1-zeta)),0)**3
!  with
!       alpha(z)=(1-gamma/3)/zeta
!       zeta(z)=1-gamma*min((1 - z), 1/3)
!  thus, for 1/3 <= z .le. 2/3,  f  is just a cubic polynomial on
!  the interval i. otherwise, it has one additional knot, at
!         tau(i)+zeta*dtau(i) .
!  as  z  approaches  1, h(.,z) has an increasingly sharp bend  near 1,
!  thus allowing  f  to turn rapidly near the additional knot.
!     in terms of f(j)=gtau(j) and
!       fsecnd(j)= second derivative of  f  at  tau(j),
!  the coefficients for (*) are given as
!       a=f(i)-d
!       b=(f(i+1)-f(i)) - (c - d)
!       c=fsecnd(i+1)*dtau(i)**2/hsecnd(1,z)
!       d=fsecnd(i)*dtau(i)**2/hsecnd(1,1-z)
!  hence can be computed once fsecnd(i),i=1,...,ntau, is fixed.
!
!    F is automatically continuous and has a continuous second derivative
!    (except when z=0 or 1 for some i). we determine fscnd(.) from
!    the requirement that also the first derivative of F be continuous.
!
!    In addition, we require that the third derivative be continuous
!    across TAU(2) and across TAU(NTAU-1).  This leads to a strictly
!    diagonally dominant tridiagonal linear system for the fsecnd(i)
!    which we solve by Gauss elimination without pivoting.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) TAU(NTAU), the sequence of data points.
!    TAU must be strictly increasing.
!
!    Input, real ( kind = 8 ) GTAU(NTAU), the corresponding sequence of
!    function values.
!
!    Input, integer NTAU, the number of data points.  NTAU must be at least 4.
!
!    Input, gamma  indicates whether additional flexibility is desired.
!        =0., no additional knots
!         in (0.,3.), under certain conditions on the given data at
!                points i-1, i, i+1, and i+2, a knot is added in the
!                i-th interval, i=2,...,ntau-2. see description of meth-
!                od below. the interpolant gets rounded with increasing
!                gamma. a value of  2.5  for gamma is typical.
!          in (3.,6.), same , except that knots might also be added in
!                intervals in which an inflection point would be permit-
!                ted.  a value of  5.5  for gamma is typical.
!
!    Output, break, coef, l, k  give the pp-representation of the interpolant.
!          specifically, for break(i) <= x .le. break(i+1), the
!        interpolant has the form
!  f(x)=coef(1,i) +dx(coef(2,i) +(dx/2)(coef(3,i) +(dx/3)coef(4,i)))
!        with  dx=x-break(i) and i=1,...,l .
!
!    Output, iflag =1, ok
!        =2, input was incorrect. a printout specifying the mistake
!            was made.
!            workspace
!
!    Output, s     is required, of size (ntau,6). the individual columns of this
!        array contain the following quantities mentioned in the write-
!        up and below.
!     s(.,1)=dtau = tau(.+1)-tau
!     s(.,2)=diag = diagonal in linear system
!     s(.,3)=u = upper diagonal in linear system
!     s(.,4)=r = right side for linear system (initially)
!          = fsecnd = solution of linear system , namely the second
!                       derivatives of interpolant at  tau
!     s(.,5)=z = indicator of additional knots
!     s(.,6)=1/hsecnd(1,x) with x = z or = 1-z. see below.
!
  implicit none

  integer ntau

  real ( kind = 8 ) alph
  real ( kind = 8 ) alpha
  real ( kind = 8 ) break(*)
  real ( kind = 8 ) c
  real ( kind = 8 ) coef(4,*)
  real ( kind = 8 ) d
  real ( kind = 8 ) del
  real ( kind = 8 ) denom
  real ( kind = 8 ) divdif
  real ( kind = 8 ) entry
  real ( kind = 8 ) entry3
  real ( kind = 8 ) factor
  real ( kind = 8 ) factr2
  real ( kind = 8 ) gam
  real ( kind = 8 ) gamma
  real ( kind = 8 ) gtau(ntau)
  integer i
  integer iflag
  integer k
  integer l
  integer method
  real ( kind = 8 ) onemg3
  real ( kind = 8 ) onemzt
  real ( kind = 8 ) ratio
  real ( kind = 8 ) s(ntau,6)
  real ( kind = 8 ) sixth
  real ( kind = 8 ) tau(ntau)
  real ( kind = 8 ) temp
  real ( kind = 8 ) x
  real ( kind = 8 ) z
  real ( kind = 8 ) zeta
  real ( kind = 8 ) zt2

  alph(x) = min ( 1.0D+00, onemg3 / x )
!
!  There must be at least 4 interpolation points.
!
  if ( ntau < 4 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TAUTSP - Fatal error!'
    write ( *, '(a)' ) '  Input NTAU must be at least 4.'
    write ( *, '(a,i6)' ) '  NTAU = ', ntau
    iflag = 2
    stop
  end if
!
!  Construct delta tau and first and second (divided) differences of data.
!
  do i = 1, ntau-1

    s(i,1) = tau(i+1)-tau(i)

    if ( s(i,1) <= 0.0D+00 ) then
      write(*,30)i,tau(i),tau(i+1)
   30     format (' point ',i3,' and the next',2e15.6,' are disordered')
      iflag=2
      return
    end if

    s(i+1,4) = ( gtau(i+1) - gtau(i) ) / s(i,1)
  end do

  do i = 2, ntau-1
    s(i,4) = s(i+1,4)-s(i,4)
  end do
!
!  Construct system of equations for second derivatives at  tau. at each
!  interior data point, there is one continuity equation, at the first
!  and the last interior data point there is an additional one for a
!  total of NTAU equations in  ntau  unknowns.
!
  i = 2
  s(2,2) = s(1,1) / 3.0D+00
  sixth = 1.0D+00 / 6.0D+00
  method = 2
  gam = gamma

  if ( gam <= 0.0D+00 ) then
    method=1
  end if

  if ( 3.0D+00 < gam ) then
    method = 3
    gam = gam - 3.0D+00
  end if

  onemg3 = 1.0D+00 - gam / 3.0D+00
!
!  loop over i
!
   70 continue
!
!  Construct z(i) and zeta(i)
!
  z = 0.5D+00

  if ( method == 1) then
    go to 100
  end if

  if ( method == 3) then
    go to 90
  end if

  if ( s(i,4)*s(i+1,4) < 0.0D+00 ) then
    go to 100
  end if

   90 continue

  temp = abs ( s(i+1,4) )
  denom = abs ( s(i,4) ) +temp

  if ( denom /= 0.0D+00 ) then
    z = temp/denom
    if ( abs ( z - 0.5D+00 ) <= sixth ) then
      z=0.5D+00
    end if
  end if

  100 continue

  s(i,5) = z
!
!  Set up part of the i-th equation which depends on the i-th interval.
!
  if ( z < 0.5D+00 ) then

    zeta = gam*z
    onemzt = 1.0D+00 - zeta
    zt2 = zeta**2
    alpha = alph(onemzt)
    factor = zeta/(alpha*(zt2-1.0D+00 ) + 1.0D+00 )
    s(i,6) = zeta*factor / 6.0D+00
    s(i,2) = s(i,2) + s(i,1) &
      * ( ( 1.0D+00 - alpha * onemzt ) * factor / 2.0D+00-s(i,6))
!
!  If z=0 and the previous z = 1, then d(i) = 0. since then
!  also u(i-1)=l(i+1) = 0, its value does not matter. reset
!  d(i)=1 to insure nonzero pivot in elimination.
!
    if ( s(i,2) <= 0.0D+00 ) then
      s(i,2) = 1.0D+00
    end if

    s(i,3)=s(i,1) / 6.0D+00
  else if ( z - 0.5D+00 == 0.0D+00 ) then
    s(i,2)=s(i,2)+s(i,1) / 3.0D+00
    s(i,3)=s(i,1) / 6.0D+00
  else if ( 0.0D+00 < z - 0.5D+00 ) then
    onemzt = gam*(1.0D+00 - z)
    zeta = 1.0D+00 - onemzt
    alpha = alph(zeta)
    factor = onemzt/(1.0D+00 - alpha * zeta * ( 1.0D+00 + onemzt ) )
    s(i,6) = onemzt*factor / 6.0D+00
    s(i,2) = s(i,2)+s(i,1) / 3.0D+00
    s(i,3) = s(i,6)*s(i,1)
  end if

  if ( 2 < i ) then
    go to 190
  end if

  s(1,5) = 0.5D+00
!
!  The first two equations enforce continuity of the first and of
!  the third derivative across tau(2).
!
  s(1,2)=s(1,1) / 6.0D+00
  s(1,3)=s(2,2)
  entry3=s(2,3)
  if ( z-0.5D+00) 150, 160, 170

  150 continue

  factr2 = zeta * ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 ) &
    / ( alpha * ( zeta * zt2 - 1.0D+00 ) + 1.0D+00 )

  ratio=factr2*s(2,1)/s(1,2)
  s(2,2)=factr2*s(2,1)+s(1,1)
  s(2,3)=-factr2*s(1,1)
  go to 180

  160 continue

  ratio=s(2,1)/s(1,2)
  s(2,2)=s(2,1)+s(1,1)
  s(2,3)=-s(1,1)
  go to 180

  170 continue

  ratio=s(2,1)/s(1,2)
  s(2,2)=s(2,1)+s(1,1)
  s(2,3)=-s(1,1)*6.0D+00 * alpha * s(2,6)
!
!  At this point, the first two equations read
!              diag(1)*x1+u(1)*x2 + entry3*x3=r(2)
!       -ratio*diag(1)*x1+diag(2)*x2 + u(2)*x3=0.0
!  Eliminate first unknown from second equation
!
  180 continue

  s(2,2)=ratio*s(1,3)+s(2,2)
  s(2,3)=ratio*entry3+s(2,3)
  s(1,4)=s(2,4)
  s(2,4)=ratio*s(1,4)
  go to 200

  190 continue
!
!  The i-th equation enforces continuity of the first derivative
!  across tau(i). it has been set up in statements 35 up to 40
!  and 21 up to 25 and reads now
!    -ratio*diag(i-1)*xi-1+diag(i)*xi + u(i)*xi+1=r(i) .
!  eliminate (i-1)st unknown from this equation
!
  s(i,2)=ratio*s(i-1,3)+s(i,2)
  s(i,4)=ratio*s(i-1,4)+s(i,4)
!
!  Set up the part of the next equation which depends on the
!  i-th interval.
!
  200 continue

  if ( z- 0.5D+00 ) 210, 220, 230

  210 continue
  ratio = -s(i,6) * s(i,1) / s(i,2)
  s(i+1,2)=s(i,1) / 3.0D+00
  go to 240

  220 continue
  ratio=-(s(i,1) / 6.0D+00 ) / s(i,2)
  s(i+1,2)=s(i,1) / 3.0D+00
  go to 240

  230 continue
  ratio=-( s(i,1) / 6.0D+00 )/s(i,2)
  s(i+1,2)=s(i,1)*((1.0D+00-zeta*alpha) * factor / 2.0D+00 - s(i,6) )
!
!  end of i loop
!
  240 continue

  i=i+1
  if ( i < ntau-1) then
    go to 70
  end if

  s(i,5) = 0.5D+00
!
!  The last two equations enforce continuity of third derivative and
!  of first derivative across  tau(ntau-1).
!
  entry=ratio*s(i-1,3)+s(i,2)+s(i,1)/3.0D+00
  s(i+1,2)=s(i,1)/6.0D+00
  s(i+1,4)=ratio*s(i-1,4)+s(i,4)
  if ( z- 0.5D+00 ) 250, 260, 270

  250 continue

  ratio = s(i,1) * 6.0D+00 * s(i-1,6) * alpha / s(i-1,2)
  s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
  s(i,3)=-s(i-1,1)
  go to 280

  260 continue

  ratio=s(i,1)/s(i-1,2)
  s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
  s(i,3)=-s(i-1,1)
  go to 280

  270 continue

  factr2=onemzt*(alpha*(onemzt**2-1.0D+00)+1.0D+00) &
    /(alpha*(onemzt**3-1.0D+00)+1.0D+00)

  ratio = factr2*s(i,1) / s(i-1,2)
  s(i,2)=ratio*s(i-1,3)+factr2*s(i-1,1)+s(i,1)
  s(i,3)=-factr2*s(i-1,1)
!
!  At this point, the last two equations read:
!
!           diag(i)*xi+     u(i)*xi+1=r(i)
!    -ratio*diag(i)*xi+diag(i+1)*xi+1=r(i+1)
!
!  Eliminate XI from the last equation.
!
  280 continue

  s(i,4)=ratio*s(i-1,4)
  ratio=-entry/s(i,2)
  s(i+1,2)=ratio*s(i,3)+s(i+1,2)
  s(i+1,4)=ratio*s(i,4)+s(i+1,4)
!
!  Back substitution.
!
  s(ntau,4) = s(ntau,4) / s(ntau,2)

  290 continue

  s(i,4)=(s(i,4)-s(i,3)*s(i+1,4))/s(i,2)
  i=i-1
  if ( 1 < i ) then
    go to 290
  end if

  s(1,4)=(s(1,4)-s(1,3)*s(2,4)-entry3*s(3,4))/s(1,2)
!
!  Construct polynomial pieces.
!
  break(1)=tau(1)
  l=1

  do i=1, ntau-1
    coef(1,l)=gtau(i)
    coef(3,l)=s(i,4)
    divdif=(gtau(i+1)-gtau(i))/s(i,1)
    z=s(i,5)
    if ( z- 0.5D+00 ) 300, 310, 320

  300   continue

    if ( z == 0.0D+00 ) go to 330
    zeta=gam*z
    onemzt=1.0D+00-zeta
    c=s(i+1,4) / 6.0D+00
    d=s(i,4)*s(i,6)
    l=l+1
    del=zeta*s(i,1)
    break(l)=tau(i)+del
    zt2=zeta**2
    alpha=alph(onemzt)
    factor=onemzt**2*alpha
    coef(1,l)=gtau(i)+divdif*del+s(i,1)**2*(d*onemzt*(factor-1.0D+00) &
      +c*zeta*(zt2-1.0D+00))
    coef(2,l)=divdif+s(i,1)*(d*(1.0D+00-3.0D+00*factor)+c*(3.0D+00*zt2-1.0D+00))
    coef(3,l)=6.0D+00*(d*alpha*onemzt+c*zeta)
    coef(4,l)=6.0D+00*(c-d*alpha)/s(i,1)
    coef(4,l-1)=coef(4,l)-6.0D+00*d*(1.0D+00-alpha)/(del*zt2)
    coef(2,l-1)=coef(2,l)-del*(coef(3,l)-(del/2.0D+00)*coef(4,l-1))
    go to 340

  310   continue

    coef(2,l) = divdif - s(i,1) * ( 2.0D+00 * s(i,4) + s(i+1,4) ) / 6.0D+00
    coef(4,l)=(s(i+1,4)-s(i,4))/s(i,1)
    go to 340

  320 continue

    onemzt=gam*(1.0D+00-z)

    if ( onemzt == 0.0D+00 ) then
      go to 330
    end if

    zeta = 1.0D+00 - onemzt
    alpha=alph(zeta)
    c=s(i+1,4)*s(i,6)
    d=s(i,4)/6.0D+00
    del=zeta*s(i,1)
    break(l+1)=tau(i)+del
    coef(2,l)=divdif-s(i,1)*(2.0D+00*d+c)
    coef(4,l)=6.0D+00*(c*alpha-d)/s(i,1)
    l=l+1
    coef(4,l)=coef(4,l-1)+6.0D+00*(1.0D+00-alpha)*c/(s(i,1)*onemzt**3)
    coef(3,l)=coef(3,l-1)+del*coef(4,l-1)
    coef(2,l)=coef(2,l-1)+del*(coef(3,l-1)+(del/2.0D+00)*coef(4,l-1))
    coef(1,l)=coef(1,l-1)+del*(coef(2,l-1)+(del/2.0D+00)*(coef(3,l-1) &
      +(del/3.0D+00)*coef(4,l-1)))
    go to 340

  330   continue

    coef(2,l) = divdif
    coef(3,l) = 0D+00
    coef(4,l) = 0.0D+00

  340   continue

    l = l + 1
    break(l) = tau(i+1)

  end do

  l = l - 1
  k = 4
  iflag = 1

  return
end