Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates Bessel functions
10 : !> \note
11 : !> Functions adapted from netlib
12 : !> \par History
13 : !> March-2006: Bessel Transform (JGH)
14 : !> \author JGH (10-02-2001)
15 : ! **************************************************************************************************
16 : MODULE bessel_lib
17 :
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: fac,&
20 : pio2
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 : PUBLIC :: bessk0, bessk1, bessel0
27 :
28 : CONTAINS
29 :
30 : ! **************************************************************************************************
31 : !> \brief ...
32 : !> \param x must be positive
33 : !> \return ...
34 : ! **************************************************************************************************
35 187200 : ELEMENTAL FUNCTION bessk0(x)
36 :
37 : REAL(KIND=dp), INTENT(IN) :: x
38 : REAL(KIND=dp) :: bessk0
39 :
40 : REAL(KIND=dp), PARAMETER :: p1 = -0.57721566_dp, p2 = 0.42278420_dp, p3 = 0.23069756_dp, &
41 : p4 = 0.3488590e-1_dp, p5 = 0.262698e-2_dp, p6 = 0.10750e-3_dp, p7 = 0.74e-5_dp, &
42 : q1 = 1.25331414_dp, q2 = -0.7832358e-1_dp, q3 = 0.2189568e-1_dp, q4 = -0.1062446e-1_dp, &
43 : q5 = 0.587872e-2_dp, q6 = -0.251540e-2_dp, q7 = 0.53208e-3_dp
44 :
45 : REAL(KIND=dp) :: y
46 :
47 187200 : IF (x < 2.0_dp) THEN
48 0 : y = x*x/4.0_dp
49 : bessk0 = (-LOG(x/2.0_dp)*bessi0(x)) + (p1 + y* &
50 0 : (p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))))
51 : ELSE
52 187200 : y = (2.0_dp/x)
53 : bessk0 = (EXP(-x)/SQRT(x))*(q1 + y*(q2 + y* &
54 187200 : (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
55 : END IF
56 :
57 187200 : END FUNCTION bessk0
58 :
59 : ! **************************************************************************************************
60 : !> \brief ...
61 : !> \param x must be positive
62 : !> \return ...
63 : ! **************************************************************************************************
64 187200 : ELEMENTAL FUNCTION bessk1(x)
65 :
66 : REAL(KIND=dp), INTENT(IN) :: x
67 : REAL(KIND=dp) :: bessk1
68 :
69 : REAL(KIND=dp), PARAMETER :: p1 = 1.0_dp, p2 = 0.15443144_dp, p3 = -0.67278579_dp, &
70 : p4 = -0.18156897_dp, p5 = -0.1919402e-1_dp, p6 = -0.110404e-2_dp, p7 = -0.4686e-4_dp, &
71 : q1 = 1.25331414_dp, q2 = 0.23498619_dp, q3 = -0.3655620e-1_dp, q4 = 0.1504268e-1_dp, &
72 : q5 = -0.780353e-2_dp, q6 = 0.325614e-2_dp, q7 = -0.68245e-3_dp
73 :
74 : REAL(KIND=dp) :: y
75 :
76 187200 : IF (x < 2.0_dp) THEN
77 0 : y = x*x/4.0_dp
78 : bessk1 = (LOG(x/2.0_dp)*bessi1(x)) + (1.0_dp/x)* &
79 : (p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y* &
80 0 : (p6 + y*p7))))))
81 : ELSE
82 187200 : y = 2.0_dp/x
83 : bessk1 = (EXP(-x)/SQRT(x))*(q1 + y*(q2 + y* &
84 187200 : (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
85 : END IF
86 :
87 187200 : END FUNCTION bessk1
88 :
89 : ! **************************************************************************************************
90 : !> \brief ...
91 : !> \param x ...
92 : !> \return ...
93 : ! **************************************************************************************************
94 0 : ELEMENTAL FUNCTION bessi0(x)
95 :
96 : REAL(KIND=dp), INTENT(IN) :: x
97 : REAL(KIND=dp) :: bessi0
98 :
99 : REAL(KIND=dp), PARAMETER :: p1 = 1.0_dp, p2 = 3.5156229_dp, p3 = 3.0899424_dp, &
100 : p4 = 1.2067492_dp, p5 = 0.2659732_dp, p6 = 0.360768e-1_dp, p7 = 0.45813e-2_dp, &
101 : q1 = 0.39894228_dp, q2 = 0.1328592e-1_dp, q3 = 0.225319e-2_dp, q4 = -0.157565e-2_dp, &
102 : q5 = 0.916281e-2_dp, q6 = -0.2057706e-1_dp, q7 = 0.2635537e-1_dp, q8 = -0.1647633e-1_dp, &
103 : q9 = 0.392377e-2_dp
104 :
105 : REAL(KIND=dp) :: ax, y
106 :
107 0 : IF (ABS(x) < 3.75_dp) THEN
108 0 : y = (x/3.75_dp)**2
109 : bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
110 0 : (p5 + y*(p6 + y*p7)))))
111 : ELSE
112 0 : ax = ABS(x)
113 0 : y = 3.75_dp/ax
114 : bessi0 = (EXP(ax)/SQRT(ax))*(q1 + y*(q2 + y* &
115 : (q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y* &
116 0 : (q8 + y*q9))))))))
117 : END IF
118 :
119 0 : END FUNCTION bessi0
120 :
121 : ! **************************************************************************************************
122 : !> \brief ...
123 : !> \param x ...
124 : !> \return ...
125 : ! **************************************************************************************************
126 0 : ELEMENTAL FUNCTION bessi1(x)
127 :
128 : REAL(KIND=dp), INTENT(IN) :: x
129 : REAL(KIND=dp) :: bessi1
130 :
131 : REAL(KIND=dp), PARAMETER :: p1 = 0.5_dp, p2 = 0.87890594_dp, p3 = 0.51498869_dp, &
132 : p4 = 0.15084934e0_dp, p5 = 0.2658733e-1_dp, p6 = 0.301532e-2_dp, p7 = 0.32411e-3_dp, &
133 : q1 = 0.39894228_dp, q2 = -0.3988024e-1_dp, q3 = -0.362018e-2_dp, q4 = 0.163801e-2_dp, &
134 : q5 = -0.1031555e-1_dp, q6 = 0.2282967e-1_dp, q7 = -0.2895312e-1_dp, q8 = 0.1787654e-1_dp, &
135 : q9 = -0.420059e-2_dp
136 :
137 : REAL(KIND=dp) :: ax, y
138 :
139 0 : IF (ABS(x) < 3.75_dp) THEN
140 0 : y = (x/3.75_dp)**2
141 : bessi1 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
142 0 : (p5 + y*(p6 + y*p7)))))
143 : ELSE
144 0 : ax = ABS(x)
145 0 : y = 3.75_dp/ax
146 : bessi1 = (EXP(ax)/SQRT(ax))*(q1 + y*(q2 + y* &
147 : (q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y* &
148 0 : (q8 + y*q9))))))))
149 0 : IF (x < 0.0_dp) bessi1 = -bessi1
150 : END IF
151 :
152 0 : END FUNCTION bessi1
153 :
154 : ! **************************************************************************************************
155 : !> \brief ...
156 : !> \param x ...
157 : !> \param l ...
158 : !> \return ...
159 : ! **************************************************************************************************
160 801600 : ELEMENTAL IMPURE FUNCTION bessel0(x, l)
161 : !
162 : ! Calculates spherical Bessel functions
163 : ! Abramowitz & Stegun using Formulas 10.1.2, 10.1.8, 10.1.9
164 : ! Adapted from P. Bloechl
165 : !
166 : REAL(KIND=dp), INTENT(IN) :: x
167 : INTEGER, INTENT(IN) :: l
168 : REAL(KIND=dp) :: bessel0
169 :
170 : REAL(KIND=dp), PARAMETER :: tol = 1.e-12_dp
171 :
172 : INTEGER :: i, ii, il, isvar, k
173 : REAL(KIND=dp) :: arg, fact, xsq
174 : REAL(KIND=dp), DIMENSION(4) :: trig
175 :
176 801600 : IF ((x*x - 0.5_dp) > l) THEN
177 390913 : arg = x - pio2*REAL(l, KIND=dp)
178 390913 : trig(1) = SIN(arg)/x
179 390913 : trig(2) = COS(arg)/x
180 390913 : trig(3) = -trig(1)
181 390913 : trig(4) = -trig(2)
182 390913 : bessel0 = trig(1)
183 390913 : IF (l /= 0) THEN
184 283552 : xsq = 0.5_dp/x
185 283552 : fact = 1._dp
186 843252 : DO k = 1, l
187 559700 : ii = MOD(k, 4) + 1
188 559700 : fact = fac(k + l)/fac(k)/fac(l - k)*xsq**k
189 843252 : bessel0 = bessel0 + fact*trig(ii)
190 : END DO
191 : END IF
192 : ELSE
193 : ! Taylor expansion for small arguments
194 : isvar = 1
195 1053387 : DO il = 1, l
196 1053387 : isvar = isvar*(2*il + 1)
197 : END DO
198 410687 : IF (l /= 0._dp) THEN
199 317648 : fact = x**l/REAL(isvar, KIND=dp)
200 : ELSE
201 93039 : fact = 1._dp/REAL(isvar, KIND=dp)
202 : END IF
203 410687 : bessel0 = fact
204 410687 : xsq = -0.5_dp*x*x
205 410687 : isvar = 2*l + 1
206 1193216 : DO i = 1, 1000
207 1193216 : isvar = isvar + 2
208 1193216 : fact = fact*xsq/REAL(i*isvar, KIND=dp)
209 1193216 : bessel0 = bessel0 + fact
210 1193216 : IF (ABS(fact) < tol) EXIT
211 : END DO
212 410687 : IF (ABS(fact) > tol) CPABORT("BESSEL0 NOT CONVERGED")
213 : END IF
214 :
215 801600 : END FUNCTION bessel0
216 :
217 : END MODULE bessel_lib
218 :
|
