c
c
c
      real function rbfrac(a,b,x)
c implementation of Lentz's continued fraction 
c based on ideas from Numerical Recipes (Press et. al., 1986)
      implicit none
      real a,b,x,zero,one,emmax,em,em2,s1,s2,s3,d2m,d2m1,
     $   p2m,p2m1,p2m0,q2m,q2m1,q2m0,fac,frac,fold,tol
      parameter (zero=0.,one=1.,emmax=100d0,tol=3e-7)
      s1=a+b
      s2=a+one
      s3=a-one
      p2m0=one
      q2m0=one
      frac=one
      fac= one-s1*x/s2
      em=zero
c  Abramovitz & Stegun, 1971 (26.5.8) continued fraction
 1      em=em+one
        em2=em+em
c  d2m fraction
        d2m=em*(b-em)*x/((s3+em2)*(a+em2))
        p2m=frac+d2m*p2m0
        q2m=fac+d2m*q2m0
c  d2m1 fraction
        d2m1=-(a+em)*(s1+em)*x/((a+em2)*(s2+em2))
        p2m1=p2m+d2m1*frac
        q2m1=q2m+d2m1*fac
c  d2m and d2m1 fractions combined
        p2m0=p2m/q2m1
        q2m0=q2m/q2m1
        fold=frac
        frac=p2m1/q2m1
        fac=one
        if(em.lt.emmax.and.
     $     abs(fold-frac).gt.tol*abs(frac)) go to 1
      if (em.ge.emmax) then
        pause 
     $  'riembi: A or B too big for EMMAX iterations to converge'
      endif
      rbfrac=frac
      end
c