besin.f90
2.28 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
real(4) function besin(n,x,factor,big)
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
!
integer :: n
real(4) :: x
integer, optional :: factor
real(4), optional :: big
integer :: tfactor
real(4) :: tbig,tsmall
real(4) :: two_on_x,binm1,bin,binp1,absx
integer :: i,start
real(4), external :: besi0,besi1
! 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
besin=besi0(x)
return
end if
if (n==1) then
besin=besi1(x)
return
end if
if (n < 0) then
write(*,*) "Error in besin, n must be >= 0"
stop
end if
absx=abs(x)
if (absx == 0.) then
besin=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(float(n*tfactor))))/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
besin=besin*tsmall
end if
if (i==n) besin=binp1
end do
besin=besin*besi0(x)/bin
end if
! if n is odd and x <0
if ((x<0.) .and. (mod(n,2)==1)) besin=-besin
end function