c Subroutine to calculate the bcs gap function in the pure case for 3He-A 
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-MacLaurin correction term
c input:
c   temp= temperature/T_c
c output:
c   gap= weak-coupling energy gap/(kB T_c)

      subroutine bcsgapa(temp,gap)
      implicit none
      double precision temp,gap,pii2,tm,tm2,tn,help1,help2
      double precision g,gder,eps,sq,y2,y3,atm,at,y,change,tol
      double precision intc,derc
      integer iter,i,m
      parameter(pii2=2.*3.1415926535897932,m=10,tol=1d-9)
c     m = number of term used in the summation (do not use too large value)
c     tol = the required tolerance
      iter=0
      if(temp.ge.1.)then
        gap=0.
        goto 200
      end if
c     initial guess for the gap (y=gap/(2*pi))
      y=2*sqrt(1.-temp)/pii2
      tm=temp*dble(m)
      tm2=tm**2
c     Newton iteration
10    iter=iter+1
      g = 0d0
      gder = 0d0
      y2=y**2
      y3=y**3
c     calculate the sums in G(y) and dG(y)/dy
      do 100 i=m,1,-1
        eps=dble(i)-0.5d0
        tn=temp*eps
        at=atan(y/tn)
        g=g+1d0/eps-0.75d0*temp*(tn*y+(y2-tn**2)*at)/y3
        help1=tn*y*(3*tn**2+y2)/(tn**2+y2)
        help2=(y2-3*tn**2)*at
        gder=gder+0.75d0*temp*(help1+help2)/y**4
100   continue
      atm=atan(y/tm)
      sq=tm2+y2
c     Euler-MacLaurin corrections for G(y)
      help1=tm2*((4d0*m**2-1d0)*temp**2+4*y2)/(16*y2*sq)
      help2=tm*((1d0-4d0*m**2)*temp**2+12*y2)*atm/(16*y3)
      g=g-5d0/6d0-1d0/(24d0*m**2)+help1+help2+log(sqrt(sq)/dble(m))
c     Euler-MacLaurin corrections for G'(y)
      help1=tm2*((3*tm2+5*y2)*temp**2/(sq**2)-12d0)/(16*y3)
      help2=3*tm*((4d0*m**2-1)*temp**2-4*y2)*atm/(16*y**4)
      gder=gder+1d0/y+help1+help2
c     iterate the gap
      change=g/gder
      y=y-change
      if ((abs(change).ge.tol/pii2) .and. (iter.le.100000)) goto 10
      gap=pii2*y
200   continue
c      write(6,1) gap,iter,temp
c1     format(' gap=',f13.10,' iter=',i5,' for temp=',f8.6)
      return
      end