(git:374b731)
Loading...
Searching...
No Matches
bessel_lib.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 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! **************************************************************************************************
16MODULE bessel_lib
17
18 USE kinds, ONLY: dp
19 USE mathconstants, ONLY: fac,&
20 pi
21#include "../base/base_uses.f90"
22
23 IMPLICIT NONE
24
25 PRIVATE
26 PUBLIC :: bessk0, bessk1, bessel0
27
28CONTAINS
29
30! **************************************************************************************************
31!> \brief ...
32!> \param x must be positive
33!> \return ...
34! **************************************************************************************************
35 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 IF (x < 2.0_dp) THEN
48 y = x*x/4.0_dp
49 bessk0 = (-log(x/2.0_dp)*bessi0(x)) + (p1 + y* &
50 (p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))))
51 ELSE
52 y = (2.0_dp/x)
53 bessk0 = (exp(-x)/sqrt(x))*(q1 + y*(q2 + y* &
54 (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
55 END IF
56
57 END FUNCTION bessk0
58
59! **************************************************************************************************
60!> \brief ...
61!> \param x must be positive
62!> \return ...
63! **************************************************************************************************
64 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 IF (x < 2.0_dp) THEN
77 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 (p6 + y*p7))))))
81 ELSE
82 y = 2.0_dp/x
83 bessk1 = (exp(-x)/sqrt(x))*(q1 + y*(q2 + y* &
84 (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
85 END IF
86
87 END FUNCTION bessk1
88
89! **************************************************************************************************
90!> \brief ...
91!> \param x ...
92!> \return ...
93! **************************************************************************************************
94 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 IF (abs(x) < 3.75_dp) THEN
108 y = (x/3.75_dp)**2
109 bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
110 (p5 + y*(p6 + y*p7)))))
111 ELSE
112 ax = abs(x)
113 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 (q8 + y*q9))))))))
117 END IF
118
119 END FUNCTION bessi0
120
121! **************************************************************************************************
122!> \brief ...
123!> \param x ...
124!> \return ...
125! **************************************************************************************************
126 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 IF (abs(x) < 3.75_dp) THEN
140 y = (x/3.75_dp)**2
141 bessi1 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
142 (p5 + y*(p6 + y*p7)))))
143 ELSE
144 ax = abs(x)
145 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 (q8 + y*q9))))))))
149 IF (x < 0.0_dp) bessi1 = -bessi1
150 END IF
151
152 END FUNCTION bessi1
153
154! **************************************************************************************************
155!> \brief ...
156!> \param x ...
157!> \param l ...
158!> \return ...
159! **************************************************************************************************
160 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 IF (x > sqrt(real(l, kind=dp) + 0.5_dp)) THEN
177 arg = x - 0.5_dp*real(l, kind=dp)*pi
178 trig(1) = sin(arg)/x
179 trig(2) = cos(arg)/x
180 trig(3) = -trig(1)
181 trig(4) = -trig(2)
182 bessel0 = trig(1)
183 IF (l /= 0) THEN
184 xsq = 0.5_dp/x
185 fact = 1._dp
186 DO k = 1, l
187 ii = mod(k, 4) + 1
188 fact = fac(k + l)/fac(k)/fac(l - k)*xsq**k
189 bessel0 = bessel0 + fact*trig(ii)
190 END DO
191 END IF
192 ELSE
193 ! Taylor expansion for small arguments
194 isvar = 1
195 DO il = 1, l
196 isvar = isvar*(2*il + 1)
197 END DO
198 IF (l /= 0._dp) THEN
199 fact = x**l/real(isvar, kind=dp)
200 ELSE
201 fact = 1._dp/real(isvar, kind=dp)
202 END IF
203 bessel0 = fact
204 xsq = -0.5_dp*x*x
205 isvar = 2*l + 1
206 DO i = 1, 1000
207 isvar = isvar + 2
208 fact = fact*xsq/real(i*isvar, kind=dp)
209 bessel0 = bessel0 + fact
210 IF (abs(fact) < tol) EXIT
211 END DO
212 IF (abs(fact) > tol) cpabort("BESSEL0 NOT CONVERGED")
213 END IF
214
215 END FUNCTION bessel0
216
217END MODULE bessel_lib
218
bessel_lib
Calculates Bessel functions.
Definition bessel_lib.F:16
bessel_lib::bessk1
elemental real(kind=dp) function, public bessk1(x)
...
Definition bessel_lib.F:65
bessel_lib::bessel0
elemental impure real(kind=dp) function, public bessel0(x, l)
...
Definition bessel_lib.F:161
bessel_lib::bessk0
elemental real(kind=dp) function, public bessk0(x)
...
Definition bessel_lib.F:36
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37