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 cumulative(n,nu,k,eps1,eps2,k0,sit,z,g,delta_g)
      implicit none
      include 'max.inc'
c
c     This routine will compute the cumulative g(k0) in
c     the specified spectral interval
c
c     Inputs:
c       + n: size of "nu" and "k" arrays
c       + nu: array of n nu values
c       + eps1,eps2: accuracy criterion (from KSPECTRUM)
c       + k: array of corresponding ka values
c       + k0: value of k for which g has to be computed
c       + sit: spline interpolation type (1:linear, 2:cubic)
c       + z: values of the n coefficients "z" that are needed for cubic spline interpolation
c
c     Outputs:
c       + g: value of g(k0)
c       + delta_g: uncertainty over g(k0)
c
      integer n,i,in,sit
      double precision delta_nu
      double precision nu(1:kmax),k(1:kmax)
      double precision z(1:kmax)
      double precision eps1,eps2
      double precision k0,nu1,nu2,sum_nu,sum_dnu,g,delta_g
      double precision k1p,k1m,k2p,k2m,nua,nub,dnug,dnud
      double precision a,b,c,d
      integer nr,nuf
      double precision x(1:3)
      double precision epsg
      double precision k_min,k_max
      logical is_nan
c     Debug
c      double precision S
c     Debug
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine cumulative'

      call minmax(n,k,k_min,k_max)
      if (k0.lt.k_min) then
         g=0.0D+0
         goto 666
      endif
      if (k0.ge.k_max) then
         g=1.0D+0
         goto 666
      endif
c     Debug
c      write(*,*) 'ok00'
c     Debug
c     Computation of g(k0)
      sum_nu=0.0D+0
      sum_dnu=0.0D+0
      in=0
      delta_nu=nu(n)-nu(1)
      epsg=1.0D-10
c     Debug
c      write(*,*) 'ok01'
c     Debug

      if (n.lt.2) then
         call error(label)
         write(*,*) 'Number of points:',n
         write(*,*) 'is not enough !'
         stop
      endif
c     Debug
c      write(*,*) 'ok02'
c     Debug

      if (n.gt.kmax) then
         call error(label)
         write(*,*) 'n=',n
         write(*,*) 'is greater than kmax=',kmax
         write(*,*) 'please adjust its value in max.inc'
         write(*,*) 'and recompile from scratch'
         stop
      endif
c     Debug
c      write(*,*) 'ok03'
c     Debug
      
      if (k(1).le.k0) then
         nu1=nu(1)
c     Debug
c         write(*,*) 'nu1=',nu1
c     Debug
c     Debug
         call test_nan(nu1,is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'nu1=',nu1
            write(*,*) '@1'
            stop
         endif
c     Debug
         in=1
         dnug=0.0D+0
      endif
c     Debug
c      write(*,*) 'ok04, n=',n
c     Debug

      do i=2,n
c     Debug
c      write(*,*) 'i=',i
c     Debug
         if ((k(i-1).gt.k0).and.(k(i).le.k0)) then ! descending intersection
c     Debug
c            write(*,*) 'descending intersection: in'
c     Debug
            if (sit.eq.1) then
               nu1=nu(i-1)+(nu(i)-nu(i-1))/(k(i)-k(i-1))*(k0-k(i-1))
c     Debug
c               write(*,*) 'nu1=',nu1
c     Debug
c     Debug
               call test_nan(nu1,is_nan)
               if (is_nan) then
                  call error(label)
                  write(*,*) 'nu1=',nu1
                  write(*,*) '@2'
                  stop
               endif
c     Debug

            else
c     Debug
c               write(*,*) 'ok000'
c               write(*,*) 'nu(',i-1,')=',nu(i-1)
c               write(*,*) 'nu(',i,')=',nu(i)
c               write(*,*) 'k0=',k0
c               write(*,*) 'k(',i-1,')=',k(i-1)
c               write(*,*) 'k(',i,')=',k(i)
c               write(*,*) 'z(',i-1,')=',z(i-1)
c               write(*,*) 'z(',i,')=',z(i)
c     Debug
               call abcd_coeffs(nu(i-1),nu(i),k0,k(i-1),k(i),
     &              z(i-1),z(i),a,b,c,d)
c     Debug
c               write(*,*) 'ok001'
c               write(*,*) 'a=',a
c               write(*,*) 'b=',b
c               write(*,*) 'c=',c
c               write(*,*) 'd=',d
c     Debug
               call findroots_3rd(a,b,c,d,nr,x)
c     Debug
c               write(*,*) 'ok002'
c     Debug
               if (nr.eq.0) then
                  call error(label)
                  write(*,*) 'Number of solutions=',nr
                  stop
               endif
c     Debug
c               write(*,*) 'ok003'
c     Debug
               call identify_solution(nr,x,nu(i-1),nu(i),
     &              k0,k(i-1),k(i),z(i-1),z(i),nu1,nuf)
c     Debug
c               write(*,*) 'nu1=',nu1
c     Debug
c     Debug
               call test_nan(nu1,is_nan)
               if (is_nan) then
                  call error(label)
                  write(*,*) 'nu1=',nu1
                  write(*,*) '@3'
                  stop
               endif
c     Debug
c     Debug
c               write(*,*) 'ok004'
c     Debug
               if (nuf.eq.0) then
                  call dichotomy(nu(i-1),nu(i),k(i-1),k(i),k0,
     &                 z(i-1),z(i),nu1)
               endif
c     Debug
c               write(*,*) 'ok005'
c     Debug
c     Debug
               if ((nu1.lt.nu(i-1)).or.(nu1.gt.nu(i))) then
                  call error(label)
                  write(*,*) 'nu1=',nu1,' is out of bounds:'
                  write(*,*) 'nu(',i-1,')=',nu(i-1)
                  write(*,*) 'nu(',i,')=',nu(i)
                  write(*,*) 'nuf=',nuf
                  stop
               endif
c     Debug
               
            endif
            k1m=k(i-1)*(1.0D+0-eps1)
            k1p=k(i-1)*(1.0D+0+eps1)
            k2m=k(i)*(1.0D+0-eps1)
            k2p=k(i)*(1.0D+0+eps1)
            nua=nu(i-1)+(k0/(1.0D+0+eps2)-k1p)
     &           *(nu(i)-nu(i-1))/(k2p-k1p)
            nub=nu(i-1)+(k0/(1.0D+0-eps2)-k1m)
     &           *(nu(i)-nu(i-1))/(k2m-k1m)
c     Debug
            if (nua.lt.nub) then
               call error(label)
               write(*,*) 'nua=',nua
               write(*,*) 'nub=',nub
               write(*,*) 'nu(',i-1,')=',nu(i-1),' nu(',i,')=',nu(i)
               write(*,*) 'k(',i-1,')=',k(i-1),' k(',i,')=',k(i)
               write(*,*) 'k1m=',k1m,' k1p=',k1p
               write(*,*) 'k2m=',k2m,' k2p=',k2p
               write(*,*) 'k0=',k0
               write(*,*) 'eps1=',eps1,' eps2=',eps2
               stop
            endif
c     Debug
            dnug=nua-nub
            if (in.eq.1) then
               call error(label)
               write(*,*) 'error 1: nu1 already fixed, i=',i
               write(*,*) 'k0=',k0,' k(',i-1,')=',k(i-1)
     &              ,' k(',i,')=',k(i)
               stop
            else
               in=1
            endif
c     Debug
c            write(*,*) 'descending intersection: out'
c     Debug
         endif ! descending intersection

         if ((k(i-1).le.k0).and.(k(i).gt.k0)) then ! ascending intersection
c     Debug
c            write(*,*) 'ascending intersection: in'
c     Debug
            if (sit.eq.1) then
               nu2=nu(i-1)+(nu(i)-nu(i-1))/(k(i)-k(i-1))*(k0-k(i-1))
c     Debug
               call test_nan(nu2,is_nan)
               if (is_nan) then
                  call error(label)
                  write(*,*) 'nu2=',nu2
                  write(*,*) '@1'
                  stop
               endif
c     Debug

            else
               call abcd_coeffs(nu(i-1),nu(i),k0,k(i-1),k(i),
     &              z(i-1),z(i),a,b,c,d)
               call findroots_3rd(a,b,c,d,nr,x)
               if (nr.eq.0) then
                  call error(label)
                  write(*,*) 'Number of solutions=',nr
                  stop
               endif
               call identify_solution(nr,x,nu(i-1),nu(i),
     &              k0,k(i-1),k(i),z(i-1),z(i),nu2,nuf)
c     Debug
               call test_nan(nu1,is_nan)
               if (is_nan) then
                  call error(label)
                  write(*,*) 'nu2=',nu2
                  write(*,*) '@3'
                  stop
               endif
c     Debug
               if (nuf.eq.0) then
                  call dichotomy(nu(i-1),nu(i),k(i-1),k(i),k0,
     &                 z(i-1),z(i),nu2)
               endif
c     Debug
               if ((nu2.lt.nu(i-1)).or.(nu2.gt.nu(i))) then
                  call error(label)
                  write(*,*) 'nu2=',nu2,' is out of bounds:'
                  write(*,*) 'nu(',i-1,')=',nu(i-1)
                  write(*,*) 'nu(',i,')=',nu(i)
                  stop
               endif
c     Debug
               
            endif
            k1m=k(i-1)*(1.0D+0-eps1)
            k1p=k(i-1)*(1.0D+0+eps1)
            k2m=k(i)*(1.0D+0-eps1)
            k2p=k(i)*(1.0D+0+eps1)
            nua=nu(i-1)+(k0/(1.0D+0+eps2)-k1p)
     &           *(nu(i)-nu(i-1))/(k2p-k1p)
            nub=nu(i-1)+(k0/(1.0D+0-eps2)-k1m)
     &           *(nu(i)-nu(i-1))/(k2m-k1m)
c     Debug
            if (nub.lt.nua) then
               call error(label)
               write(*,*) 'nua=',nua
               write(*,*) 'nub=',nub
               write(*,*) 'nu(',i-1,')=',nu(i-1),' nu(',i,')=',nu(i)
               write(*,*) 'k1m=',k1m,' k1p=',k1p
               write(*,*) 'k2m=',k2m,' k2p=',k2p
               write(*,*) 'k0=',k0
               write(*,*) 'eps1=',eps1,' eps2=',eps2
               stop
            endif
c     Debug
            dnud=nub-nua
            if (in.eq.0) then
               call error(label)
               write(*,*) 'error 1: nu1 not fixed'
               stop
            else
               in=0
            endif
            if (in.eq.0) then
               sum_nu=sum_nu+nu2-nu1
               sum_dnu=sum_dnu+dnud+dnug
c     Debug
               call test_nan(sum_nu,is_nan)
               if (is_nan) then
                  call error(label)
                  write(*,*) '@1'
                  write(*,*) 'sum_nu=',sum_nu
                  write(*,*) 'nu2-nu1=',nu2-nu1
                  write(*,*) 'nu1=',nu1
                  write(*,*) 'nu2=',nu2
                  stop
               endif
c     Debug
            endif

c     Debug
c            write(*,*) 'ascending intersection: out'
c     Debug
         endif ! ascending intersection

      enddo                     ! i=1,n
c     Debug
c      write(*,*) 'ok05'
c     Debug

      if (k(n).le.k0) then
         if (in.eq.0) then
            call error(label)
            write(*,*) 'error 2: nu1 not fixed'
            write(*,*) 'n=',n
            write(*,*) 'k(',n,')=',k(n)
            write(*,*) 'k0=',k0
            write(*,*) 'k(',n-1,')=',k(n-1)
            write(*,*) 'nu(',n,')=',nu(n)
            write(*,*) 'nu(',n-1,')=',nu(n-1)
            write(*,*) 'sit=',sit
            stop
         else
            in=0
         endif
         if (in.eq.0) then
            nu2=nu(n)
c     Debug
            call test_nan(nu2,is_nan)
            if (is_nan) then
               call error(label)
               write(*,*) 'nu2=',nu2
               write(*,*) '@2'
               stop
            endif
c     Debug

            sum_nu=sum_nu+nu2-nu1
            sum_dnu=sum_dnu+dnug
c     Debug
            call test_nan(sum_nu,is_nan)
            if (is_nan) then
               call error(label)
               write(*,*) '@2'
               write(*,*) 'sum_nu=',sum_nu
               write(*,*) 'nu2-nu1=',nu2-nu1
               write(*,*) 'nu1=',nu1
               write(*,*) 'nu2=',nu2
               stop
            endif
c     Debug
         endif
      endif
      
      g=sum_nu/delta_nu
      delta_g=sum_dnu/delta_nu

      call test_nan(g,is_nan)
      if (is_nan) then
         call error(label)
         write(*,*) 'g=',g
         write(*,*) 'sum_nu=',sum_nu
         write(*,*) 'delta_nu=',delta_nu
         stop
      endif

 666  continue
      return
      end