-
to increase accuracy of computed values compared to dbesri and dbesrj
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@60 b657c933-2333-4658-acf2-d3c7c2708721
dbesin.f90
2.43 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
function dbesin(n,x,factor,big)
use fvn_common
implicit none
! This function compute the rank n Bessel J function
! using recurrence relation :
! In+1(x)=-2n/x * In(x) + In-1(x)
!
! Two optional parameters :
! factor : an integer that is used in Miller's algorithm to determine the
! starting point of iteration. Default value is 40, an increase of this value
! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
! big : a real that determine the threshold for taking anti overflow counter measure
! default value is 1e10
!
! 2009 : increasing factor defult value to 150
!
real(dp_kind) :: dbesin
integer :: n
real(dp_kind) :: x
integer, optional :: factor
real(dp_kind), optional :: big
integer :: tfactor
real(dp_kind) :: tbig,tsmall
real(dp_kind) :: two_on_x,binm1,bin,binp1,absx
integer :: i,start
real(dp_kind), external :: dbesi0,dbesi1
! Initialization of optional parameters
tfactor=150
if(present(factor)) tfactor=factor
tbig=1e10
if(present(big)) tbig=big
tsmall=1./tbig
if (n==0) then
dbesin=dbesi0(x)
return
end if
if (n==1) then
dbesin=dbesi1(x)
return
end if
if (n < 0) then
write(*,*) "Error in dbesin, n must be >= 0"
stop
end if
absx=abs(x)
if (absx == 0.) then
dbesin=0.
else
! We use Miller's Algorithm
! as upward reccurence is unstable.
! This is adapted from Numerical Recipes
! Principle : use of downward recurrence from an arbitrary
! higher than n value with an arbitrary seed,
! and then use the normalization formula :
! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
! call to besi0
two_on_x=2./absx
start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start
binp1=0.
bin=1.
do i=start,1,-1
! begin downward rec
binm1=two_on_x*bin*i+binp1
binp1=bin
bin=binm1
! Action to prevent overflow
if (abs(bin) > tbig) then
bin=bin*tsmall
binp1=binp1*tsmall
dbesin=dbesin*tsmall
end if
if (i==n) dbesin=binp1
end do
dbesin=dbesin*dbesi0(x)/bin
end if
! if n is odd and x <0
if ((x<0.) .and. (mod(n,2)==1)) dbesin=-dbesin
end function