30 lines
914 B
Fortran
30 lines
914 B
Fortran
FUNCTION rc(x,y)
|
|
REAL rc,x,y,ERRTOL,TINY,SQRTNY,BIG,TNBG,COMP1,COMP2,THIRD,C1,C2,
|
|
*C3,C4
|
|
PARAMETER (ERRTOL=.04,TINY=1.69e-38,SQRTNY=1.3e-19,BIG=3.E37,
|
|
*TNBG=TINY*BIG,COMP1=2.236/SQRTNY,COMP2=TNBG*TNBG/25.,THIRD=1./3.,
|
|
*C1=.3,C2=1./7.,C3=.375,C4=9./22.)
|
|
REAL alamb,ave,s,w,xt,yt
|
|
if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.TINY.or.(x+
|
|
*abs(y)).gt.BIG.or.(y.lt.-COMP1.and.x.gt.0..and.x.lt.COMP2))pause
|
|
*'invalid arguments in rc'
|
|
if(y.gt.0.)then
|
|
xt=x
|
|
yt=y
|
|
w=1.
|
|
else
|
|
xt=x-y
|
|
yt=-y
|
|
w=sqrt(x)/sqrt(xt)
|
|
endif
|
|
1 continue
|
|
alamb=2.*sqrt(xt)*sqrt(yt)+yt
|
|
xt=.25*(xt+alamb)
|
|
yt=.25*(yt+alamb)
|
|
ave=THIRD*(xt+yt+yt)
|
|
s=(yt-ave)/ave
|
|
if(abs(s).gt.ERRTOL)goto 1
|
|
rc=w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave)
|
|
return
|
|
END
|