besjn.f90
2.74 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
real(4) function besjn(n,x,factor,big)
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
!
integer :: n
real(4) :: x
integer, optional :: factor
real(4), optional :: big
integer :: tfactor
real(4) :: tbig,tsmall,som
real(4),external :: besj0,besj1
real(4) :: two_on_x,bjnm1,bjn,bjnp1,absx
integer :: i,start
logical :: iseven
! Initialization of optional parameters
tfactor=40
if(present(factor)) tfactor=factor
tbig=1e10
if(present(big)) tbig=big
tsmall=1./tbig
if (n==0) then
besjn=besj0(x)
return
end if
if (n==1) then
besjn=besj1(x)
return
end if
if (n < 0) then
write(*,*) "Error in besjn, n must be >= 0"
stop
end if
absx=abs(x)
if (absx == 0.) then
besjn=0.
else if (absx > float(n)) then
! For x > n upward reccurence is stable
two_on_x=2./absx
bjnm1=besj0(absx)
bjn=besj1(absx)
do i=1,n-1
bjnp1=two_on_x*bjn*i-bjnm1
bjnm1=bjn
bjn=bjnp1
end do
besjn=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(float(n*tfactor))))/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
besjn=besjn*tsmall
som=som*tsmall
end if
if (iseven) then
som=som+bjn
end if
iseven= .not. iseven
if (i==n) besjn=bjnp1
end do
som=2.*som-bjn
besjn=besjn/som
end if
! if n is odd and x <0
if ((x<0.) .and. (mod(n,2)==1)) besjn=-besjn
end function