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 acceleration(interface,band,n,nu,k,eps1,eps2,
     &     Np,k_max,sit,z,kac,gac)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute the "gac" acceleration array
c     
c     Inputs:
c       + interface: interface index
c       + band: band index
c       + n: size of "nu" and "k" arrays
c       + nu: array of n nu values
c       + k: array of corresponding ka values
c       + eps1,eps2: accuracy criterion (from KSPECTRUM)
c       + Np: number of points in the acceleration array
c       + k_max: maximum value of k
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       + kac: values of k for the Np acceleration points
c       + gac: values of g(k) for the Np acceleration points
c
c     I/O
      integer interface,band
      integer n,Np,sit
      double precision nu(1:kmax),k(1:kmax)
      double precision eps1,eps2
      double precision z(1:kmax)
      double precision k_max
      double precision gac(0:Npmax),kac(0:Npmax)
c     temp
      integer i
      double precision delta_g
      character*(Nchar_mx) db_file
      logical is_nan
c     label
      integer strlen
      character*(Nchar_mx) label
      label='acceleration'

      do i=0,Np-1
         kac(i)=k_max/dble(Np-1)*dble(i)
      enddo
      
      do i=0,Np-1
         call cumulative(n,nu,k,eps1,eps2,kac(i),sit,z,gac(i),delta_g)
      enddo
c     Debug
c      do i=0,Np-1
c         write(*,*) 'gac(',kac(i),')=',gac(i)
c      enddo
c      stop
c     Debug

c     check the final value is equal to 1
c      if (dabs(gac(Np-1)).ne.1.0D+0) then
c         write(*,*) 'gac(',Np-1,')',gac(Np-1)
c         do i=0,Np-1
c            gac(i)=gac(i)/gac(Np-1)
c         enddo ! i
c      endif

      do i=0,Np-1
         call test_nan(gac(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'Before normalization'
            write(*,*) 'gac(',i,')=',gac(i)
            stop
         endif
      enddo ! i

      if (dabs(dabs(gac(Np-1))-1.0D+0).gt.1.0D-1) then
         write(*,*) 'Warning from routine ',
     &        label(1:strlen(label)),' :'
         write(*,*) 'gac(',Np-1,')=',gac(Np-1)
c         if (dabs(dabs(gac(Np-1)-1.0D+0)).lt.1.0D-2) then
         if (gac(Np-1).eq.0.0D+0) then
            write(*,*) 'could not normalize -> error !'
            write(*,*) 'interface=',interface
            write(*,*) 'band=',band
            db_file='./results/debug.txt'
            open(10,file=db_file(1:strlen(db_file)))
            write(10,*) 'interface=',interface
            write(10,*) 'band=',band
            write(10,*)
            write(10,*) 'Np=',Np
            do i=0,Np-1
               write(10,*) 'kac(',i,')=',kac(i),' gac(',i,')=',gac(i)
            enddo ! i
            write(10,*)
            write(10,*) 'n=',n
            do i=1,n
               write(10,*) 'nu(',i,')=',nu(i),' k(',i,')=',k(i)
            enddo ! i
            close(10)
            write(*,*) 'Debugging information written to file:'
            write(*,*) db_file(1:strlen(db_file))
            stop            
         else if (gac(Np-1).ne.1.0D+0) then
            do i=0,Np-1
               gac(i)=gac(i)/gac(Np-1)
            enddo               ! i
            write(*,*) 'normalization OK'
         endif
      endif

      do i=0,Np-1
         call test_nan(gac(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'After normalization'
            write(*,*) 'gac(',i,')=',gac(i)
            stop
         endif
      enddo ! i

      if (gac(Np-1).ne.1.00000000000000000D+0) then
         gac(Np-1)=1.00000000000000000D+0
      endif

      return
      end