Blame view
fvn_fnlib/r9chm.f
1.15 KB
38581db0c git-svn-id: https... |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 |
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 |