*       ************* FOR SUBROUTINE GAP_B ************************
*       * t=T/T_c0                                                *
*       * deltaB=deltaB/(k_B*T_c0), z(m)=(eps_m+ib)/(k_B*T_c0)    *
*       * N=Number of terms, mu and sigma are impurity parameters *
*       * mu=coherence length/mean free path                      *
*       * required 0 < sigma < 1                                  *
*       ***********************************************************

	SUBROUTINE GAP_B(N, t, deltaB, z, mu, sigma)

	IMPLICIT NONE
	INTEGER          N, i, idime
        parameter(idime=52)
	DOUBLE PRECISION pi
	PARAMETER        (pi=3.141592653589793d0)
	DOUBLE PRECISION HH, HHder, FF(idime), FFder(idime), GG, GGder 
	DOUBLE PRECISION t, mu, sigma, deltaB, z(idime)
	DOUBLE PRECISION y, L, Lder, f, fder, g, gder
	DOUBLE PRECISION hlp1, hlp2, tmp, tmpder
	DOUBLE PRECISION alpha, alphader, at1, at2, at
	DOUBLE PRECISION err, tol
	INTRINSIC        SQRT, ATAN

*       rough estimates for the solution
	deltaB = 1
	y = deltaB/(2*pi) 
*       exact when delta_B=0 or mu=0
	DO 50 i = 1, N
	   z(i) = t*(i-0.5d0)+0.5*mu
 50	CONTINUE
*       lower limit for the integration 
	L = t*N+0.5*mu                    
	HH = 1.0d0
	err = 1.0d0
*       tol is how close to zero the equations must be
	tol = 1.0d-10

*       Newton iteration for delta_B (y) and z(i)
 60	CONTINUE

*       calculate the lower limit for the integration 
*       (for a given y and z_m)
 65	CONTINUE
	f = SQRT(L**2+y**2)
	g = L**2+(1.0-sigma)*y**2
	HH = L-0.5*mu*L*f/g-t*N
	IF (ABS(HH) .GT. tol) THEN
	   HHder = 1.0d0-0.5*mu*(f/g+L**2/(f*g)-2*L**2*f/(g**2)) 
	   L = L - HH/HHder
	   GOTO 65
	END IF
*       end of lower limit calculation
        
	err = ABS(HH)
*       iterate the z(i) variables
	DO 80 i = 1, N
	   hlp1 = SQRT(z(i)**2+y**2)
	   hlp2 = z(i)**2+(1-sigma)*y**2
	   FF(i) = z(i)-0.5*mu*z(i)*hlp1/hlp2-t*(i-0.5d0)
	   FFder(i) = 1.0d0-0.5*mu*(hlp1/hlp2+z(i)**2/(hlp1*hlp2)
     &	              -2*z(i)**2*hlp1/hlp2**2)
	   z(i) = z(i)-FF(i)/FFder(i)
	   IF (ABS(FF(i)) .GT. err) THEN
	      err = ABS(FF(i))
	   END IF
 80	CONTINUE

*       iterate delta_B (y)
	f = SQRT(L**2+y**2)
	g = L**2+(1.0-sigma)*y**2
	hlp1 = 0.5*mu*(L*y*g-2*(1.0-sigma)*y*L*f**2)
	hlp2 = f*g**2-0.5*mu*(f**2*g+L**2*g-2*L**2*f**2)
	Lder = hlp1/hlp2
	fder = (L*Lder+y)/f
	gder = 2*L*Lder+2*(1.0-sigma)*y
	alpha = 1.0-0.5*mu*(f/g+L**2/(f*g)-2*L**2*f/g**2)
	alphader = 0.5*mu*((fder*g-f*gder)/g**2
     &       +(2*L*Lder*f*g-L**2*(fder*g+f*gder))/(f*g)**2
     &       -(2*(2*L*Lder*f+L**2*fder)*g**2-4*L**2*f*g*gder)/g**4)
	at1 = ATAN(L/y)
	at2 = SQRT(1.0-sigma)*ATAN(L/(y*SQRT(1.0-sigma)))
	at = (at1-at2)/sigma
*	IF (sigma .GT. 0.9) THEN
*	   at = ATAN(L/y)/sigma 
*	ELSEIF (sigma .LT. 0.1) THEN
*          at = 0.5*(ATAN(L/y)-L*y/f**2)
*	END IF
	tmp = LOG((L+f)/(2*N))+0.25*pi*mu/(y*(1.0+SQRT(1.0-sigma)))
     &     -0.5*mu*(L/g+at/y)+(-1.0d0/N**2+t**2*L/(f**3*alpha))/24.0
	tmpder = (Lder+fder)/(L+f)
     &     -0.25*pi*mu/(y**2*(1.0+SQRT(1.0-sigma)))
     &	   -0.5*mu*((Lder*g-L*gder)/g**2-at/(y**2)
     &     +((y*Lder-L)/f**2-(1.0-sigma)*(y*Lder-L)/g)/(sigma*y))
     &     +t**2*((Lder*f-3*L*fder)/(alpha*f**4)
     &     -L*alphader/(f**3*alpha**2))/24.0
	GG = 0d0
	GGder = 0d0
	DO 110 i = N, 1, -1
	   GG = GG + (1.0d0/(i-0.5d0)-t/SQRT(z(i)**2+y**2))
	   GGder = GGder + t*y/SQRT(z(i)**2+y**2)**3  
 110	CONTINUE
	GG = GG+tmp
	IF ((err .GT. tol) .OR. (ABS(GG) .GT. tol)) THEN
	   GGder = GGder+tmpder
	   y = y-GG/GGder
	   GOTO 60
	END IF
*       end of Newton iteration

	DO 130 i = 1, N
	   z(i) = z(i)*2*pi
 130	CONTINUE
	deltaB = 2*pi*y
  
	RETURN

	END