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