-
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
dbesjn.f90
2.89 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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function dbesjn(n,x,factor,big)
use fvn_common
implicit none
! This function compute the rank n Bessel J function
! using recurrence relation :
! Jn+1(x)=2n/x * Jn(x) - Jn-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) :: dbesjn
integer :: n
real(dp_kind) :: x
integer, optional :: factor
real(dp_kind), optional :: big
integer :: tfactor
real(dp_kind) :: tbig,tsmall,som
real(dp_kind),external :: dbesj0,dbesj1
real(dp_kind) :: two_on_x,bjnm1,bjn,bjnp1,absx
integer :: i,start
logical :: iseven
! Initialization of optional parameters
tfactor=150
if(present(factor)) tfactor=factor
tbig=1d10
if(present(big)) tbig=big
tsmall=1./tbig
if (n==0) then
dbesjn=dbesj0(x)
return
end if
if (n==1) then
dbesjn=dbesj1(x)
return
end if
if (n < 0) then
write(*,*) "Error in dbesjn, n must be >= 0"
stop
end if
absx=abs(x)
if (absx == 0.) then
dbesjn=0.
else if (absx > real(n,dp_kind)) then
! For x > n upward reccurence is stable
two_on_x=2./absx
bjnm1=dbesj0(absx)
bjn=dbesj1(absx)
do i=1,n-1
bjnp1=two_on_x*bjn*i-bjnm1
bjnm1=bjn
bjn=bjnp1
end do
dbesjn=bjnp1
else
! For x <= n 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=J0+2J2+2J4+2J6+....
two_on_x=2./absx
start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start
som=0.
iseven=.false.
bjnp1=0.
bjn=1.
do i=start,1,-1
! begin downward rec
bjnm1=two_on_x*bjn*i-bjnp1
bjnp1=bjn
bjn=bjnm1
! Action to prevent overflow
if (abs(bjn) > tbig) then
bjn=bjn*tsmall
bjnp1=bjnp1*tsmall
dbesjn=dbesjn*tsmall
som=som*tsmall
end if
if (iseven) then
som=som+bjn
end if
iseven= .not. iseven
if (i==n) dbesjn=bjnp1
end do
som=2.*som-bjn
dbesjn=dbesjn/som
end if
! if n is odd and x <0
if ((x<0.) .and. (mod(n,2)==1)) dbesjn=-dbesjn
end function