Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/GKCoulomb/Tljpmf/src/T5_mod.f90~
3c0a5500b5baimproved_sugama

Last Change
Raw File

T5_mod.f90~
View Options

MODULE T5
!
! contains routines to compute T5 coefficients
!
USE FMZM
USE PREC_CONST
!
IMPLICIT NONE
PRIVATE
INTEGER, PARAMETER :: cp = dp !real(dp), enforced through the code
interface CBinomial
module procedure CBinomial_cp
module procedure CBinomial_i
end interface CBinomial
public :: T5zero, Bruteforce_T5
CONTAINS
integer function T5zero( l, k, m, p, j ) result( return_value )
implicit none
integer, intent(in) :: l, k, m, p, j
integer :: off1, off2, D
return_value = 0
off1 = mod(l-m,2)
off2 = mod(off1+1,2)
D = (l-m)/2 ! Floor of l-m, Delta in script
if( mod(l+m+p,2) /= 0 .or. &
& D+k<0 .or. &x
& 2*D+k+(1-off2-p)/2<0 .or. &
& 2*off2+2*k+l-m-p-2*j<0 .or. &
& (m==0 .and. T4zero(l,k,p,j)) .or. &
& .not. (l-m>=max(0,p+2*j-2*k) .and. mod(l-m+p,2)==0) ) then
return_value = 1
end if
end function T5zero
logical function T4zero( l, k, p, j ) result( return_value )
implicit none
integer, intent(in) :: l, k, p, j
integer :: tmp
return_value = .false.
tmp = 1 - modulo(l-p,2)
if( tmp == 0 .or. l+2*k/=p+2*j ) return_value = .true.
end function T4zero
type(FM) function Bruteforce_T5( l, k, m, p, j ) result( return_value )
use FMZM
implicit none
integer, intent(in) :: l, k, m, p, j
real(cp) :: rl, rk, rm, rp, rj
type(FM), save :: fmPI, fmsqrtPI
type(FM), save :: num(10), den(10)
integer :: off1, off2, r=0, i=0, s=0, eta=0, D=0
real(cp) :: rr, ri, rs, reta, rD
CALL FM_ENTER_USER_FUNCTION( return_value )
num = TO_FM(1)
den = TO_FM(1)
return_value = TO_FM(0)
off1 = mod(l-m,2)
off2 = mod(off1+1,2)
D = (l-m)/2 ! Floor of l-m, Delta in script
rD= real(D,kind=cp)
rl = real(l, kind=cp)
rk = real(k, kind=cp)
rm = real(m, kind=cp)
rp = real(p, kind=cp)
rj = real(j, kind=cp)
call FM_PI(fmPI)
fmsqrtPI = sqrt(fmPI)
num(1) = TO_FM(1 - 2*mod(m+j,2)) * TO_FM(2)**l
den(1) = CFactorial( rp )
!do r = max(0,p-off1-3*D+2*j-k), k
do r = max(0,p-off1-3*D+2*j-k,j-D+(p-off2-1)/2), k
rr = real( r, kind=cp )
num(2) = TO_FM(1 - 2*mod(r,2))
num(3) = CBinomial( rl+rk+.5_cp, rk-rr )
den(2) = CFactorial( rr )
! NO ETA SUMS
do i = max(0,(p-off1)/2-r-D+j), min(D,r+2*D-j-(p-off1)/2)
if(r<p-off1-3*D+2*j-k) exit
ri = real(i,kind=cp)
num(4) = CBinomial( l, 2*i + m + off1 )
num(5) = CBinomial( ri + .5_cp*(rl+rm-real(off2,kind=cp)), rl )
num(6) = CFactorial( 2._cp*ri + rm + real(off1,kind=cp) )
den(3) = CFactorial( 2._cp*ri + real(off1,kind=cp) )
do s = max(0,(p-off1)/2-i), r+D-i-j
rs = real(s,kind=cp)
num(7) = CBinomial( r+D-i, s )
num(8) = CBinomial( r+D-i-s, j )
num(9) = CFactorial( real(off1,kind=cp) + 2._cp*(ri+rs) )
num(10)= CFactorial( rr+rD-ri-rs )
den(4) = TO_FM(2)**(off1 + 2*(i+s))
den(5) = CFactorial( ri + rs - .5_cp*(rp-real(off1,kind=cp)) )
! IIF off1+2*(i+s) >= p
! IIF r+D-i-s >= j
!if( r+D-i-s >= j .and. off1 + 2*(i+s) >= p ) then
return_value = return_value + product(num)/product(den)
!end if
end do !s
end do !i
!ETA SUMS
if( r < j-D+(p-1-off2)/2 ) cycle
do i = 0, D-off2
ri = real(i,kind=cp)
num(4) = fmsqrtPI * CBinomial( l, 2*i+m+off2 ) / TO_FM(2)
num(5) = CBinomial( ri + .5_cp*(rl+rm-real(off1,kind=cp)), rl )
num(6) = CFactorial( 2._cp*ri + rm + real(off2,kind=cp) )
den(3) = CFactorial( 2._cp*ri + real(off2,kind=cp) )
do s = max(0,(p-1-off2)/2-i), min(r+D-i,2*(r+D)-j-i+(1+off2-p)/2)
rs = real(s,kind=cp)
num(7) = CBinomial( r+D-i, s )
do eta = max(0,j+i+s-r-D), min(D+r+(1-off2-p)/2,i+s+(1+off2-p)/2)
reta = real(eta,kind=cp)
num(8) = CBinomial( r+D+eta-i-s, j )
num(9) = CFactorial( real(1+off2,kind=cp) + 2._cp*(ri+rs-reta) )
num(10)= CFactorial( rr+rD+reta-ri-rs )
den(4) = CFactorial( .5_cp - reta )
den(5) = TO_FM(2)**(1+off2+2*(i+s-eta))
den(6) = CFactorial( ri+rs-reta-.5_cp*rp + .5_cp*real(1+off2,kind=cp) )
! IIF r+D+eta-i-s >= j
!if( r+D+eta-i-s >= j ) then
return_value = return_value + product(num)/product(den)
!end if
end do !eta
end do !s
end do !i
den(6) = TO_FM(1)
end do !r
call FM_EXIT_USER_FUNCTION( return_value )
end function Bruteforce_T5
type(FM) function CFactorial( val ) result( return_value )
use FMZM
implicit none
integer :: ival, iloop, istr
integer, parameter :: LOWFACT=-200, HIGHFACT=400
type(FM), save :: fval, bval, rval, step
real(cp), intent(in) :: val
type(FM), save :: PI, sqrtPI
type(FM), save :: store(LOWFACT:HIGHFACT)
logical, save :: is_init(LOWFACT:HIGHFACT) = .false., isPI = .false.
call FM_ENTER_USER_FUNCTION(return_value)
if( .not. isPI ) then
call FM_PI(PI)
sqrtPI = sqrt(PI)
isPI = .true.
end if
return_value = TO_FM(0)
bval = TO_FM(1)
fval = TO_FM(val)
step = sign(bval,fval)
rval = step
istr = nint(val*2._cp)
! If we're dealing with half-integers, recursion relation with sqrt(PI)
if( modulo(val,1._cp) > 0.25_cp ) then
bval = sqrtPI
rval = rval / TO_FM(2)
end if
if( istr<0 ) istr = (istr-1)/2
if( istr < LOWFACT .or. istr > HIGHFACT ) then
!! print *, "Out of range factorial."
!! return_value = GAMMA( val + 1._cp )
return_value = FACTORIAL(fval)
call FM_EXIT_USER_FUNCTION( return_value )
return
end if
if(.not. is_init(istr)) then
store(istr) = TO_FM(1)
ival = abs(nint(val+0.25_cp))
do iloop = 1, ival
store(istr) = store(istr) * rval
rval = rval + step
end do
if( val < 0._cp ) store(istr) = TO_FM(1) / store(istr)
store(istr) = store(istr) * bval
is_init(istr) = .true.
end if
return_value = store(istr)
call FM_EXIT_USER_FUNCTION(return_value)
end function CFactorial
type(FM) function CBinomial_cp( n, k ) result( return_value )
use FMZM
implicit none
real(cp), intent(in) :: n, k
type(FM), save :: f1, f2, f3
! See http://mathworld.wolfram.com/BinomialCoefficient.html
call FM_ENTER_USER_FUNCTION(return_value)
return_value = TO_FM(0)
if( mod(n,1._cp) < .1_cp .and. mod(k,1._cp) < .1_cp ) then
return_value = CBinomial_i( nint(n), nint(k) )
else
f1 = CFactorial( n )
f2 = CFactorial( k )
f3 = CFactorial( n-k )
return_value = f1 / (f2*f3)
end if
call FM_EXIT_USER_FUNCTION(return_value)
end function CBinomial_cp
type(FM) function CBinomial_i( n, k ) result( return_value )
use FMZM
implicit none
integer, intent(in) :: n, k
real(cp) :: rn, rk
type(FM), save :: f1, f2, f3
! See http://mathworld.wolfram.com/BinomialCoefficient.html
call FM_ENTER_USER_FUNCTION(return_value)
return_value = TO_FM(0)
rn = real(n, kind=cp)
rk = real(k, kind=cp)
if( k>=0 .and. n>=k ) then
f1 = CFactorial(rn)
f2 = CFactorial(rk)
f3 = CFactorial(rn-rk)
return_value = f1 / (f2*f3)
else if( n<0 ) then
if( k>=0 ) then
f1 = CFactorial( -rn+rk-1._cp )
f2 = CFactorial( rk )
f3 = CFactorial( -rn-1._cp )
return_value = TO_FM(-1)**k * f1 / ( f2*f3 )
else if( k<=n ) then
f1 = CFactorial( -rk-1._cp )
f2 = CFactorial( rn-rk )
f3 = CFactorial( -rn-1._cp )
return_value = TO_FM(-1)**(n-k) * f1 / (f2*f3)
end if
end if
call FM_EXIT_USER_FUNCTION(return_value)
end function CBinomial_i
END MODULE T5
c4science · Help