21 #include "../base/base_uses.f90"
37 REAL(kind=
dp),
INTENT(IN) :: x
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
49 bessk0 = (-log(x/2.0_dp)*bessi0(x)) + (p1 + y* &
50 (p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))))
53 bessk0 = (exp(-x)/sqrt(x))*(q1 + y*(q2 + y* &
54 (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
66 REAL(kind=
dp),
INTENT(IN) :: x
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
78 bessk1 = (log(x/2.0_dp)*bessi1(x)) + (1.0_dp/x)* &
79 (p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y* &
83 bessk1 = (exp(-x)/sqrt(x))*(q1 + y*(q2 + y* &
84 (q3 + y*(q4 + y*(q5 + y*(q6 + y*q7))))))
94 ELEMENTAL FUNCTION bessi0(x)
96 REAL(kind=
dp),
INTENT(IN) :: x
97 REAL(kind=
dp) :: bessi0
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, &
105 REAL(kind=
dp) :: ax, y
107 IF (abs(x) < 3.75_dp)
THEN
109 bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
110 (p5 + y*(p6 + y*p7)))))
114 bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y* &
115 (q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y* &
126 ELEMENTAL FUNCTION bessi1(x)
128 REAL(kind=
dp),
INTENT(IN) :: x
129 REAL(kind=
dp) :: bessi1
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, &
137 REAL(kind=
dp) :: ax, y
139 IF (abs(x) < 3.75_dp)
THEN
141 bessi1 = p1 + y*(p2 + y*(p3 + y*(p4 + y* &
142 (p5 + y*(p6 + y*p7)))))
146 bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y* &
147 (q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y* &
149 IF (x < 0.0_dp) bessi1 = -bessi1
166 REAL(kind=
dp),
INTENT(IN) :: x
167 INTEGER,
INTENT(IN) :: l
170 REAL(kind=
dp),
PARAMETER :: tol = 1.e-12_dp
172 INTEGER :: i, ii, il, isvar, k
173 REAL(kind=
dp) :: arg, fact, xsq
174 REAL(kind=
dp),
DIMENSION(4) :: trig
176 IF (x > sqrt(real(l, kind=
dp) + 0.5_dp))
THEN
177 arg = x - 0.5_dp*real(l, kind=
dp)*
pi
188 fact =
fac(k + l)/
fac(k)/
fac(l - k)*xsq**k
196 isvar = isvar*(2*il + 1)
199 fact = x**l/real(isvar, kind=
dp)
201 fact = 1._dp/real(isvar, kind=
dp)
208 fact = fact*xsq/real(i*isvar, kind=
dp)
210 IF (abs(fact) < tol)
EXIT
212 IF (abs(fact) > tol) cpabort(
"BESSEL0 NOT CONVERGED")
Calculates Bessel functions.
elemental real(kind=dp) function, public bessk1(x)
...
elemental impure real(kind=dp) function, public bessel0(x, l)
...
elemental real(kind=dp) function, public bessk0(x)
...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), dimension(0:maxfac), parameter, public fac