c     kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum
c     Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet
c     
c     This file must be used under the terms of the CeCILL license.
c     This source file is licensed as described in the file COPYING, which
c     you should have received as part of this distribution.  The terms
c     are also available at
c     http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
c     
      subroutine lorentz_tabulation(eps_max,x_min,x_max,file_out)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c     
c     Purpose: to tabulate Lorentz profile
      integer a0,a25,a50,a75,a100
      integer strlen
      integer ndiv
      integer it,Nit_max
      double precision x1,x1p,x2,x,pp
      double precision x_min,x_max,dx0
      double precision emax,eps_max
      double precision epsm,xs
      double precision apc
      character*(Nchar_mx) file_out
      
      ndiv=50
      dx0=2.0D+0
      Nit_max=100
      pp=0.99D+0
      
      a0=0
      a25=0
      a50=0
      a75=0
      a100=0
      apc=0

      open(11,file=file_out(1:strlen(file_out)))

      x1p=x_min
      do while (x1p.lt.x_max)
         x1=x1p
c     set variable emax
         if (dabs(x1).le.1.0D+0) then
            emax=eps_max/1.5D+1
         else
            emax=eps_max
         endif

         epsm=0.0D+0
         it=0
         x=x1+dx0
         xs=x1
         do it=1,Nit_max
            call find_epsml(x1,x,epsm)
            if (epsm.lt.emax*pp) then
               xs=x
               x=x+x-x1
            endif
            if (epsm.gt.emax) then
               x=(xs+x)/2.0D+0
            endif
            if ((epsm.gt.emax*pp).and.(epsm.lt.emax)) then
               x2=x
               goto 123
            endif
c     Debug
c            write(*,*) 'it=',it,' x1=',x1,' x=',x,' xs=',xs,
c     &           ' epsm=',epsm
c     Debug
         enddo ! it
         if (it.ge.Nit_max) then
            write(*,*) 'error: Nit_max reached'
            stop
         endif

 123     continue
         x1p=x1+(x2-x1)/ndiv
         dx0=1.1D+0*(x2-x1)
         
         if ((x1.lt.0.0D+0)
     &        .and.(x2.gt.0.0D+0)
     &        .and.(epsm.lt.emax)) then
            x2=0.0D+0
            if (x1p.gt.0.0D+0) then
               x1p=0.0D+0
            endif
         endif
         if ((x1p.gt.x_max)
     &        .and.(epsm.lt.emax)) then
            x1p=x_max
         endif
c     Debug
c         write(*,*) x1,x1p,x2!,epsm
c     Debug
         write(11,50) x1,x1p,x2
         apc=(x1p-x_min)/(x_max-x_min)
         call print_tabulation_progress(apc,a0,a25,a50,a75,a100)
      enddo

      close(11)

      write(*,*) file_out(1:strlen(file_out)),
     &     ' was successfully generated' 

      return
      end




      subroutine find_epsml(x1,x,epsm)
      implicit none
      include 'max.inc'
      double precision fL,fs
      double precision x,x1,epsm,xl
      double precision eps(1:Nmax_vt)
      integer N,i,imax

c     Debug
      if (x.lt.x1) then
         write(*,*) 'This is routine find_epsml'
         write(*,*) 'x1=',x1,' x=',x
      endif
c     Debug
      N=nint((x-x1)*1.0D+4)
      if (N.gt.Nmax_vt) then
         do while (N.gt.Nmax_vt)
            N=N/10
         enddo
      endif
      
      do i=1,N
         xl=x1+(x-x1)/(N-1)*(i-1)
         call fstarL(x1,x,xl,fs)
         eps(i)=dabs(fL(xl)-fs)/fL(xl)
      enddo
      call find_max_vt(eps,N,epsm,imax)

      return
      end
   

      
      double precision function fL(x)
      implicit none
      double precision x

      fL=1.0D+0/(1.0D+0+x**2)

      return
      end


      subroutine fstarL(x1,x2,x,fs)
      implicit none
      double precision fL
      double precision x1,x2,x,fs,f1,f2
      
      f1=fL(x1)
      f2=fL(x2)
      fs=f1+(x-x1)*(f2-f1)/(x2-x1)
      
      return
      end