c **********************************************************************
      subroutine DEint(FUNC,a,b,n,sum)
c **********************************************************************
c
c     This subroutine calculate the integral
c       sum = integral_a^b FUNC(x) dx
c     by using the Double Exponential formula (Takahashi-Mori).
c
c       sum= integral_{-infinity}^infinity FUNC(phi(t)) d(phi)/dt dt
c       phi = (b-a)/2*tanh(pi/2*sinh(t))+(b+a)/2
c
c **********************************************************************
      implicit real*8 (a-h,o-z)
      external FUNC
c
      tmax=4.7d0
      dt=tmax/n
      pi=4.0d0*atan(1.0d0)
c
      sum = 0.0d0
      do i = -n,n
         t=i*dt
         s=sinh(t)
         phi=(b-a)/2*tanh(pi/2*s)+(b+a)/2
         dphidt=pi/2*cosh(t)/(cosh(pi/2*s))**2*(b-a)/2
         sum=sum+FUNC(phi)*dphidt
      end do
         sum = sum*dt
c
      end