FUNCTION rd(x,y,z) REAL rd,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6 PARAMETER (ERRTOL=.05,TINY=1.e-25,BIG=4.5E21,C1=3./14.,C2=1./6., *C3=9./22.,C4=3./26.,C5=.25*C3,C6=1.5*C4) REAL alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty, *sqrtz,sum,xt,yt,zt if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y, *z).gt.BIG)pause 'invalid arguments in rd' xt=x yt=y zt=z sum=0. fac=1. 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt) sqrtz=sqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz sum=sum+fac/(sqrtz*(zt+alamb)) fac=.25*fac xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=.2*(xt+yt+3.*zt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1 ea=delx*dely eb=delz*delz ec=ea-eb ed=ea-6.*eb ee=ed+ec+ec rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3* *ec+delz*C4*ea)))/(ave*sqrt(ave)) return END C (C) Copr. 1986-92 Numerical Recipes Software .