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 analyse_g(band,nk,nu,k,epsilon1,epsilon2,
     &     k_min,k_max,sit,z,ns,g_step)
      implicit none
      include 'max.inc'
c
c     Purpose: to detect steps in the cumulated function of the provided spectrum
c
c     Inputs:
c       + band: index of the narrowband spectral interval
c       + nk: number of points in the k(nu) spectrum
c       + nu: values of nu
c       + k: associated values of k
c       + epsilon1,epsilon2: accuracy criterion (from KSPECTRUM)
c       + k_min: minimum value in the "k" array
c       + k_max: maximum value in the "k" array
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       + ns: number of steps
c       + g_step: values of "g" for the "ns" steps
c
      integer strlen
      integer nk,Np,i,j,Nps,ns,i0,j0,fd,je,sit
      double precision epsilon1,epsilon2
      double precision nu(1:kmax),k(1:kmax)
      double precision z(1:kmax)
      double precision k_min,k_max
      double precision kp(1:Npmax),gp(1:Npmax),delta_g
      double precision g_step(1:Nsmx)
      double precision v1,v2,vmin,eps1,eps2,epsm1,epsm2
      integer band
      double precision klo,khi
c     Debug
      double precision dk,fp,g2,dg,delta_g2
c     Debug
c     label
      character*(Nchar_mx) label
      label='subroutine analyse_g'

c     Debug
c      open(20,file='./results/k.txt')
c      do i=1,nk
c         write(20,*) nu(i),k(i)
c      enddo ! i
c      close(20)
c     Debug

c     Debug
      open(21,file='./results/f_function.txt')
      open(22,file='./results/g_function.txt')
c     Debug
c     Computing function g(k) for Np points
      Np=10000
      if (Np.gt.Npmax) then
         call error(label)
         write(*,*) 'Np=',Np
         write(*,*) 'while Npmax=',Npmax
         stop
      endif
c     Debug
c      write(*,*) 'Np=',Np
c     Debug
      klo=k_min
      khi=k_max
c     ------------------------------------------------------------------
c     This is the part you might want to modify
c      klo=1.0D-7
c      khi=1.0D-6
c     ------------------------------------------------------------------
      do i=1,Np
         kp(i)=klo+(khi-klo)/(Np-1)*(i-1)
         call cumulative(nk,nu,k,epsilon1,epsilon2,kp(i),sit,z,
     &        gp(i),delta_g)
         dk=1.0D-4
         call cumulative(nk,nu,k,epsilon1,epsilon2,kp(i)+dk,sit,z,
     &        g2,delta_g2)
         dg=g2-gp(i)
         fp=dg/dk
         write(21,*) kp(i),fp
         write(22,*) kp(i),gp(i)
      enddo ! i
      close(21)
      close(22)

c     Detecting steps in function g(k)
      Nps=10
c      epsm1=2.0D-4
c      epsm2=5.0D-4
      epsm1=5.0D-4
      epsm2=1.0D-3
      ns=0
      i0=1

 111  continue
c     Debug
c      write(*,*) 'i0=',i0,' Np-Nps=',Np-Nps
c     Debug
      if (i0.gt.Np-Nps) then
         goto 333
      endif
      do i=i0,Np-Nps
         v1=gp(i)
         v2=gp(i+Nps)
         call find_min(v1,v2,vmin)
c     Debug
c         write(*,*) 'vmin=',vmin
c     Debug
         eps1=dabs(v2-v1)/dabs(vmin)
         if (eps1.lt.epsm1) then
            ns=ns+1
            j0=i
            v1=gp(j0)
            fd=0
            goto 112
         endif
      enddo                     ! i
      goto 333

 112  continue
      do j=j0,Np
         v2=gp(j)
         call find_min(v1,v2,vmin)
         eps2=dabs(v2-v1)/dabs(vmin)
         if (eps2.gt.epsm2) then
            fd=1
            je=j-1
            goto 222
         endif
      enddo                     ! j
 222  continue
      if (fd.eq.1) then
         g_step(ns)=(gp(j0)+gp(je))/2.0D+0
         i0=je+1
         goto 111
      else
         g_step(ns)=(gp(j0)+gp(Np))/2.0D+0
         goto 333
      endif

 333  continue

      return
      end