r9chm.f 1.15 KB
function r9chm (a, c, x)
      real a, c, x, aa(4), bb(4), ct1, ct2, eps,
     1  g1, g2, g3, x2, xn0, xn1, xn2, xn3
      external r1mach
c
      data eps /0.0/
c
      if (eps.eq.0.0) eps = 4.0*r1mach(4)
c
 10   aa(1) = 1.0
      bb(1) = 1.0
c
      x2 = -0.5*x
      bb(2) = 1.0 + (1.0+a)*x2/c
      aa(2) = bb(2) + a*x/c
c
      ct2 = x2/(1.0+c)
      bb(3) = 1.0 + (2.0+bb(2))*(2.0+a)*ct2/3.0
      aa(3) = bb(3) + (1.0+ct2)*a*x/c
c
      do 30 i=3,300
        xn0 = i
        xn1 = i - 1
        xn2 = i - 2
        xn3 = i - 3
        ct1 = 2*i - 3
c
        ct2 = x2/(ct1*(c+xn1))
        g1 = 1.0 + ct2*(xn2-a)
        ct2 = ct2*(a+xn1)/(c+xn2)
        g2 = ct2*(c-xn1+(a+xn0)*x2/(ct1+2.0))
        g3 = ct2*x2*x2*(a+xn2)*(a-xn2)/(ct1*(ct1-2.0)*(c+xn3))
c
        bb(4) = g1*bb(3) + g2*bb(2) + g3*bb(1)
        aa(4) = g1*aa(3) + g2*aa(2) + g3*aa(1)
c
        r9chm = aa(4)/bb(4)
        if (abs(r9chm-aa(1)/bb(1)).lt.eps*abs(r9chm)) go to 40
c
        do 20 j=1,3
          aa(j) = aa(j+1)
          bb(j) = bb(j+1)
 20     continue
 30   continue
      call seteru (40hr9chm   no convergence in 300 iterations, 40,
     1    1, 2)
c
 40   return
      end