(git:0de0cc2)
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 ! **************************************************************************************************
16 MODULE 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 
28 CONTAINS
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 
217 END MODULE bessel_lib
218 
Calculates Bessel functions.
Definition: bessel_lib.F:16
elemental real(kind=dp) function, public bessk1(x)
...
Definition: bessel_lib.F:65
elemental impure real(kind=dp) function, public bessel0(x, l)
...
Definition: bessel_lib.F:161
elemental real(kind=dp) function, public bessk0(x)
...
Definition: bessel_lib.F:36
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37