lxsd / fvn

Blame view

fvn_fnlib/dbesjn.f90 2.76 KB
edit raw normal view history
e1aefab23   daniau   git-svn-id: https... 
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(8) function dbesjn(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(8) :: x
      integer, optional :: factor
      real(8), optional :: big
  
      integer :: tfactor
      real(8) :: tbig,tsmall,som
      real(8),external :: dbesj0,dbesj1
      real(8) :: two_on_x,bjnm1,bjn,bjnp1,absx
      integer :: i,start
      logical :: iseven
  
      ! Initialization of optional parameters
      tfactor=40
      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 > float(n)) 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(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
                  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