From e1aefab23388470bc445506e2c9ccf682803f108 Mon Sep 17 00:00:00 2001 From: daniau Date: Mon, 18 Feb 2008 09:30:30 +0000 Subject: [PATCH] git-svn-id: https://lxsd.femto-st.fr/svn/fvn@38 b657c933-2333-4658-acf2-d3c7c2708721 --- fvn_fnlib/Makefile | 4 +- fvn_fnlib/besin.f90 | 79 +++++++++++++++++++++++++++++++++++++++ fvn_fnlib/besjn.f90 | 98 +++++++++++++++++++++++++++++++++++++++++++++++++ fvn_fnlib/beskn.f90 | 41 +++++++++++++++++++++ fvn_fnlib/besyn.f90 | 41 +++++++++++++++++++++ fvn_fnlib/dbesin.f90 | 80 ++++++++++++++++++++++++++++++++++++++++ fvn_fnlib/dbesjn.f90 | 98 +++++++++++++++++++++++++++++++++++++++++++++++++ fvn_fnlib/dbeskn.f90 | 41 +++++++++++++++++++++ fvn_fnlib/dbesyn.f90 | 40 ++++++++++++++++++++ fvn_fnlib/fvn_fnlib.f90 | 56 ++++++++++++++++++++++++++++ 10 files changed, 577 insertions(+), 1 deletion(-) create mode 100644 fvn_fnlib/besin.f90 create mode 100644 fvn_fnlib/besjn.f90 create mode 100644 fvn_fnlib/beskn.f90 create mode 100644 fvn_fnlib/besyn.f90 create mode 100644 fvn_fnlib/dbesin.f90 create mode 100644 fvn_fnlib/dbesjn.f90 create mode 100644 fvn_fnlib/dbeskn.f90 create mode 100644 fvn_fnlib/dbesyn.f90 diff --git a/fvn_fnlib/Makefile b/fvn_fnlib/Makefile index c9def1c..988177a 100644 --- a/fvn_fnlib/Makefile +++ b/fvn_fnlib/Makefile @@ -61,7 +61,9 @@ zasin.o zatan2.o zatanh.o zatan.o \ zbeta.o zcbrt.o zcosh.o zcot.o \ zexprl.o zgamma.o zgamr.o zlbeta.o \ zlngam.o zlnrel.o zlog10.o zpsi.o \ -zsinh.o ztanh.o ztan.o +zsinh.o ztanh.o ztan.o besyn.o \ +besjn.o dbesyn.o dbesjn.o beskn.o \ +besin.o dbeskn.o dbesin.o lib:$(library) diff --git a/fvn_fnlib/besin.f90 b/fvn_fnlib/besin.f90 new file mode 100644 index 0000000..5b6db4f --- /dev/null +++ b/fvn_fnlib/besin.f90 @@ -0,0 +1,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 diff --git a/fvn_fnlib/besjn.f90 b/fvn_fnlib/besjn.f90 new file mode 100644 index 0000000..b4ebb60 --- /dev/null +++ b/fvn_fnlib/besjn.f90 @@ -0,0 +1,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 diff --git a/fvn_fnlib/beskn.f90 b/fvn_fnlib/beskn.f90 new file mode 100644 index 0000000..cf5f906 --- /dev/null +++ b/fvn_fnlib/beskn.f90 @@ -0,0 +1,41 @@ +real(4) function beskn(n,x) + implicit none + ! This function compute the rank n Bessel Y function + ! using recurrence relation : + ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x) + ! + integer :: n + real(4) :: x + + real(4),external :: besk0,besk1 + real(4) :: two_on_x,bknm1,bkn,bktmp + integer :: i + + if (n==0) then + beskn=besk0(x) + return + end if + if (n==1) then + beskn=besk1(x) + return + end if + + if (n < 0) then + write(*,*) "Error in beskn, n must be >= 0" + stop + end if + if (x <= 0.) then + write(*,*) "Error in beskn, x must be strictly positive" + end if + + two_on_x=2./x + bknm1=besk0(x) + bkn=besk1(x) + + do i=1,n-1 + bktmp=two_on_x*bkn*i+bknm1 + bknm1=bkn + bkn=bktmp + end do + beskn=bktmp +end function diff --git a/fvn_fnlib/besyn.f90 b/fvn_fnlib/besyn.f90 new file mode 100644 index 0000000..9188969 --- /dev/null +++ b/fvn_fnlib/besyn.f90 @@ -0,0 +1,41 @@ +real(4) function besyn(n,x) + implicit none + ! This function compute the rank n Bessel Y function + ! using recurrence relation : + ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x) + ! + integer :: n + real(4) :: x + + real(4),external :: besy0,besy1 + real(4) :: two_on_x,bynm1,byn,bytmp + integer :: i + + if (n==0) then + besyn=besy0(x) + return + end if + if (n==1) then + besyn=besy1(x) + return + end if + + if (n < 0) then + write(*,*) "Error in besyn, n must be >= 0" + stop + end if + if (x <= 0.) then + write(*,*) "Error in besyn, x must be strictly positive" + end if + + two_on_x=2./x + bynm1=besy0(x) + byn=besy1(x) + + do i=1,n-1 + bytmp=two_on_x*byn*i-bynm1 + bynm1=byn + byn=bytmp + end do + besyn=bytmp +end function diff --git a/fvn_fnlib/dbesin.f90 b/fvn_fnlib/dbesin.f90 new file mode 100644 index 0000000..a71321f --- /dev/null +++ b/fvn_fnlib/dbesin.f90 @@ -0,0 +1,80 @@ +real(8) function dbesin(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(8) :: x + integer, optional :: factor + real(8), optional :: big + + integer :: tfactor + real(8) :: tbig,tsmall + real(8) :: two_on_x,binm1,bin,binp1,absx + integer :: i,start + real(8), external :: dbesi0,dbesi1 + + ! 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 + dbesin=dbesi0(x) + return + end if + if (n==1) then + dbesin=dbesi1(x) + return + end if + + if (n < 0) then + write(*,*) "Error in dbesin, n must be >= 0" + stop + end if + + absx=abs(x) + if (absx == 0.) then + dbesin=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 + dbesin=dbesin*tsmall + end if + if (i==n) dbesin=binp1 + end do + dbesin=dbesin*dbesi0(x)/bin + end if + ! if n is odd and x <0 + if ((x<0.) .and. (mod(n,2)==1)) dbesin=-dbesin + +end function diff --git a/fvn_fnlib/dbesjn.f90 b/fvn_fnlib/dbesjn.f90 new file mode 100644 index 0000000..da26cd1 --- /dev/null +++ b/fvn_fnlib/dbesjn.f90 @@ -0,0 +1,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 diff --git a/fvn_fnlib/dbeskn.f90 b/fvn_fnlib/dbeskn.f90 new file mode 100644 index 0000000..083c428 --- /dev/null +++ b/fvn_fnlib/dbeskn.f90 @@ -0,0 +1,41 @@ +real(8) function dbeskn(n,x) + implicit none + ! This function compute the rank n Bessel Y function + ! using recurrence relation : + ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x) + ! + integer :: n + real(8) :: x + + real(8),external :: dbesk0,dbesk1 + real(8) :: two_on_x,bknm1,bkn,bktmp + integer :: i + + if (n==0) then + dbeskn=dbesk0(x) + return + end if + if (n==1) then + dbeskn=dbesk1(x) + return + end if + + if (n < 0) then + write(*,*) "Error in dbeskn, n must be >= 0" + stop + end if + if (x <= 0.) then + write(*,*) "Error in dbeskn, x must be strictly positive" + end if + + two_on_x=2./x + bknm1=dbesk0(x) + bkn=dbesk1(x) + + do i=1,n-1 + bktmp=two_on_x*bkn*i+bknm1 + bknm1=bkn + bkn=bktmp + end do + dbeskn=bktmp +end function diff --git a/fvn_fnlib/dbesyn.f90 b/fvn_fnlib/dbesyn.f90 new file mode 100644 index 0000000..41b133f --- /dev/null +++ b/fvn_fnlib/dbesyn.f90 @@ -0,0 +1,40 @@ +real(8) function dbesyn(n,x) + implicit none + ! This function compute the rank n Bessel Y function + ! using recurrence relation : + ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x) + ! + integer :: n + real(8) :: x + + real(8),external :: dbesy0,dbesy1 + real(8) :: two_on_x,bynm1,byn,bytmp + integer :: i + + if (n==0) then + dbesyn=dbesy0(x) + return + end if + if (n==1) then + dbesyn=dbesy1(x) + return + end if + if (n < 0) then + write(*,*) "Error in dbesyn, n must be >= 0" + stop + end if + if (x <= 0.) then + write(*,*) "Error in dbesyn, x must be strictly positive" + end if + + two_on_x=2./x + bynm1=dbesy0(x) + byn=dbesy1(x) + + do i=1,n-1 + bytmp=two_on_x*byn*i-bynm1 + bynm1=byn + byn=bytmp + end do + dbesyn=bytmp +end function diff --git a/fvn_fnlib/fvn_fnlib.f90 b/fvn_fnlib/fvn_fnlib.f90 index bb6dee3..20079a1 100644 --- a/fvn_fnlib/fvn_fnlib.f90 +++ b/fvn_fnlib/fvn_fnlib.f90 @@ -733,6 +733,62 @@ interface bsk1e end function dbsk1e end interface bsk1e +! nth order J +interface bsjn + real(4) function besjn(n,x,factor,big) + integer(4) :: n + real(4) :: x + integer(4), optional :: factor + real(4), optional :: big + end function besjn + real(8) function dbesjn(n,x,factor,big) + integer(4) :: n + real(8) :: x + integer(4), optional :: factor + real(8), optional :: big + end function dbesjn +end interface bsjn + +! nth order Y +interface bsyn + real(4) function besyn(n,x) + integer(4) :: n + real(4) :: x + end function besyn + real(8) function dbesyn(n,x) + integer(4) :: n + real(8) :: x + end function dbesyn +end interface bsyn + +! nth order I +interface bsin + real(4) function besin(n,x,factor,big) + integer(4) :: n + real(4) :: x + integer(4), optional :: factor + real(4), optional :: big + end function besin + real(8) function dbesin(n,x,factor,big) + integer(4) :: n + real(8) :: x + integer(4), optional :: factor + real(8), optional :: big + end function dbesin +end interface bsin + +! nth order K +interface bskn + real(4) function beskn(n,x) + integer(4) :: n + real(4) :: x + end function beskn + real(8) function dbeskn(n,x) + integer(4) :: n + real(8) :: x + end function dbeskn +end interface bskn + !!!!!!!!!!!!!!!!!!!!! ! MISSING BSJNS ! MISSING BSINS -- 2.16.4