function besin(n,x,factor,big)

      use fvn_common

      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

      !
      real(sp_kind) :: besin

      integer :: n

      real(sp_kind) :: x

      integer, optional :: factor

      real(sp_kind), optional :: big

  
      integer :: tfactor

      real(sp_kind) :: tbig,tsmall
      real(sp_kind) :: two_on_x,binm1,bin,binp1,absx

      integer :: i,start

      real(sp_kind), 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(real(n*tfactor,sp_kind))))/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