57 lines
1.8 KiB
Fortran
57 lines
1.8 KiB
Fortran
FUNCTION rj(x,y,z,p)
|
|
REAL rj,p,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6,C7,C8
|
|
PARAMETER (ERRTOL=.05,TINY=2.5e-13,BIG=9.E11,C1=3./14.,C2=1./3.,
|
|
*C3=3./22.,C4=3./26.,C5=.75*C3,C6=1.5*C4,C7=.5*C2,C8=C3+C3)
|
|
CU USES rc,rf
|
|
REAL a,alamb,alpha,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee,
|
|
*fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt,rc,rf
|
|
if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z,abs(p)).lt.TINY.or.max(x,y,
|
|
*z,abs(p)).gt.BIG)pause 'invalid arguments in rj'
|
|
sum=0.
|
|
fac=1.
|
|
if(p.gt.0.)then
|
|
xt=x
|
|
yt=y
|
|
zt=z
|
|
pt=p
|
|
else
|
|
xt=min(x,y,z)
|
|
zt=max(x,y,z)
|
|
yt=x+y+z-xt-zt
|
|
a=1./(yt-p)
|
|
b=a*(zt-yt)*(yt-xt)
|
|
pt=yt+b
|
|
rho=xt*zt/yt
|
|
tau=p*pt/yt
|
|
rcx=rc(rho,tau)
|
|
endif
|
|
1 continue
|
|
sqrtx=sqrt(xt)
|
|
sqrty=sqrt(yt)
|
|
sqrtz=sqrt(zt)
|
|
alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
|
|
alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)**2
|
|
beta=pt*(pt+alamb)**2
|
|
sum=sum+fac*rc(alpha,beta)
|
|
fac=.25*fac
|
|
xt=.25*(xt+alamb)
|
|
yt=.25*(yt+alamb)
|
|
zt=.25*(zt+alamb)
|
|
pt=.25*(pt+alamb)
|
|
ave=.2*(xt+yt+zt+pt+pt)
|
|
delx=(ave-xt)/ave
|
|
dely=(ave-yt)/ave
|
|
delz=(ave-zt)/ave
|
|
delp=(ave-pt)/ave
|
|
if(max(abs(delx),abs(dely),abs(delz),abs(delp)).gt.ERRTOL)goto 1
|
|
ea=delx*(dely+delz)+dely*delz
|
|
eb=delx*dely*delz
|
|
ec=delp**2
|
|
ed=ea-3.*ec
|
|
ee=eb+2.*delp*(ea-ec)
|
|
rj=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))+
|
|
*delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave))
|
|
if (p.le.0.) rj=a*(b*rj+3.*(rcx-rf(xt,yt,zt)))
|
|
return
|
|
END
|