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 g_inversion(lay,band,quad,g0,Np,gac,kac,epsilon,
     &     n,nu,k,eps1,eps2,sit,z,k0,delta_k0,err_f)
      implicit none
      include 'max.inc'
c
c     Dichotomy procedure for inversion of g values
c     
c     Inputs:
c     + g0: value of g to invert
c     + Np: number of points in acceleration arrays
c     + gac and kac: acceleration arrays
c     + epsilon: max error over the value of k0
c     + n: number of nu and ka values
c     + nu and k : nu and ka arrays (n values)
c     + eps1 & eps2: accuracy criterion from KSPECTRUM
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     Ouput:
c     + k0: value of ka corresponding to the value of g provided
c     + err_f: error over k0
c
      integer n,Np,i,j,p,niter,band,quad,lay
      double precision g0,epsilon,k0,delta_k0,err_f,err_x
      double precision nu(1:kmax),k(1:kmax)
      double precision eps1,eps2
      integer sit
      double precision z(1:kmax)
      double precision gac(0:Npmax),kac(0:Npmax),delta_g
      double precision k1,k2,g1,g2,km,gm,f,k_estf,dk
      double precision g_1,g_2,gref
      double precision last(1:5),moy,std
      double precision epsilon_x,epsilon_f
      integer cl,d
      logical pf,keep_going,solution_found,cl_changed,no_variation_of_x
      logical is_nan
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine g_inversion'

      d=0
      epsilon_f=epsilon ! maximum variation of f (g)
      epsilon_x=1.0D-3 ! maximum variation of x (k)
c     Debug
c      if ((band.eq.73).and.(quad.eq.1)) then
c         write(*,*) 'g0=',g0
c         d=1
c      endif
c      do i=0,Np
c         write(*,*) 'gac(',i,')=',gac(i),kac(i)
c      enddo
c     Debug
c     values of g0 lower than gac(0):
c      if (g0.le.gac(0)) then
c         km=0.0D+0
c         err=0.0D+0
c         goto 321
c      endif

      pf=.false.
c     acceleration
      do i=1,Np-1
c     Debug
c         write(*,*) 'i=',i
c         write(*,*) gac(i-1),g0,gac(i),pf
c         if (i.gt.20) then
c            stop
c         endif
c     Debug
         if ((g0.ge.gac(i-1)).and.(g0.le.gac(i))) then
            p=i
            pf=.true.
c     Debug
            if (d.eq.1) then
               write(*,*) 'ok'
            endif
c     Debug
            goto 123
         endif
      enddo
 123  continue
      if (.not.pf) then
         call error(label)
         write(*,*) 'acceleration interval was not found'
         write(*,*) 'g0=',g0
         do i=1,Np
            write(*,*) gac(i-1),gac(i)
         enddo ! i
c         write(*,*) gac(Np-1),gac(Np)
         stop
      endif
c     Debug
      if (d.eq.1) then
         write(*,*) 'gac(',p-1,')=',gac(p-1),' g0=',g0,
     &        ' gac(',p,')=',gac(p)
      endif
c     Debug
c     dichotomy procedure
      k1=kac(p-1)
      k2=kac(p)
      g1=gac(p-1)
      g2=gac(p)
c     Debug
      if (d.eq.1) then
         write(*,*) 'k1=',k1,' k2=',k2
         write(*,*) 'g1=',g1,' g2=',g2
         write(*,*) 'g0=',g0
         call cumulative(n,nu,k,eps1,eps2,k1,sit,z,gm,delta_g)
         write(*,*) 'g(k1)=',gm
         call cumulative(n,nu,k,eps1,eps2,k2,sit,z,gm,delta_g)
         write(*,*) 'g(k2)=',gm
         write(*,*) 'epsilon_f=',epsilon_f
         write(*,*) 'epsilon_x=',epsilon_x
      endif
c     Debug
      err_x=k2-k1
      err_f=1.0D+0
      niter=0
      cl=0
      do j=1,5
         last(j)=dble(j)
      enddo
      keep_going=.true.
      solution_found=.false.
      cl_changed=.false.
      no_variation_of_x=.false.
      do while (keep_going)
         niter=niter+1
         km=(k1+k2)/2.0D+0
         call cumulative(n,nu,k,eps1,eps2,km,sit,z,gm,delta_g)
c     
         if (niter.le.5) then
            last(niter)=gm
         else
            do j=5,2,-1
               last(j-1)=last(j)
            enddo
            last(5)=gm
         endif
         call std5(last,moy,std,cl)

c     Debug
c         if (cl.eq.1) then
c            write(*,*) 'niter=',niter
c            do j=1,5
c               write(*,*) 'last(',j,')=',last(j)
c            enddo
c            write(*,*) 'moy=',moy
c         endif
c     Debug

c     
         if (g0.eq.0.0D+0) then
            gref=gm
         else
            gref=g0
         endif
         err_f=dabs(gm-g0)/dabs(gref)
         err_x=dabs(g2-g1)
c     Debug
         if (d.eq.1) then
            write(*,*) 'niter=',niter,' k1=',k1,' k2=',k2,' gm=',gm
            write(*,*) 'gm=',gm,' g0=',g0,' err_f=',err_f,' cl=',cl
         endif
c     Debug
         if (gm.gt.g0) then
            k2=km
            g2=gm
         else
            k1=km
            g1=gm
         endif
         if (niter.gt.Niter_max) then
            write(*,*) 'g_inversion failed after',niter,' iterations'
            write(*,*) 'lay=',lay,' band=',band,' quad=',quad
            stop
         endif

         if (err_f.lt.epsilon_f) then
            keep_going=.false.
            solution_found=.true.
         endif
         if (cl.ne.0) then
            keep_going=.false.
            cl_changed=.true.
         endif
         if (err_x.lt.epsilon_x) then
            keep_going=.false.
            no_variation_of_x=.true.
         endif

      enddo

 321  continue
      k0=km

c     Debug
      if (d.eq.1) then
         write(*,*) 'solution_found=',solution_found
         write(*,*) 'cl_changed=',cl_changed
         write(*,*) 'no_variation_of_x=',no_variation_of_x
         write(*,*) 'k0=',k0
      endif
c     Debug

c     Debug
      call test_nan(k0,is_nan)
      if (is_nan) then
         call error(label)
         write(*,*) 'k0=',k0
         stop
      endif
c     Debug

c     computation of delta_k0
      dk=1.0D-5
      k_estf=k0
      call cumulative(n,nu,k,eps1,eps2,k_estf-dk,sit,z,g_1,delta_g)
      call cumulative(n,nu,k,eps1,eps2,k_estf+dk,sit,z,g_2,delta_g)
      f=(g_2-g_1)/(2.0D+0*dk)
      delta_k0=delta_g/f
c     Debug
      if (d.eq.1) then
         stop
      endif
c     Debug

      return
      end