c Subroutine to calculate the bcs gap function in the pure case (SC or 3He-B)
c Essentially an exact calculation, see (Ajop6,10)
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= the energy gap/T_c
c main program for testing: see file mhsm.f
      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

c laskee matsubara cuf-off energian kun on annettu tarkkuus
c accur, jolla PUHTAAN nesteen gap pitŠŠ olla on oikein.
c input:
c   temp= temperature/Tc
c   accur= allowed error in the gap/Tc
c output:
c   nfreq= number of positive matsubara energies needed to achieve
c          the desired accuracy (assuming that all higher energies are
c          simply omitted from the sum in the gap equation)
      subroutine cutoff(temp,accur,nfreq)
      implicit none
      integer nfreq
      double precision temp, gap, fprime, err
      double precision pi, zeta37, eps, sqterm
      double precision acc, accur,small
      parameter (pi = 3.141592654,small=.000001)
      parameter (zeta37 = 1.2020569*7.)
      if(accur.ge.0.5)then
        nfreq=int(accur+small)
      else
        acc = accur/(2*pi*temp)
        gap = 1.7638*sqrt(1. - temp)/(2*pi*temp)
        fprime = 0.
        err = zeta37
        nfreq = 0
80      nfreq = nfreq + 1
        eps = float(nfreq) - 0.5
        err = err - 1./eps**3
        sqterm = sqrt(eps**2 + gap**2)
c seuraava 2:nen tulee koska virhe itse asiassa on gap**2*err/2
        fprime = fprime + 2./sqterm**3
        if (gap*err/fprime .gt. acc) goto 80
      end if
      return
      end