c**********************************************************************

      subroutine integral(function,l_min,l_max,answer)

c**********************************************************************

        integer n,m
	parameter(n=3)
        
	real xgauss(n),wgauss(n)
	real l_min,l_max,crit
        real bound
        
	complex function,re1,re2,answer

        common /int_data/ bound
      
        external function
        
      	m=2

        crit=bound+1.0
      	call GAULEG(0.0,1.0,xgauss,wgauss,n)
	call cGs1d(function,l_min,l_max,2,re1,xgauss,wgauss,n)
       
        do while (crit .gt. bound)


          call cGs1d(function,l_min,l_max,2**m,re2,xgauss,wgauss,n)
          if (cabs(re1) .lt. 1.0e-20) then
          
            crit=cabs(re1-re2)/1.0e-20

          else
          
             crit=cabs((re1-re2))/cabs(re1)
            
          endif
        
          m=m+1
          re1=re2

        enddo
        
        answer=re2

c        print*,'int m =',m

	return
	end

c**********************************************************************
      SUBROUTINE cGs1d(Fn, ax, bx, Nx, res, Xgauss, Wgauss, Ngauss )
c**********************************************************************
c
c  This routine computes the gaussian integration of a complex 
c  one-dimensional function Fn over the interval [ax,bx].
c
c----------------------------------------------------------------------
c  Parameters:
c        (ax,bx) - lower and upper limits of the x-dependence of Fn
c             Nx - Number of subintervals for the x-dependence
c         Ngauss - Number of gaussian points to use for each 
c                  subinterval (provided via routine GAULEG)
c         Xgauss - Number of gaussian abcissas to use for each 
c                  subinterval (provided via routine GAULEG)
c         Wgauss - Number of gaussian weights to use for each 
c                  subinterval (provided via routine GAULEG)
c**********************************************************************
                                              
      INTEGER Ngauss, Nx
      COMPLEX Fn, res
      REAL ax, bx, Xgauss(Ngauss), Wgauss(Ngauss)
c...declaration of additional passed parameters
              
      INTEGER I, P
      REAL delx
      COMPLEX sum

      delx = (bx-ax)/Nx

      sum = CMPLX(0.0,0.0)
      DO P=1,Ngauss
        
        DO I=0,Nx-1
          
          sum = sum + Wgauss(P)
     -                * Fn( ( Xgauss(P) + I ) * delx + ax)

          
        ENDDO

      ENDDO

      res = sum * delx

      RETURN
      END




c**********************************************************************
      SUBROUTINE GAULEG(X1,X2,X,W,N)
c**********************************************************************
c  This is a numerical recipes routine for computing the
c  weights and abscissas for Gauss-Legendre quadrature.
c**********************************************************************
      DOUBLE PRECISION   EPS, XM, XL, Z, pi, P1, P2, P3, PP, Z1
      REAL X1,X2,X(N),W(N)
      INTEGER I, J, M, N
      PARAMETER (EPS=3.D-14)
      DATA pi /3.141592741012573D0/

      M=(N+1)/2
      XM=0.5D0*(X2+X1)
      XL=0.5D0*(X2-X1)
      DO 12 I=1,M
        Z=COS(pi*(I-.25D0)/(N+.5D0))
1       CONTINUE
          P1=1.D0
          P2=0.D0
          DO 11 J=1,N
            P3=P2
            P2=P1
            P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
11        CONTINUE
          PP=N*(Z*P1-P2)/(Z*Z-1.D0)
          Z1=Z
          Z=Z1-P1/PP
        IF(ABS(Z-Z1).GT.EPS)GO TO 1
        X(I)=XM-XL*Z
        X(N+1-I)=XM+XL*Z
        W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
        W(N+1-I)=W(I)
12    CONTINUE
      RETURN
      END