c Subroutine to calculate the bcs gap function in the pure case 
c (s-wave superconductor or 3He-B)
c Essentially an exact calculation, see 
c http://boojum.hut.fi/research/theory/qc/bcsgap.html
c The idea is to sum exactly the first, say 10, terms in the matsubara
c sum, and estimate the rest by an integral and a Euler-McLaurin correction term
c input:
c   temp= temperature/T_c
c output:
c   gap= weak-coupling energy gap/T_c
      subroutine bcsgap(temp,gap)
      implicit none
      double precision temp,gap,pii2,tn
      double precision fprime,eps,sq,gaps,f,x
      integer iter,ifreq,nsplit
      parameter(pii2=2.*3.1415926535897932,nsplit=10)
      iter=0
      if(temp.ge.1.)then
        gap=0.
        goto 200
      end if
      gap=1.7638*sqrt(1.-temp)/pii2
10    iter=iter+1
      gaps=gap**2
      tn=temp*float(nsplit)
      sq=sqrt(tn**2+gaps)
      f=log((tn+sq)/(2.*float(nsplit)))
     x -(1./float(nsplit)**2-float(nsplit)*(temp/sq)**3)/24.
      fprime=gap/(sq*(tn+sq))-tn*temp**2*gap/(8*sq**5)
      do 100 ifreq=1,nsplit
        eps=float(ifreq)-0.5
        sq=sqrt((temp*eps)**2+gaps)
        f=f+1./eps-temp/sq
        fprime=fprime+temp*gap/sq**3
100   continue
      x=f/fprime
      gap=gap-x
      if(abs(x).ge.0.00000001/pii2)goto 10
      gap=pii2*gap
200   continue
c     write(6,1)gap,iter,temp
c1    format(' gap=',f8.6,' iter=',i3,' for temp=',f8.6)
      return
      end