SFEMaNS  version 5.3
Reference documentation for SFEMaNS
bessel.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond, Raphael Laguerre, Caroline Nore, Copyright 2005
3 !
4 MODULE bessel
5 
6  PUBLIC :: bessk, bessk0, bessk1
7  PUBLIC :: bessi, bessi0, bessi1
8  PUBLIC :: bessj0, bessj1
9  PRIVATE
10 
11 CONTAINS
12  ! ----------------------------------------------------------------------
13 
14  FUNCTION bessk(N,X) RESULT(vv)
15  REAL *8 X,vv,TOX,BK,BKM,BKP
16  INTEGER :: N, J
17  ! ------------------------------------------------------------------------
18  ! CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE
19  ! D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0. ON UTILISE ICI LA
20  ! FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1.
21  !
22  ! THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD
23  ! KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE
24  ! CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1.
25  ! ------------------------------------------------------------------------
26  ! REFERENCE:
27  ! C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
28  ! MATHEMATICAL TABLES, VOL.5, 1962.
29  ! ------------------------------------------------------------------------
30  IF (n.EQ.0) THEN
31  vv = bessk0(x)
32  RETURN
33  ENDIF
34  IF (n.EQ.1) THEN
35  vv = bessk1(x)
36  RETURN
37  ENDIF
38  IF (x.EQ.0.d0) THEN
39  vv = 1.d30
40  RETURN
41  ENDIF
42  tox = 2.d0/x
43  bk = bessk1(x)
44  bkm = bessk0(x)
45  DO 11 j=1,n-1
46  bkp = bkm+dfloat(j)*tox*bk
47  bkm = bk
48  bk = bkp
49 11 CONTINUE
50  vv = bk
51  RETURN
52  END FUNCTION bessk
53 
54  ! ----------------------------------------------------------------------
55  FUNCTION bessk0(X) RESULT(vv)
56  ! CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
57  ! POUR TOUT X REEL NON NUL.
58  !
59  ! CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF
60  ! ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
61  ! ----------------------------------------------------------------------
62  REAL*8 X,vv,Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7
63  DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0, &
64  0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
65  DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1, &
66  -0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
67  IF(x.EQ.0.d0) THEN
68  vv=1.d30
69  RETURN
70  ENDIF
71  IF(x.LE.2.d0) THEN
72  y=x*x/4.d0
73  ax=-log(x/2.d0)*bessi0(x)
74  vv=ax+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
75  ELSE
76  y=(2.d0/x)
77  ax=exp(-x)/sqrt(x)
78  vv=ax*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))))
79  ENDIF
80  RETURN
81  END FUNCTION bessk0
82  ! ----------------------------------------------------------------------
83  FUNCTION bessk1(X) RESULT(vv)
84  ! CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
85  ! POUR TOUT X REEL POSITF NON NUL.
86  !
87  ! CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF
88  ! ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
89  ! ----------------------------------------------------------------------
90  REAL*8 X,vv,Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7
91  DATA p1,p2,p3,p4,p5,p6,p7/1.d0,0.15443144d0,-0.67278579d0, &
92  -0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
93  DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1, &
94  0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
95  IF(x.EQ.0.d0) THEN
96  vv=1.d32
97  RETURN
98  ENDIF
99  IF(x.LE.2.d0) THEN
100  y=x*x/4.d0
101  ax=log(x/2.d0)*bessi1(x)
102  vv=ax+(1.d0/x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
103  ELSE
104  y=(2.d0/x)
105  ax=exp(-x)/sqrt(x)
106  vv=ax*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))))
107  ENDIF
108  RETURN
109  END FUNCTION bessk1
110 
111 
112  FUNCTION bessi(N,X) RESULT(vv)
113  !
114  ! This subroutine calculates the first kind modified Bessel function
115  ! of integer order N, for any REAL X. We use here the classical
116  ! recursion formula, when X > N. For X < N, the Miller's algorithm
117  ! is used to avoid overflows.
118  ! REFERENCE:
119  ! C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
120  ! MATHEMATICAL TABLES, VOL.5, 1962.
121  !
122  IMPLICIT NONE
123  INTEGER, PARAMETER :: IACC = 40
124  REAL(KIND=8), PARAMETER :: BIGNO = 1.d10, bigni = 1.d-10
125  REAL(KIND=8) :: X,vv,TOX,BIM,BI,BIP
126  INTEGER :: n, m, j
127  IF (n.EQ.0) THEN
128  vv = bessi0(x)
129  RETURN
130  ENDIF
131  IF (n.EQ.1) THEN
132  vv = bessi1(x)
133  RETURN
134  ENDIF
135  IF(x.EQ.0.d0) THEN
136  vv=0.d0
137  RETURN
138  ENDIF
139  tox = 2.d0/x
140  bip = 0.d0
141  bi = 1.d0
142  vv = 0.d0
143  m = 2*((n+int(sqrt(float(iacc*n)))))
144  DO 12 j = m,1,-1
145  bim = bip+dfloat(j)*tox*bi
146  bip = bi
147  bi = bim
148  IF (abs(bi).GT.bigno) THEN
149  bi = bi*bigni
150  bip = bip*bigni
151  vv = vv*bigni
152  ENDIF
153  IF (j.EQ.n) vv = bip
154 12 CONTINUE
155  vv = vv*bessi0(x)/bi
156  RETURN
157  END FUNCTION bessi
158  ! ----------------------------------------------------------------------
159  ! ----------------------------------------------------------------------
160  ! Auxiliary Bessel functions for N=0, N=1
161  FUNCTION bessi0(X) RESULT(vv)
162  IMPLICIT NONE
163  REAL(KIND=8) :: X,vv,Y,P1,P2,P3,P4,P5,P6,P7, &
164  Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
165  DATA p1,p2,p3,p4,p5,p6,p7/1.d0,3.5156229d0,3.0899424d0,1.2067429d0, &
166  0.2659732d0,0.360768d-1,0.45813d-2/
167  DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, &
168  0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1, &
169  0.2635537d-1,-0.1647633d-1,0.392377d-2/
170  IF(abs(x).LT.3.75d0) THEN
171  y=(x/3.75d0)**2
172  vv=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
173  ELSE
174  ax=abs(x)
175  y=3.75d0/ax
176  bx=exp(ax)/sqrt(ax)
177  ax=q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))
178  vv=ax*bx
179  ENDIF
180  RETURN
181  END FUNCTION bessi0
182  ! ----------------------------------------------------------------------
183  FUNCTION bessi1(X) RESULT(vv)
184  IMPLICIT NONE
185  REAL(KIND=8) :: X,vv,Y,P1,P2,P3,P4,P5,P6,P7, &
186  Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
187  DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0, &
188  0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
189  DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, &
190  -0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1, &
191  -0.2895312d-1,0.1787654d-1,-0.420059d-2/
192  IF(abs(x).LT.3.75d0) THEN
193  y=(x/3.75d0)**2
194  vv=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
195  ELSE
196  ax=abs(x)
197  y=3.75d0/ax
198  bx=exp(ax)/sqrt(ax)
199  ax=q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))
200  vv=ax*bx
201  ENDIF
202  RETURN
203  END FUNCTION bessi1
204  !-----------------------------------------------------------------------
205  FUNCTION bessj0(x) RESULT(vv)
206  IMPLICIT NONE
207  REAL(KIND=8) :: x, y, ax, xx,vv , z
208  REAL(KIND=8) :: p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6
209  DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4, &
210  -.2073370639d-5,.2093887211d-6/
211  DATA q1,q2,q3,q4,q5/-.1562499995d-1, &
212  .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
213 
214  DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0, &
215  651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/
216  DATA s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0, &
217  59272.64853d0,267.8532712d0,1.d0/
218  !
219  IF(abs(x).LT.8.d0)THEN
220  y=x**2
221  vv=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y* &
222  (s4+y*(s5+y*s6)))))
223  ELSE
224  ax=abs(x)
225  z=8./ax
226  y=z**2
227  xx=ax-.785398164
228  vv=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* &
229  p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
230  ENDIF
231  RETURN
232  END FUNCTION bessj0
233  !-----------------------------------------------------------------------
234 
235  FUNCTION bessj1(x) RESULT(vv)
236  IMPLICIT NONE
237  REAL(KIND=8) :: x, y, ax, xx,vv , z
238  REAL(KIND=8) :: p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6
239  !April 18 2008, JLG
240  REAL(KIND=8) :: one
241  DATA one/1.d0/
242  !April 18 2008
243  DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0, &
244  242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/
245  DATA s1,s2,s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0, &
246  99447.43394d0,376.9991397d0,1.d0/
247  DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4, &
248  .2457520174d-5,-.240337019d-6/
249  DATA q1,q2,q3,q4,q5/.04687499995d0, &
250  -.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
251  !
252  IF(abs(x).LT.8.)THEN
253  y=x**2
254  vv=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+ &
255  y*(s4+y*(s5+y*s6)))))
256  ELSE
257  ax=abs(x)
258  z=8./ax
259  y=z**2
260  xx=ax-2.356194491
261  !April 18 2008: 1.d0 -> one
262  vv=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y* &
263  p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(one,x)
264  !April 18 2008
265  ENDIF
266  RETURN
267  END FUNCTION bessj1
268 
269 
270  END MODULE bessel
bessel::bessk0
real *8 function, public bessk0(X)
Definition: bessel.f90:56
bessel
Definition: bessel.f90:4
bessel::bessi0
real(kind=8) function, public bessi0(X)
Definition: bessel.f90:162
bessel::bessj0
real(kind=8) function, public bessj0(x)
Definition: bessel.f90:206
bessel::bessi1
real(kind=8) function, public bessi1(X)
Definition: bessel.f90:184
bessel::bessk
real *8 function, public bessk(N, X)
Definition: bessel.f90:15
bessel::bessk1
real *8 function, public bessk1(X)
Definition: bessel.f90:84
bessel::bessi
real(kind=8) function, public bessi(N, X)
Definition: bessel.f90:113
bessel::bessj1
real(kind=8) function, public bessj1(x)
Definition: bessel.f90:236