-
-
Notifications
You must be signed in to change notification settings - Fork 3
/
libboys.f90
114 lines (106 loc) · 5.3 KB
/
libboys.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
!
! libboys - a FORTRAN library to numerically calculate the Boys function
! Copyright (C) 2014-2016 Michael Böhme <boehme.mic@gmail.com>
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MDRCHANTABILITY or FITNDSS FOR A PARTICULAR PURPOSD. See the GNU
! Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this program; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
!
! ***************************************************************************************
! * Boys function F_n(x)
! ***************************************************************************************
function Boys_func(n, x) result(res)
use libboys_data
implicit none
!
integer, intent(in) :: n
double precision, intent(in) :: x
double precision :: res, sfac, dx, dxi
integer :: i, j
double precision, dimension(6) :: n_fac_dble = (/ 1d0, 2d0, 6d0, 24d0, 120d0, 720d0 /)
double precision, dimension(31) :: n_fac2_dble = (/ &
1d0, 1d0, 1d0, 2d0, 3d0, 8d0, 15d0, 48d0, 105d0, 384d0, 945d0, 3840d0, &
10395d0,46080d0,135135d0,645120d0,2027025d0,10321920d0,34459425d0,185794560d0,654729075d0,3715891200d0,13749310575d0,81749606400d0, &
316234143225d0,1961990553600d0,7905853580625d0,51011754393600d0,213458046676875d0,1428329123020800d0,6190283353629375d0 /)
double precision, parameter :: Pi = 3.1415926535897932d0
double precision, parameter :: eps = 1d-14
integer, parameter :: MAX_RECURSION = 6
double precision :: epsrel
!
res = 0d0
if ( n .eq. 0 ) then
if ( x .lt. eps ) then
res = 1d0
else
res = dsqrt( Pi / (4d0*x) ) * derf(dsqrt(x))
end if
else
if ( n .gt. libBoysMaxN ) then
write(*,*) "libboys error: not implemented!"
return
end if
if ( n .lt. 0 ) then
write(*,*) "libboys error: n < 0!"
return
end if
if ( x .lt. eps ) then
res = 1d0 / ( 2d0*dble(n) + 1d0 )
else if ( x .gt. 50d0 ) then
res = n_fac2_dble(2*n-1 +2) / 2d0**(n+1) * dsqrt(Pi/x**(2*n+1))
else
if ( x .ge. 10d0 ) then
j = int((x-9.95d0)*10d0) + 1
dx = BoysFuncValuesL(j, 1) - x
dxi = dx
res = BoysFuncValuesL(j, n + 2)
epsrel = res * eps
do i = 1, MAX_RECURSION
sfac = BoysFuncValuesL(j, n + 2 + i) * dxi / n_fac_dble(i)
res = res + sfac
if ( abs(sfac) .lt. epsrel ) then
return
end if
dxi = dxi * dx
end do
else if ( x .ge. 5d0 ) then
j = int((x-4.975d0)*20d0) + 1
dx = BoysFuncValuesM(j, 1) - x
dxi = dx
res = BoysFuncValuesM(j, n + 2)
epsrel = res * eps
do i = 1, MAX_RECURSION
sfac = BoysFuncValuesM(j, n + 2 + i) * dxi / n_fac_dble(i)
res = res + sfac
if ( abs(sfac) .lt. epsrel ) then
return
end if
dxi = dxi * dx
end do
else
j = int(x*40d0+0.5d0) + 1
dx = BoysFuncValuesS(j, 1) - x
dxi = dx
res = BoysFuncValuesS(j, n + 2)
epsrel = res * eps
do i = 1, MAX_RECURSION
sfac = BoysFuncValuesS(j, n + 2 + i) * dxi / n_fac_dble(i)
res = res + sfac
if ( abs(sfac) .lt. epsrel ) then
return
end if
dxi = dxi * dx
end do
end if
end if
end if
end function