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