c subroutine to calculate the Cross coefficients for 3He-A
c This calculation is 6 decimal correct at any temperature
c input:
c   temp= temperature/T_c
c output:
c   gap= the maximum A phase energy gap/ (kB T_c)
c   alpha,beta,gamma = the Cross coefficients, see for example
c      M.C. Cross, J. Low Temp. Phys. {\bf 21}, 525 (1975)
c for method of calculation see
c   http://boojum.hut.fi/research/theory/qc/bcsgap.html
c subprograms:
c   bcsgapa

      subroutine crossc(temp,gap,alpha,beta,gamma)
      implicit none
      double precision temp,gap,alpha,beta,gamma,intc,derc,eps
      double precision pii2,y,y2,y3,tm,tm2,tn,at,atm,sq,help1,help2
      integer i,m
      parameter(pii2=2.*3.1415926535897932,m=10)

      call bcsgapa(temp,gap)
      y=gap/pii2
      y2=y**2
      y3=y**3
      alpha=0d0
      beta=0d0
      gamma=0d0
      do 100 i=m,1,-1
        eps=dble(i)-0.5d0
        tn=temp*eps
        at=atan(y/tn)
        sq=y**2+tn**2
        alpha=alpha+(tn*y*(3*tn**2+y2)+(y**4-2*y2*tn**2-3*tn**4)*at)/sq
        beta=beta+(3*tn**2+y2)*at-3*tn*y
        gamma=gamma+3d0/eps+4.5d0*temp*(tn*y-sq*at)/y3
100   continue
      alpha=3*temp*alpha/(16*y3)
      beta=3*temp*beta/(4*y3)
      tm=temp*dble(m)
      tm2=tm**2
      sq=y2+tm2
      atm=atan(y/tm)
c     Euler-MacLaurin corrections
      help1=((3*tm2+5*y2)*temp**2/sq**2-12d0)*tm2/(64*y2)
      help2=3*tm*((4d0*m**2-1d0)*temp**2-4*y2)*atm/(64*y3)
      alpha=alpha+0.25d0+help1+help2

      help1=((4d0*m**2-1)*temp**2+4*y2)*(3*y*tm2+2*y3-3*tm*sq*atm)
      beta=beta+help1/(16*sq*y3)

      help1=0.375*tm*((4*m**2-1d0)*temp**2+12*y2)*atm/y3-1d0/(8*m**2)
      help2=0.375*(1d0-4d0*m**2)*temp**2/y2+1.5*log(sq/tm2)-4d0
      gamma=gamma+help1+help2

c      write(*,*) 'alpha = ',alpha,' beta = ',beta,' gamma = ',gamma
      return
      end