9 real,
dimension(0:20) :: fac(0:20)
11 End Module rdf_variables
29 SUBROUTINE rdf_int ( dense, mseg, radius, rdf )
35 REAL,
INTENT(IN) :: dense
36 REAL,
INTENT(IN) :: mseg
37 REAL,
INTENT(IN) :: radius
38 REAL,
INTENT(OUT) :: rdf
42 REAL :: rg, omega, radind, cfun, facul
53 omega = ( mseg - 1.0 ) / mseg
58 radind = radius - 1.0 -
REAL(n)
59 rg = rg + (-12.0*dense)**n *( (1.0+dense/2.0-omega*(1.0-dense)) &
60 *cfun(dense,omega,2,n,(n+1),radind) &
61 +(1.0+dense*2.0-(1.0+dense/2.0)*omega +(1.0-dense)/2.0*omega**2) &
62 *cfun(dense,omega,1,n,(n+1),radind) &
63 +((1.0+2.0*dense)*omega-(1.0-dense)*omega**2 ) &
64 *cfun(dense,omega,0,n,(n+1),radind) &
65 -((1.0+dense/2.0)*omega-(1.0-dense)/2.0*omega**2 ) &
66 *cfun(dense,omega,1,n,(n+1),(radind-1.0)) &
67 -((1.0+2.0*dense)*omega-(1.0-dense)*omega**2 ) &
68 *cfun(dense,omega,0,n,(n+1),(radind-1.0)) )
75 END SUBROUTINE rdf_int
82 REAL FUNCTION lnfac(xx)
87 REAL,
INTENT(IN) :: xx
94 DATA cof,stp/76.18009172947146,-86.50532032941677, &
95 24.01409824083091,-1.231739572450155,.1208650973866179e-2, &
96 -.5395239384953e-5,2.5066282746310005/
102 tmp = ( x + 0.5 ) * log( tmp ) - tmp
103 ser = 1.000000000190015
106 ser = ser + cof(j) / y
108 lnfac = tmp + log( stp * ser / x )
118 REAL FUNCTION facul(n)
123 INTEGER,
INTENT(IN) :: n
127 INTEGER,
SAVE :: ntop
128 REAL,
SAVE :: table(165)
131 DATA ntop,table(1) / 0, 1.0 /
134 write (*,*)
' Negative factorial in FUNCTION fac' 136 ELSE IF ( n <= ntop )
THEN 138 ELSE IF ( n <= 165 )
THEN 141 table(j+1)=mult*table(j)
146 facul = log( lnfac(
REAL(n)+1.0 ) )
157 COMPLEX FUNCTION anroot(dens,omeg,m)
162 REAL,
INTENT(IN) :: dens
163 REAL,
INTENT(IN OUT) :: omeg
164 INTEGER,
INTENT(IN) :: m
167 REAL :: coeff(3),qq,rr,aa,bb,ohm, pi
169 COMPLEX,
SAVE :: anr0,anr1,anr2
171 DATA roh2,anr0,anr1,anr2 / 0.0, (0.0, 0.0), (0.0, 0.0), (0.0, 0.0) /
176 IF (dens /= roh2)
THEN 178 coeff(1)=(6.0*dens-(1.0-dens)*omeg) /(1.0-dens)
179 coeff(2)=(18.0*dens**2 - 6.0*dens*(1.0-dens)*omeg) /((1.0-dens)**2)
180 coeff(3)=(-12.0*dens*(1.0+2.0*dens) +12.0*dens*(1.0-dens)*omeg) /((1.0-dens)**2)
181 qq = ( coeff(1)**2 - 3.0*coeff(2) ) / 9.0
182 rr = ( 2.0*coeff(1)**3 - 9.0*coeff(1)*coeff(2) + 27.0*coeff(3) ) / 54.0
184 IF ( rr**2 < qq**3 )
THEN 186 ohm = acos( rr/qq**(3.0/2.0) )
187 anr0=-2.0* sqrt(qq)* cos( ohm/3.0 )-coeff(1)/3.0
188 anr1=-2.0* sqrt(qq)* cos((ohm+2.0*pi)/3.0)-coeff(1)/3.0
189 anr2=-2.0* sqrt(qq)* cos((ohm-2.0*pi)/3.0)-coeff(1)/3.0
193 aa = - ( abs(rr) + sqrt( rr**2 -qq**3 ) )**(1.0/3.0)
194 aa = rr/ abs(rr) * aa
195 IF ( aa /= 0.0 )
THEN 201 anr0 = cmplx( aa+bb -coeff(1)/3.0 , 0.0 )
202 anr1 = cmplx((-(aa+bb)/2.0 -coeff(1)/3.0), (sqrt(3.0)/2.0*(aa-bb)))
203 anr2 = cmplx((-(aa+bb)/2.0 -coeff(1)/3.0), -(sqrt(3.0)/2.0*(aa-bb)))
211 ELSE IF (m == 1)
THEN 213 ELSE IF (m == 2)
THEN 216 WRITE (*,*)
' Only 3 roots do exist. Asked for: ',m+1
248 COMPLEX FUNCTION afun (dense,om,an1,an2,ak,aalpha)
254 INTEGER :: i,j,an1,an2,ak
256 COMPLEX :: aalpha, af
263 IF ((i+j) <= an2 .AND. (i+j+j) >= (ak-an1))
THEN 265 af = af + fac(an2)*fac(i+an1+j+j) &
266 /(fac(i)*fac(j)*fac(an2-i-j)*fac(i+j+j+an1-ak)) &
267 * (1.0+dense/2.0-(1.0-dense)/2.0*om)**
REAL(i) &
268 * (1.0+2.0*dense-(1.0-dense)*om)**
REAL(an2-i-j) &
269 * ((1.0-dense)**2 *om/(12.0*dense))**
REAL(j) &
270 * aalpha**
REAL(an1+i+j+j-ak)
286 COMPLEX FUNCTION bfun (dense,om,bn1,bn2,bn3,bi,bal)
292 INTEGER :: k,l,bn1,bn2,bn3,bi,bal,beta,gamma
293 REAL :: dense, om, imsum,factor
294 COMPLEX :: talpha, anroot, afun, bf,comfra
300 ELSE IF (bal == 1)
THEN 303 ELSE IF (bal == 2)
THEN 307 WRITE(*,*)
' no proper value for bal in BFUN !!' 317 talpha = anroot(dense,om,bal)
327 factor=(-1.0)**
REAL(bn3-bi-k)*fac(bn3-1+l)*fac(bn3+bn3-1-bi-k-l) &
328 /( fac(k)*fac(l)*fac(bn3-bi-k-l)*fac(bn3-1)*fac(bn3-1) )
329 comfra= afun(dense,om,bn1,bn2,k,talpha) /( (talpha &
330 - anroot(dense,om,beta))**
REAL(bn3+l) *(talpha &
331 - anroot(dense,om,gamma))**
REAL(bn3+bn3-bi-k-l) )
332 imsum = imsum + factor*aimag(comfra)
333 bf = bf + factor * comfra
338 bfun = 1.0 / (1.0-dense)**
REAL(bn3+bn3) * bf
348 REAL FUNCTION cfun (dense,om,cn1,cn2,cn3,rad)
354 INTEGER :: cn1,cn2,cn3, i, alph
355 REAL :: dense,om,rad,fact
356 COMPLEX :: anroot,bfun,cf,root,swing,compl
359 cf = cmplx( 0.0, 0.0 )
360 IF (abs(rad) < 1.e-10) rad = 0.0
365 root = anroot(dense,om,alph)
367 compl =bfun(dense,om,cn1,cn2,cn3,i,alph)*cexp(swing)
368 fact = rad**
REAL(i-1) / fac(i-1)
369 cf = cf + compl * fact
375 IF ( (cn1+cn2) < (cn3+cn3+cn3) )
THEN 377 ELSE IF ( (cn1+cn2) == (cn3+cn3+cn3) )
THEN 379 cfun = cf + (1.0+dense/2.0)**
REAL(cn2) /((1.0-dense)**
REAL(cn3+cn3))
380 WRITE (*,*)
'in Schleife',cfun
385 WRITE(*,*)
' no proper case in CFUN !!'