git-svn-id: https://lxsd.femto-st.fr/svn/fvn@38 b657c933-2333-4658-acf2-d3c7c2708721

daniau
1 parent f26a262db0
Showing 10 changed files with 577 additions and 1 deletions Side-by-side Diff
fvn_fnlib/Makefile
fvn_fnlib/besin.f90
fvn_fnlib/besjn.f90
fvn_fnlib/beskn.f90
fvn_fnlib/besyn.f90
fvn_fnlib/dbesin.f90
fvn_fnlib/dbesjn.f90
fvn_fnlib/dbeskn.f90
fvn_fnlib/dbesyn.f90
fvn_fnlib/fvn_fnlib.f90
@@ -61,7 +61,9 @@
 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)
  
+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
+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
+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
+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
+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
+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
+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
+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
@@ -733,6 +733,62 @@
       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
...	...	@@ -61,7 +61,9 @@
61	61	zbeta.o zcbrt.o zcosh.o zcot.o \
62	62	zexprl.o zgamma.o zgamr.o zlbeta.o \
63	63	zlngam.o zlnrel.o zlog10.o zpsi.o \
64		-zsinh.o ztanh.o ztan.o
	64	+zsinh.o ztanh.o ztan.o besyn.o \
	65	+besjn.o dbesyn.o dbesjn.o beskn.o \
	66	+besin.o dbeskn.o dbesin.o
65	67
66	68	lib:$(library)
67	69
	1	+real(4) function besin(n,x,factor,big)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel J function
	4	+ ! using recurrence relation :
	5	+ ! In+1(x)=-2n/x * In(x) + In-1(x)
	6	+ !
	7	+ ! Two optional parameters :
	8	+ ! factor : an integer that is used in Miller's algorithm to determine the
	9	+ ! starting point of iteration. Default value is 40, an increase of this value
	10	+ ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
	11	+ ! big : a real that determine the threshold for taking anti overflow counter measure
	12	+ ! default value is 1e10
	13	+ !
	14	+ integer :: n
	15	+ real(4) :: x
	16	+ integer, optional :: factor
	17	+ real(4), optional :: big
	18	+
	19	+ integer :: tfactor
	20	+ real(4) :: tbig,tsmall
	21	+ real(4) :: two_on_x,binm1,bin,binp1,absx
	22	+ integer :: i,start
	23	+ real(4), external :: besi0,besi1
	24	+
	25	+ ! Initialization of optional parameters
	26	+ tfactor=40
	27	+ if(present(factor)) tfactor=factor
	28	+ tbig=1e10
	29	+ if(present(big)) tbig=big
	30	+ tsmall=1./tbig
	31	+
	32	+ if (n==0) then
	33	+ besin=besi0(x)
	34	+ return
	35	+ end if
	36	+ if (n==1) then
	37	+ besin=besi1(x)
	38	+ return
	39	+ end if
	40	+ if (n < 0) then
	41	+ write(,) "Error in besin, n must be >= 0"
	42	+ stop
	43	+ end if
	44	+
	45	+ absx=abs(x)
	46	+ if (absx == 0.) then
	47	+ besin=0.
	48	+ else
	49	+ ! We use Miller's Algorithm
	50	+ ! as upward reccurence is unstable.
	51	+ ! This is adapted from Numerical Recipes
	52	+ ! Principle : use of downward recurrence from an arbitrary
	53	+ ! higher than n value with an arbitrary seed,
	54	+ ! and then use the normalization formula :
	55	+ ! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
	56	+ ! call to besi0
	57	+ two_on_x=2./absx
	58	+ start=2((n+int(sqrt(float(ntfactor))))/2) ! even start
	59	+ binp1=0.
	60	+ bin=1.
	61	+ do i=start,1,-1
	62	+ ! begin downward rec
	63	+ binm1=two_on_xbini+binp1
	64	+ binp1=bin
	65	+ bin=binm1
	66	+ ! Action to prevent overflow
	67	+ if (abs(bin) > tbig) then
	68	+ bin=bin*tsmall
	69	+ binp1=binp1*tsmall
	70	+ besin=besin*tsmall
	71	+ end if
	72	+ if (i==n) besin=binp1
	73	+ end do
	74	+ besin=besin*besi0(x)/bin
	75	+ end if
	76	+ ! if n is odd and x <0
	77	+ if ((x<0.) .and. (mod(n,2)==1)) besin=-besin
	78	+
	79	+end function
	1	+real(4) function besjn(n,x,factor,big)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel J function
	4	+ ! using recurrence relation :
	5	+ ! Jn+1(x)=2n/x * Jn(x) - Jn-1(x)
	6	+ !
	7	+ ! Two optional parameters :
	8	+ ! factor : an integer that is used in Miller's algorithm to determine the
	9	+ ! starting point of iteration. Default value is 40, an increase of this value
	10	+ ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
	11	+ ! big : a real that determine the threshold for taking anti overflow counter measure
	12	+ ! default value is 1e10
	13	+ !
	14	+ integer :: n
	15	+ real(4) :: x
	16	+ integer, optional :: factor
	17	+ real(4), optional :: big
	18	+
	19	+ integer :: tfactor
	20	+ real(4) :: tbig,tsmall,som
	21	+ real(4),external :: besj0,besj1
	22	+ real(4) :: two_on_x,bjnm1,bjn,bjnp1,absx
	23	+ integer :: i,start
	24	+ logical :: iseven
	25	+
	26	+ ! Initialization of optional parameters
	27	+ tfactor=40
	28	+ if(present(factor)) tfactor=factor
	29	+ tbig=1e10
	30	+ if(present(big)) tbig=big
	31	+ tsmall=1./tbig
	32	+
	33	+ if (n==0) then
	34	+ besjn=besj0(x)
	35	+ return
	36	+ end if
	37	+ if (n==1) then
	38	+ besjn=besj1(x)
	39	+ return
	40	+ end if
	41	+ if (n < 0) then
	42	+ write(,) "Error in besjn, n must be >= 0"
	43	+ stop
	44	+ end if
	45	+
	46	+ absx=abs(x)
	47	+ if (absx == 0.) then
	48	+ besjn=0.
	49	+ else if (absx > float(n)) then
	50	+ ! For x > n upward reccurence is stable
	51	+ two_on_x=2./absx
	52	+ bjnm1=besj0(absx)
	53	+ bjn=besj1(absx)
	54	+ do i=1,n-1
	55	+ bjnp1=two_on_xbjni-bjnm1
	56	+ bjnm1=bjn
	57	+ bjn=bjnp1
	58	+ end do
	59	+ besjn=bjnp1
	60	+ else
	61	+ ! For x <= n we use Miller's Algorithm
	62	+ ! as upward reccurence is unstable.
	63	+ ! This is adapted from Numerical Recipes
	64	+ ! Principle : use of downward recurrence from an arbitrary
	65	+ ! higher than n value with an arbitrary seed,
	66	+ ! and then use the normalization formula :
	67	+ ! 1=J0+2J2+2J4+2J6+....
	68	+ two_on_x=2./absx
	69	+ start=2((n+int(sqrt(float(ntfactor))))/2) ! even start
	70	+ som=0.
	71	+ iseven=.false.
	72	+ bjnp1=0.
	73	+ bjn=1.
	74	+ do i=start,1,-1
	75	+ ! begin downward rec
	76	+ bjnm1=two_on_xbjni-bjnp1
	77	+ bjnp1=bjn
	78	+ bjn=bjnm1
	79	+ ! Action to prevent overflow
	80	+ if (abs(bjn) > tbig) then
	81	+ bjn=bjn*tsmall
	82	+ bjnp1=bjnp1*tsmall
	83	+ besjn=besjn*tsmall
	84	+ som=som*tsmall
	85	+ end if
	86	+ if (iseven) then
	87	+ som=som+bjn
	88	+ end if
	89	+ iseven= .not. iseven
	90	+ if (i==n) besjn=bjnp1
	91	+ end do
	92	+ som=2.*som-bjn
	93	+ besjn=besjn/som
	94	+ end if
	95	+ ! if n is odd and x <0
	96	+ if ((x<0.) .and. (mod(n,2)==1)) besjn=-besjn
	97	+
	98	+end function
	1	+real(4) function beskn(n,x)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel Y function
	4	+ ! using recurrence relation :
	5	+ ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x)
	6	+ !
	7	+ integer :: n
	8	+ real(4) :: x
	9	+
	10	+ real(4),external :: besk0,besk1
	11	+ real(4) :: two_on_x,bknm1,bkn,bktmp
	12	+ integer :: i
	13	+
	14	+ if (n==0) then
	15	+ beskn=besk0(x)
	16	+ return
	17	+ end if
	18	+ if (n==1) then
	19	+ beskn=besk1(x)
	20	+ return
	21	+ end if
	22	+
	23	+ if (n < 0) then
	24	+ write(,) "Error in beskn, n must be >= 0"
	25	+ stop
	26	+ end if
	27	+ if (x <= 0.) then
	28	+ write(,) "Error in beskn, x must be strictly positive"
	29	+ end if
	30	+
	31	+ two_on_x=2./x
	32	+ bknm1=besk0(x)
	33	+ bkn=besk1(x)
	34	+
	35	+ do i=1,n-1
	36	+ bktmp=two_on_xbkni+bknm1
	37	+ bknm1=bkn
	38	+ bkn=bktmp
	39	+ end do
	40	+ beskn=bktmp
	41	+end function
	1	+real(4) function besyn(n,x)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel Y function
	4	+ ! using recurrence relation :
	5	+ ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x)
	6	+ !
	7	+ integer :: n
	8	+ real(4) :: x
	9	+
	10	+ real(4),external :: besy0,besy1
	11	+ real(4) :: two_on_x,bynm1,byn,bytmp
	12	+ integer :: i
	13	+
	14	+ if (n==0) then
	15	+ besyn=besy0(x)
	16	+ return
	17	+ end if
	18	+ if (n==1) then
	19	+ besyn=besy1(x)
	20	+ return
	21	+ end if
	22	+
	23	+ if (n < 0) then
	24	+ write(,) "Error in besyn, n must be >= 0"
	25	+ stop
	26	+ end if
	27	+ if (x <= 0.) then
	28	+ write(,) "Error in besyn, x must be strictly positive"
	29	+ end if
	30	+
	31	+ two_on_x=2./x
	32	+ bynm1=besy0(x)
	33	+ byn=besy1(x)
	34	+
	35	+ do i=1,n-1
	36	+ bytmp=two_on_xbyni-bynm1
	37	+ bynm1=byn
	38	+ byn=bytmp
	39	+ end do
	40	+ besyn=bytmp
	41	+end function
	1	+real(8) function dbesin(n,x,factor,big)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel J function
	4	+ ! using recurrence relation :
	5	+ ! In+1(x)=-2n/x * In(x) + In-1(x)
	6	+ !
	7	+ ! Two optional parameters :
	8	+ ! factor : an integer that is used in Miller's algorithm to determine the
	9	+ ! starting point of iteration. Default value is 40, an increase of this value
	10	+ ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
	11	+ ! big : a real that determine the threshold for taking anti overflow counter measure
	12	+ ! default value is 1e10
	13	+ !
	14	+ integer :: n
	15	+ real(8) :: x
	16	+ integer, optional :: factor
	17	+ real(8), optional :: big
	18	+
	19	+ integer :: tfactor
	20	+ real(8) :: tbig,tsmall
	21	+ real(8) :: two_on_x,binm1,bin,binp1,absx
	22	+ integer :: i,start
	23	+ real(8), external :: dbesi0,dbesi1
	24	+
	25	+ ! Initialization of optional parameters
	26	+ tfactor=40
	27	+ if(present(factor)) tfactor=factor
	28	+ tbig=1e10
	29	+ if(present(big)) tbig=big
	30	+ tsmall=1./tbig
	31	+
	32	+ if (n==0) then
	33	+ dbesin=dbesi0(x)
	34	+ return
	35	+ end if
	36	+ if (n==1) then
	37	+ dbesin=dbesi1(x)
	38	+ return
	39	+ end if
	40	+
	41	+ if (n < 0) then
	42	+ write(,) "Error in dbesin, n must be >= 0"
	43	+ stop
	44	+ end if
	45	+
	46	+ absx=abs(x)
	47	+ if (absx == 0.) then
	48	+ dbesin=0.
	49	+ else
	50	+ ! We use Miller's Algorithm
	51	+ ! as upward reccurence is unstable.
	52	+ ! This is adapted from Numerical Recipes
	53	+ ! Principle : use of downward recurrence from an arbitrary
	54	+ ! higher than n value with an arbitrary seed,
	55	+ ! and then use the normalization formula :
	56	+ ! 1=I0-2I2+2I4-2I6+.... however it is easier to use a
	57	+ ! call to besi0
	58	+ two_on_x=2./absx
	59	+ start=2((n+int(sqrt(float(ntfactor))))/2) ! even start
	60	+ binp1=0.
	61	+ bin=1.
	62	+ do i=start,1,-1
	63	+ ! begin downward rec
	64	+ binm1=two_on_xbini+binp1
	65	+ binp1=bin
	66	+ bin=binm1
	67	+ ! Action to prevent overflow
	68	+ if (abs(bin) > tbig) then
	69	+ bin=bin*tsmall
	70	+ binp1=binp1*tsmall
	71	+ dbesin=dbesin*tsmall
	72	+ end if
	73	+ if (i==n) dbesin=binp1
	74	+ end do
	75	+ dbesin=dbesin*dbesi0(x)/bin
	76	+ end if
	77	+ ! if n is odd and x <0
	78	+ if ((x<0.) .and. (mod(n,2)==1)) dbesin=-dbesin
	79	+
	80	+end function
	1	+real(8) function dbesjn(n,x,factor,big)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel J function
	4	+ ! using recurrence relation :
	5	+ ! Jn+1(x)=2n/x * Jn(x) - Jn-1(x)
	6	+ !
	7	+ ! Two optional parameters :
	8	+ ! factor : an integer that is used in Miller's algorithm to determine the
	9	+ ! starting point of iteration. Default value is 40, an increase of this value
	10	+ ! will increase accuracy. Starting point ~ nearest even integer of sqrt(factor*n)
	11	+ ! big : a real that determine the threshold for taking anti overflow counter measure
	12	+ ! default value is 1e10
	13	+ !
	14	+ integer :: n
	15	+ real(8) :: x
	16	+ integer, optional :: factor
	17	+ real(8), optional :: big
	18	+
	19	+ integer :: tfactor
	20	+ real(8) :: tbig,tsmall,som
	21	+ real(8),external :: dbesj0,dbesj1
	22	+ real(8) :: two_on_x,bjnm1,bjn,bjnp1,absx
	23	+ integer :: i,start
	24	+ logical :: iseven
	25	+
	26	+ ! Initialization of optional parameters
	27	+ tfactor=40
	28	+ if(present(factor)) tfactor=factor
	29	+ tbig=1d10
	30	+ if(present(big)) tbig=big
	31	+ tsmall=1./tbig
	32	+
	33	+ if (n==0) then
	34	+ dbesjn=dbesj0(x)
	35	+ return
	36	+ end if
	37	+ if (n==1) then
	38	+ dbesjn=dbesj1(x)
	39	+ return
	40	+ end if
	41	+ if (n < 0) then
	42	+ write(,) "Error in dbesjn, n must be >= 0"
	43	+ stop
	44	+ end if
	45	+
	46	+ absx=abs(x)
	47	+ if (absx == 0.) then
	48	+ dbesjn=0.
	49	+ else if (absx > float(n)) then
	50	+ ! For x > n upward reccurence is stable
	51	+ two_on_x=2./absx
	52	+ bjnm1=dbesj0(absx)
	53	+ bjn=dbesj1(absx)
	54	+ do i=1,n-1
	55	+ bjnp1=two_on_xbjni-bjnm1
	56	+ bjnm1=bjn
	57	+ bjn=bjnp1
	58	+ end do
	59	+ dbesjn=bjnp1
	60	+ else
	61	+ ! For x <= n we use Miller's Algorithm
	62	+ ! as upward reccurence is unstable.
	63	+ ! This is adapted from Numerical Recipes
	64	+ ! Principle : use of downward recurrence from an arbitrary
	65	+ ! higher than n value with an arbitrary seed,
	66	+ ! and then use the normalization formula :
	67	+ ! 1=J0+2J2+2J4+2J6+....
	68	+ two_on_x=2./absx
	69	+ start=2((n+int(sqrt(float(ntfactor))))/2) ! even start
	70	+ som=0.
	71	+ iseven=.false.
	72	+ bjnp1=0.
	73	+ bjn=1.
	74	+ do i=start,1,-1
	75	+ ! begin downward rec
	76	+ bjnm1=two_on_xbjni-bjnp1
	77	+ bjnp1=bjn
	78	+ bjn=bjnm1
	79	+ ! Action to prevent overflow
	80	+ if (abs(bjn) > tbig) then
	81	+ bjn=bjn*tsmall
	82	+ bjnp1=bjnp1*tsmall
	83	+ dbesjn=dbesjn*tsmall
	84	+ som=som*tsmall
	85	+ end if
	86	+ if (iseven) then
	87	+ som=som+bjn
	88	+ end if
	89	+ iseven= .not. iseven
	90	+ if (i==n) dbesjn=bjnp1
	91	+ end do
	92	+ som=2.*som-bjn
	93	+ dbesjn=dbesjn/som
	94	+ end if
	95	+ ! if n is odd and x <0
	96	+ if ((x<0.) .and. (mod(n,2)==1)) dbesjn=-dbesjn
	97	+
	98	+end function
	1	+real(8) function dbeskn(n,x)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel Y function
	4	+ ! using recurrence relation :
	5	+ ! Kn+1(x)=2n/x * Kn(x) + Kn-1(x)
	6	+ !
	7	+ integer :: n
	8	+ real(8) :: x
	9	+
	10	+ real(8),external :: dbesk0,dbesk1
	11	+ real(8) :: two_on_x,bknm1,bkn,bktmp
	12	+ integer :: i
	13	+
	14	+ if (n==0) then
	15	+ dbeskn=dbesk0(x)
	16	+ return
	17	+ end if
	18	+ if (n==1) then
	19	+ dbeskn=dbesk1(x)
	20	+ return
	21	+ end if
	22	+
	23	+ if (n < 0) then
	24	+ write(,) "Error in dbeskn, n must be >= 0"
	25	+ stop
	26	+ end if
	27	+ if (x <= 0.) then
	28	+ write(,) "Error in dbeskn, x must be strictly positive"
	29	+ end if
	30	+
	31	+ two_on_x=2./x
	32	+ bknm1=dbesk0(x)
	33	+ bkn=dbesk1(x)
	34	+
	35	+ do i=1,n-1
	36	+ bktmp=two_on_xbkni+bknm1
	37	+ bknm1=bkn
	38	+ bkn=bktmp
	39	+ end do
	40	+ dbeskn=bktmp
	41	+end function
	1	+real(8) function dbesyn(n,x)
	2	+ implicit none
	3	+ ! This function compute the rank n Bessel Y function
	4	+ ! using recurrence relation :
	5	+ ! Yn+1(x)=2n/x * Yn(x) - Yn-1(x)
	6	+ !
	7	+ integer :: n
	8	+ real(8) :: x
	9	+
	10	+ real(8),external :: dbesy0,dbesy1
	11	+ real(8) :: two_on_x,bynm1,byn,bytmp
	12	+ integer :: i
	13	+
	14	+ if (n==0) then
	15	+ dbesyn=dbesy0(x)
	16	+ return
	17	+ end if
	18	+ if (n==1) then
	19	+ dbesyn=dbesy1(x)
	20	+ return
	21	+ end if
	22	+ if (n < 0) then
	23	+ write(,) "Error in dbesyn, n must be >= 0"
	24	+ stop
	25	+ end if
	26	+ if (x <= 0.) then
	27	+ write(,) "Error in dbesyn, x must be strictly positive"
	28	+ end if
	29	+
	30	+ two_on_x=2./x
	31	+ bynm1=dbesy0(x)
	32	+ byn=dbesy1(x)
	33	+
	34	+ do i=1,n-1
	35	+ bytmp=two_on_xbyni-bynm1
	36	+ bynm1=byn
	37	+ byn=bytmp
	38	+ end do
	39	+ dbesyn=bytmp
	40	+end function
...	...	@@ -733,6 +733,62 @@
733	733	end function dbsk1e
734	734	end interface bsk1e
735	735
	736	+! nth order J
	737	+interface bsjn
	738	+ real(4) function besjn(n,x,factor,big)
	739	+ integer(4) :: n
	740	+ real(4) :: x
	741	+ integer(4), optional :: factor
	742	+ real(4), optional :: big
	743	+ end function besjn
	744	+ real(8) function dbesjn(n,x,factor,big)
	745	+ integer(4) :: n
	746	+ real(8) :: x
	747	+ integer(4), optional :: factor
	748	+ real(8), optional :: big
	749	+ end function dbesjn
	750	+end interface bsjn
	751	+
	752	+! nth order Y
	753	+interface bsyn
	754	+ real(4) function besyn(n,x)
	755	+ integer(4) :: n
	756	+ real(4) :: x
	757	+ end function besyn
	758	+ real(8) function dbesyn(n,x)
	759	+ integer(4) :: n
	760	+ real(8) :: x
	761	+ end function dbesyn
	762	+end interface bsyn
	763	+
	764	+! nth order I
	765	+interface bsin
	766	+ real(4) function besin(n,x,factor,big)
	767	+ integer(4) :: n
	768	+ real(4) :: x
	769	+ integer(4), optional :: factor
	770	+ real(4), optional :: big
	771	+ end function besin
	772	+ real(8) function dbesin(n,x,factor,big)
	773	+ integer(4) :: n
	774	+ real(8) :: x
	775	+ integer(4), optional :: factor
	776	+ real(8), optional :: big
	777	+ end function dbesin
	778	+end interface bsin
	779	+
	780	+! nth order K
	781	+interface bskn
	782	+ real(4) function beskn(n,x)
	783	+ integer(4) :: n
	784	+ real(4) :: x
	785	+ end function beskn
	786	+ real(8) function dbeskn(n,x)
	787	+ integer(4) :: n
	788	+ real(8) :: x
	789	+ end function dbeskn
	790	+end interface bskn
	791	+
736	792	!!!!!!!!!!!!!!!!!!!!!
737	793	! MISSING BSJNS
738	794	! MISSING BSINS