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