c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine cubic_spline(nk,nu,k,z,valid)
      implicit none
      include 'max.inc'
c
c     Purpose: to find the nk values of "z" coefficients that are needed to compute
c     the nk-1 cubic spline functions.
c
c     Inputs:
c       + nk: number of nu/k points
c       + nu: array of wavenumber points
c       + k: array of absorption coefficient points
c
c     Outputs:
c       + z: nk values of z coefficients
c       + valid: 1 if solution found; 0 otherwise.
c

c     input
      integer nk
      double precision nu(1:kmax)
      double precision k(1:kmax)
c     output
      double precision z(1:kmax)
      integer valid
c     temp
      integer i
      double precision alpha(1:kmax)
      double precision beta(1:kmax)
      double precision gamma(1:kmax)
      double precision delta(1:kmax)
      double precision a(1:kmax)
      double precision b(1:kmax)
      double precision c(1:kmax)
      double precision d(1:kmax)
      double precision sol(1:kmax)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine cubic_spline'

c     computing the n values of z coefficients
      z(1)=0.0D+0
      z(nk)=0.0D+0
      do i=2,nk-1
         if (i.eq.2) then
            alpha(i)=0.0D+0
         else
            alpha(i)=nu(i)-nu(i-1)
         endif
         beta(i)=2.0D+0*(nu(i+1)-nu(i-1))
         if (i.lt.nk-1) then
            gamma(i)=nu(i+1)-nu(i)
         else
            gamma(i)=0.0D+0
         endif
         delta(i)=6.0D+0*((k(i+1)-k(i))/(nu(i+1)-nu(i))
     &        -(k(i)-k(i-1))/(nu(i)-nu(i-1)))
c     Debug
c         if (i.eq.2) then
c            write(*,*) 'delta(',i,')=',delta(i)
c            write(*,*) 'k(',i-1,')=',k(i-1)
c            write(*,*) 'k(',i,')=',k(i)
c            write(*,*) 'k(',i+1,')=',k(i+1)
c            write(*,*) 'nu(',i-1,')=',nu(i-1)
c            write(*,*) 'nu(',i,')=',nu(i)
c            write(*,*) 'nu(',i+1,')=',nu(i+1)
c         endif
c     Debug
      enddo
      do i=1,nk-2
         a(i)=alpha(i+1)
         b(i)=beta(i+1)
         c(i)=gamma(i+1)
         d(i)=delta(i+1)
      enddo
      call tridiag(nk-2,a,b,c,d,sol,valid)
      if (valid.eq.0) then
         write(*,*) 'Error from routine ',label(1:strlen(label)),':'
         write(*,*) 'solution NOT found'
         stop
      endif
      do i=2,nk-1
         z(i)=sol(i-1)
      enddo
c     Debug
c      do i=2,nk-1
c         write(*,*) 'z(',i,')=',z(i)
c      enddo ! i
c     Debug

      return
      end