10
0
mirror of https://gitlab.com/scemama/eplf synced 2024-10-31 19:23:55 +01:00
eplf/src/Util.irp.f

137 lines
2.4 KiB
FortranFixed
Raw Normal View History

2009-06-10 00:41:32 +02:00
2009-09-11 17:35:23 +02:00
!recursive double precision function Boys(x,n) result(res)
! implicit none
! include 'constants.F'
!
! real, intent(in) :: x
! integer, intent(in) :: n
!
! ASSERT (x > 0.)
! if (n == 0) then
! res = sqrt(pi/(4.*x))*erf(sqrt(x))
! else
! res = (dble(2*n-1) * Boys(x,(n-1)) - exp(-x) )/(2.*x)
! endif
!
!end function
2009-06-10 00:41:32 +02:00
double precision function binom(k,n)
2011-05-25 11:52:10 +02:00
implicit none
integer, intent(in) :: k,n
double precision, save :: memo(0:20,0:20)
2012-05-31 22:43:06 +02:00
!DEC$ ATTRIBUTES ALIGN : 32 :: memo
2011-05-25 11:52:10 +02:00
integer, save :: ifirst
double precision :: fact
if (ifirst == 0) then
ifirst = 1
integer :: i,j
do j=0,20
2012-05-31 22:43:06 +02:00
do i=0,j
2011-05-25 11:52:10 +02:00
memo(i,j) = fact(j)/(fact(i)*fact(j-i))
enddo
enddo
endif
if (n<20) then
binom = memo(k,n)
else
binom=fact(n)/(fact(k)*fact(n-k))
endif
end
2009-06-10 00:41:32 +02:00
double precision function fact2(n)
implicit none
integer :: n
2011-05-25 11:52:10 +02:00
double precision, save :: memo(100)
2009-06-10 00:41:32 +02:00
integer, save :: memomax = 1
ASSERT (mod(n,2) /= 0)
2011-05-25 11:52:10 +02:00
memo(1) = 1.d0
if (n<2) then
2009-06-10 00:41:32 +02:00
fact2 = 1.d0
2011-05-25 11:52:10 +02:00
else if (n<=memomax) then
2009-06-10 00:41:32 +02:00
fact2 = memo(n)
2011-05-25 11:52:10 +02:00
else
2009-06-10 00:41:32 +02:00
2011-05-25 11:52:10 +02:00
integer :: i
do i=memomax+2,min(n,99),2
memo(i) = memo(i-2)* float(i)
enddo
memomax = min(n,99)
fact2 = memo(memomax)
2009-06-10 00:41:32 +02:00
2011-05-25 11:52:10 +02:00
do i=101,n,2
fact2 *= float(i)
enddo
endif
2009-06-10 00:41:32 +02:00
end function
double precision function fact(n)
implicit none
integer :: n
2011-05-25 11:52:10 +02:00
double precision, save :: memo(0:100)
2009-06-10 00:41:32 +02:00
integer, save :: memomax = 1
2011-05-25 11:52:10 +02:00
memo(0) = 1.d0
memo(1) = 1.d0
2009-06-10 00:41:32 +02:00
if (n<=memomax) then
2011-05-25 11:52:10 +02:00
fact = memo(n)
else
integer :: i
do i=memomax+1,min(n,100)
memo(i) = memo(i-1)*float(i)
enddo
memomax = min(n,100)
fact = memo(memomax)
do i=101,n
fact *= float(i)
enddo
2009-06-10 00:41:32 +02:00
endif
end function
2009-06-15 11:33:16 +02:00
double precision function rintgauss(n)
implicit none
include 'constants.F'
integer :: n
rintgauss = dsqrt_pi
2009-06-15 11:33:16 +02:00
if ( n == 0 ) then
return
else if ( n == 1 ) then
rintgauss = 0.
else if ( mod(n,2) == 1) then
rintgauss = 0.
else
double precision :: fact2
rintgauss = rintgauss/(2.**(n/2))
rintgauss = rintgauss * fact2(n-1)
endif
end function
double precision function goverlap(gamA,gamB,nA)
implicit none
real :: gamA, gamB
integer :: nA(3)
double precision :: gamtot
gamtot = gamA+gamB
goverlap=1.0
integer :: l
double precision :: rintgauss
do l=1,3
goverlap = goverlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
enddo
end function