lxsd / fvn

Download as
Browse Code ยป

Commit e80b2ec787b1e4e1a6c894f2c274e6814d46fd81

Authored by wdaniau
1 parent 8d883e8a1a
Exists in master and in 3 other branches geevx, geevx_2020, multi

Increase factor default value from 40 to 150 in dbesin.f90 and dbesjn.f90 

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

Showing 2 changed files with 6 additions and 2 deletions Inline Diff

fvn_fnlib/dbesin.f90
Diff comments   View file @ e80b2ec
function dbesin(n,x,factor,big) 1 1 function dbesin(n,x,factor,big)
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
! This function compute the rank n Bessel J function 4 4 ! This function compute the rank n Bessel J function
! using recurrence relation : 5 5 ! using recurrence relation :
! In+1(x)=-2n/x * In(x) + In-1(x) 6 6 ! In+1(x)=-2n/x * In(x) + In-1(x)
! 7 7 !
! Two optional parameters : 8 8 ! Two optional parameters :
! factor : an integer that is used in Miller's algorithm to determine the 9 9 ! 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 10 10 ! 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) 11 11 ! 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 12 12 ! big : a real that determine the threshold for taking anti overflow counter measure
! default value is 1e10 13 13 ! default value is 1e10
! 14 14 !
15 ! 2009 : increasing factor defult value to 150
16 !
real(dp_kind) :: dbesin 15 17 real(dp_kind) :: dbesin
integer :: n 16 18 integer :: n
real(dp_kind) :: x 17 19 real(dp_kind) :: x
integer, optional :: factor 18 20 integer, optional :: factor
real(dp_kind), optional :: big 19 21 real(dp_kind), optional :: big
20 22
integer :: tfactor 21 23 integer :: tfactor
real(dp_kind) :: tbig,tsmall 22 24 real(dp_kind) :: tbig,tsmall
real(dp_kind) :: two_on_x,binm1,bin,binp1,absx 23 25 real(dp_kind) :: two_on_x,binm1,bin,binp1,absx
integer :: i,start 24 26 integer :: i,start
real(dp_kind), external :: dbesi0,dbesi1 25 27 real(dp_kind), external :: dbesi0,dbesi1
26 28
! Initialization of optional parameters 27 29 ! Initialization of optional parameters
tfactor=40 28 30 tfactor=150
if(present(factor)) tfactor=factor 29 31 if(present(factor)) tfactor=factor
tbig=1e10 30 32 tbig=1e10
if(present(big)) tbig=big 31 33 if(present(big)) tbig=big
tsmall=1./tbig 32 34 tsmall=1./tbig
33 35
if (n==0) then 34 36 if (n==0) then
dbesin=dbesi0(x) 35 37 dbesin=dbesi0(x)
return 36 38 return
end if 37 39 end if
if (n==1) then 38 40 if (n==1) then
dbesin=dbesi1(x) 39 41 dbesin=dbesi1(x)
return 40 42 return
end if 41 43 end if
42 44
if (n < 0) then 43 45 if (n < 0) then
write(*,*) "Error in dbesin, n must be >= 0" 44 46 write(*,*) "Error in dbesin, n must be >= 0"
stop 45 47 stop
end if 46 48 end if
47 49
absx=abs(x) 48 50 absx=abs(x)
if (absx == 0.) then 49 51 if (absx == 0.) then
dbesin=0. 50 52 dbesin=0.
else 51 53 else
! We use Miller's Algorithm 52 54 ! We use Miller's Algorithm
! as upward reccurence is unstable. 53 55 ! as upward reccurence is unstable.
! This is adapted from Numerical Recipes 54 56 ! This is adapted from Numerical Recipes
! Principle : use of downward recurrence from an arbitrary 55 57 ! Principle : use of downward recurrence from an arbitrary
! higher than n value with an arbitrary seed, 56 58 ! higher than n value with an arbitrary seed,
! and then use the normalization formula : 57 59 ! and then use the normalization formula :
! 1=I0-2I2+2I4-2I6+.... however it is easier to use a 58 60 ! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
! call to besi0 59 61 ! call to besi0
two_on_x=2./absx 60 62 two_on_x=2./absx
start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start 61 63 start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start
binp1=0. 62 64 binp1=0.
bin=1. 63 65 bin=1.
do i=start,1,-1 64 66 do i=start,1,-1
! begin downward rec 65 67 ! begin downward rec
binm1=two_on_x*bin*i+binp1 66 68 binm1=two_on_x*bin*i+binp1
binp1=bin 67 69 binp1=bin
bin=binm1 68 70 bin=binm1
! Action to prevent overflow 69 71 ! Action to prevent overflow
if (abs(bin) > tbig) then 70 72 if (abs(bin) > tbig) then
bin=bin*tsmall 71 73 bin=bin*tsmall
binp1=binp1*tsmall 72 74 binp1=binp1*tsmall
dbesin=dbesin*tsmall 73 75 dbesin=dbesin*tsmall
end if 74 76 end if
if (i==n) dbesin=binp1 75 77 if (i==n) dbesin=binp1
end do 76 78 end do
dbesin=dbesin*dbesi0(x)/bin 77 79 dbesin=dbesin*dbesi0(x)/bin
fvn_fnlib/dbesjn.f90
Diff comments   View file @ e80b2ec
function dbesjn(n,x,factor,big) 1 1 function dbesjn(n,x,factor,big)
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
! This function compute the rank n Bessel J function 4 4 ! This function compute the rank n Bessel J function
! using recurrence relation : 5 5 ! using recurrence relation :
! Jn+1(x)=2n/x * Jn(x) - Jn-1(x) 6 6 ! Jn+1(x)=2n/x * Jn(x) - Jn-1(x)
! 7 7 !
! Two optional parameters : 8 8 ! Two optional parameters :
! factor : an integer that is used in Miller's algorithm to determine the 9 9 ! 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 10 10 ! 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) 11 11 ! 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 12 12 ! big : a real that determine the threshold for taking anti overflow counter measure
! default value is 1e10 13 13 ! default value is 1e10
! 14 14 !
15 ! 2009 : increasing factor defult value to 150
16 !
real(dp_kind) :: dbesjn 15 17 real(dp_kind) :: dbesjn
integer :: n 16 18 integer :: n
real(dp_kind) :: x 17 19 real(dp_kind) :: x
integer, optional :: factor 18 20 integer, optional :: factor
real(dp_kind), optional :: big 19 21 real(dp_kind), optional :: big
20 22
integer :: tfactor 21 23 integer :: tfactor
real(dp_kind) :: tbig,tsmall,som 22 24 real(dp_kind) :: tbig,tsmall,som
real(dp_kind),external :: dbesj0,dbesj1 23 25 real(dp_kind),external :: dbesj0,dbesj1
real(dp_kind) :: two_on_x,bjnm1,bjn,bjnp1,absx 24 26 real(dp_kind) :: two_on_x,bjnm1,bjn,bjnp1,absx
integer :: i,start 25 27 integer :: i,start
logical :: iseven 26 28 logical :: iseven
27 29
! Initialization of optional parameters 28 30 ! Initialization of optional parameters
tfactor=40 29 31 tfactor=150
if(present(factor)) tfactor=factor 30 32 if(present(factor)) tfactor=factor
tbig=1d10 31 33 tbig=1d10
if(present(big)) tbig=big 32 34 if(present(big)) tbig=big
tsmall=1./tbig 33 35 tsmall=1./tbig
34 36
if (n==0) then 35 37 if (n==0) then
dbesjn=dbesj0(x) 36 38 dbesjn=dbesj0(x)
return 37 39 return
end if 38 40 end if
if (n==1) then 39 41 if (n==1) then
dbesjn=dbesj1(x) 40 42 dbesjn=dbesj1(x)
return 41 43 return
end if 42 44 end if
if (n < 0) then 43 45 if (n < 0) then
write(*,*) "Error in dbesjn, n must be >= 0" 44 46 write(*,*) "Error in dbesjn, n must be >= 0"
stop 45 47 stop
end if 46 48 end if
47 49
absx=abs(x) 48 50 absx=abs(x)
if (absx == 0.) then 49 51 if (absx == 0.) then
dbesjn=0. 50 52 dbesjn=0.
else if (absx > real(n,dp_kind)) then 51 53 else if (absx > real(n,dp_kind)) then
! For x > n upward reccurence is stable 52 54 ! For x > n upward reccurence is stable
two_on_x=2./absx 53 55 two_on_x=2./absx
bjnm1=dbesj0(absx) 54 56 bjnm1=dbesj0(absx)
bjn=dbesj1(absx) 55 57 bjn=dbesj1(absx)
do i=1,n-1 56 58 do i=1,n-1
bjnp1=two_on_x*bjn*i-bjnm1 57 59 bjnp1=two_on_x*bjn*i-bjnm1
bjnm1=bjn 58 60 bjnm1=bjn
bjn=bjnp1 59 61 bjn=bjnp1
end do 60 62 end do
dbesjn=bjnp1 61 63 dbesjn=bjnp1
else 62 64 else
! For x <= n we use Miller's Algorithm 63 65 ! For x <= n we use Miller's Algorithm
! as upward reccurence is unstable. 64 66 ! as upward reccurence is unstable.
! This is adapted from Numerical Recipes 65 67 ! This is adapted from Numerical Recipes
! Principle : use of downward recurrence from an arbitrary 66 68 ! Principle : use of downward recurrence from an arbitrary
! higher than n value with an arbitrary seed, 67 69 ! higher than n value with an arbitrary seed,
! and then use the normalization formula : 68 70 ! and then use the normalization formula :
! 1=J0+2J2+2J4+2J6+.... 69 71 ! 1=J0+2J2+2J4+2J6+....
two_on_x=2./absx 70 72 two_on_x=2./absx
start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start 71 73 start=2*((n+int(sqrt(real(n*tfactor,sp_kind))))/2) ! even start
som=0. 72 74 som=0.
iseven=.false. 73 75 iseven=.false.
bjnp1=0. 74 76 bjnp1=0.
bjn=1. 75 77 bjn=1.
do i=start,1,-1 76 78 do i=start,1,-1
! begin downward rec 77 79 ! begin downward rec
bjnm1=two_on_x*bjn*i-bjnp1 78 80 bjnm1=two_on_x*bjn*i-bjnp1
bjnp1=bjn 79 81 bjnp1=bjn
bjn=bjnm1 80 82 bjn=bjnm1
! Action to prevent overflow 81 83 ! Action to prevent overflow
if (abs(bjn) > tbig) then 82 84 if (abs(bjn) > tbig) then
bjn=bjn*tsmall 83 85 bjn=bjn*tsmall
bjnp1=bjnp1*tsmall 84 86 bjnp1=bjnp1*tsmall
dbesjn=dbesjn*tsmall 85 87 dbesjn=dbesjn*tsmall
som=som*tsmall 86 88 som=som*tsmall
end if 87 89 end if
if (iseven) then 88 90 if (iseven) then
som=som+bjn 89 91 som=som+bjn
end if 90 92 end if
iseven= .not. iseven 91 93 iseven= .not. iseven
if (i==n) dbesjn=bjnp1 92 94 if (i==n) dbesjn=bjnp1
end do 93 95 end do
som=2.*som-bjn 94 96 som=2.*som-bjn
dbesjn=dbesjn/som 95 97 dbesjn=dbesjn/som
end if 96 98 end if