!        Program: Example Inx.f90
!
!        Brian Geelen.   Created: november 2001.
!
!        Example use of Module Bessf_In_x
!
	PROGRAM example_Inx

	USE Bessf_In_x

    IMPLICIT NONE
    INTEGER, PARAMETER :: i15 = SELECTED_REAL_KIND(15,300)
	REAL (KIND = i15)  :: x, Bess_In

	INTEGER, PARAMETER  :: i9 = SELECTED_INT_KIND(9)
    INTEGER (KIND = i9) :: n, loop_1

    n = -1
    x =  1.0D0

        DO loop_1 = 1, 21

            n = n + 1

	        CALL Bessf_In(x, n, Bess_In)

		    WRITE(*,100)x, n, Bess_In
	    100 FORMAT(2x, F5.1, 4x, I3, 4x, ES16.9)

	    END DO

	STOP
	END PROGRAM example_Inx

	MODULE Bessf_In_x
!
!   Module Program: Bessf_In.f90
!   FN: Bessf_In    Version 1; Created: november 2001.
!
!   Purpose: evaluation of the modifed Bessel function of the first kind,
!   In(x), where n, an integer value, over the interval 0 to 100, and
!   x, is a real positive argument, over the interval 1 to 100. Calculation
!   procedure is based on convergent power series and asymptotic series.
!
!   Accuracy: returns precision of at least 9 significant decimal places.
!
!   Reference: Geelen, B., "Accurate solution for the modified Bessel
!   function of the first kind", Advances in Engineering Software, No 23,
!   1995, pp. 105-109.
!
!   Uses apart subroutines for I0(x), I1(x) and In(x) where n = or > 2.
!
!   Brian Geelen; PB10416, 6000GK Weert, The Netherlands. Email: bwave@iae.nl
!
!   Example output result listing:
!   x        n               In(x)
!   1.0      0     1.266065878E+00
!   1.0      1     5.651591040E-01
!   1.0      2     1.357476698E-01
!   1.0      3     2.216842492E-02
!   1.0     20     3.966835986E-25
!   1.0    100     8.473674008E-189
!
!***
!***
!***
!
    IMPLICIT NONE
    PRIVATE
    PUBLIC :: Bessf_In
    CONTAINS

    SUBROUTINE Bessf_In(x, n, Bess_In)

    INTEGER, PARAMETER  :: i15 = SELECTED_REAL_KIND(15,300)
    INTEGER, PARAMETER  :: i9  = SELECTED_INT_KIND(9)

    REAL (KIND = i15), INTENT(IN)   :: x
    REAL (KIND = i15), INTENT(OUT)  :: Bess_In
    INTEGER (KIND = i9), INTENT(IN) :: n

        IF (n < 1) THEN

            CALL Bessf_I0(x, Bess_In)

        ELSEIF( n < 2) THEN

            CALL Bessf_I1(x, Bess_In)

        ELSE

            CALL Bessf_In_gt2(x, n, Bess_In)

        ENDIF

    END SUBROUTINE Bessf_In
!***
!******************************************************************************
!***
	SUBROUTINE Bessf_I0(x, Bess_I0)
!
!   Subroutine Program: Bessf_I0.f90
!   FN: Bessf_I0_ver1                  Created: november 2001.
!
!   Version 1; Brian Geelen.
!
!   Purpose: evaluation of the modifed Bessel function, I0(x),
!   where x, is a positive real positive argument, over the interval 1 to 100.
!
!   Accuracy: returns precision of at least 9 significant decimal places.
!
!   Convergent power series solution for x < 14, asymptotic series for x > 14.
!
!
!	precision and variable declarations

    IMPLICIT NONE
    INTEGER, PARAMETER  :: i15 = SELECTED_REAL_KIND(15,300)
    REAL (KIND = i15)   :: ak, px, an, series_delta, series_rslt,    &
                           odd, seq_num, ser_num, fact, b1, b2, b3

	INTEGER, PARAMETER  :: i9 = SELECTED_INT_KIND(9)
	INTEGER (KIND = i9) :: loop_1, loop_2

    REAL (KIND = i15), PARAMETER ::	e = 2.718281828459045D0,         &
	                   pi = 3.141592653589793D0, error = 0.5D-12

    REAL (KIND = i15), INTENT(IN)  :: x
    REAL (KIND = i15), INTENT(OUT) :: Bess_I0

!   select power series or asympotic series solution

	   expression_select: IF (x < 14.0) THEN              ! power series for x < 14

          ak = 0.0D0; an = 1.0D0; series_rslt = 0.0D0     ! set variables

            series_solution: DO loop_1 = 1, 50

                ak             = ak + 2.0D0
                px             = x **(ak)
                an             = an * (ak * ak)
                series_delta   = (px / an)
                series_rslt    = series_rslt + (px / an)
                Bess_I0        = 1.0D0 + series_rslt

				IF (series_delta < error) EXIT        ! exit if error term satisfied

			END DO series_solution

	   ELSE		!  asymptotic series solution for x > 14

          odd   = 1.0D0; seq_num = 1.0D0; ser_num  = 1.0D0 ! set variables
          fact  = 2.0D0; an      = 2.0D0

          b1            = (e **x) / DSQRT(2.0D0 * pi * x)
          b2            = 8.0D0 * x
          series_rslt   = 1.0D0 + (1.0D0 / b2) + (9.0D0 / (2.0D0 * (b2**2.0D0) ) )

            asymptotic_series: DO loop_2 = 1,20

                odd             = odd + 2.0D0
                seq_num         = seq_num * (odd** 2.0D0)
                ser_num         = seq_num * ( (odd + 2.0D0) **2.0D0)

                an              = an + 1.0D0
                fact            = fact * an
                b3              = (8.0D0 * x) **an

                series_delta    = (ser_num / (fact * b3) )
                series_rslt     = series_rslt + (ser_num / (fact * b3) )
                Bess_I0         = b1 * series_rslt

            END DO asymptotic_series

	   END IF expression_select

	END SUBROUTINE Bessf_I0
!***
!******************************************************************************
!***
    SUBROUTINE Bessf_I1(x, Bess_I1)
!
!   Subroutine Program: Bessf_I1
!   FN: Bessf_I1_ver1.f90              Created: november 2001.
!
!   Version 1; Brian Geelen.
!
!   Purpose: evaluation of the modified Bessel function, I1(x),
!   where x, is a real positive argument, over the interval 1 to 100
!
!   Accuracy: returns precision of at least 9 significant decimal places.
!
!   Convergent power series solution for x < 14, asymptotic series for x > 14.
!
!
!   precision and variable declarations

    IMPLICIT NONE
    INTEGER, PARAMETER  :: i15 = SELECTED_REAL_KIND(15,300)
    REAL (KIND = i15)   :: ak, an, fact1, sx, a2, fact2, series_delta, &
                           odd, seq_num, ser_num, c1, c2, c3, series

    INTEGER, PARAMETER  :: i9 = SELECTED_INT_KIND(9)
    INTEGER (KIND = i9) :: loop_1, loop_2

   	REAL (KIND = i15), PARAMETER ::	e = 2.718281828459045D0,           &
	                   pi = 3.141592653589793D0, error = 0.5D-12

    REAL (KIND = i15), INTENT(IN)  :: x
    REAL (KIND = i15), INTENT(OUT) :: Bess_I1

!   select power series or asympotic series solution

	   expression_select: IF (x < 14.0) THEN      ! power series solution for x < 14

          ak = 1.0D0; an = 0.0D0; fact1 = 1.0D0   ! set variables

          Bess_I1 = (x / 2.0D0)

            series_solution: DO loop_1 = 1,31

                ak    = ak + 2.0D0
                sx    = x **(ak)

                a2    = 2.0D0**(ak)

                an    = an + 1.0D0
                fact1 = fact1 * an
                fact2 = fact1 * (an + 1.0D0)

                Bess_I1       = ( sx / (a2 * fact1 * fact2 ) ) +  Bess_I1
                series_delta  = ( sx / (a2 * fact1 * fact2 ) )

                IF (series_delta < error) EXIT    ! exit if error term satisfied

            END DO series_solution

       ELSE      !  asymptotic series solution for x > 14

          odd    = 1.0D0; seq_num = 1.0D0; ser_num = 1.0D0 ! set variables
          an     = 2.0D0; fact1 = 2.0D0

          c1     = (e **x) / DSQRT(2.0D0 * pi * x)
		  c2     = 8.0D0 * x
		  series = 1.0D0 - (3.0D0 / c2) - (15.0D0 / (2.0D0 * (c2**2.0D0) ))

            asymptotic_series: DO loop_2 = 1,29

                odd      = odd + 2.0D0
                seq_num  = (odd * odd) * seq_num
                ser_num  = seq_num * (odd + 2.0D0) * (odd + 4.0D0)

                an       = an + 1.0D0
                fact1    = fact1 * an

                c3       = (8.0D0 * x) **an
                series   = series - ( ser_num / ( fact1 * c3 ))

                Bess_I1  = series * c1

			END DO asymptotic_series

	   END IF expression_select

	END SUBROUTINE Bessf_I1
!***
!******************************************************************************
!***
	Subroutine Bessf_In_gt2(x, n, Bess_In)
!
!   Subroutine Program: Bessf_In.f90
!   FN: Bessf_In_GT2_ver1              Created: november 2001.
!
!   Version 1; Brian Geelen.
!
!   Purpose: evaluation of the modifed Bessel function, In(x),
!   where n, a integer value, is = or > 2, and x, a positive real
!   argument, over the interval 1 to 100.
!
!   Accuracy: returns precision of at least 9 significant decimal places.
!
!   Employs a downward recurrence approach with nornmalisation to I0(x).
!	External call made to subroutine Bessf_I0
!   to solve the modifed Bessel function I0(x).
!
!
!	precision and variable declarations

	IMPLICIT NONE
	INTEGER, PARAMETER  :: i15 = SELECTED_REAL_KIND(15,300)
	REAL (KIND = i15)   :: a1, a2, a3, a4, Bess_I0

	INTEGER, PARAMETER  :: i9 = SELECTED_INT_KIND(9)
	INTEGER (KIND = i9) :: j1, m, k, iexp

	REAL (KIND = i15), INTENT(IN)  :: x
	REAL (KIND = i15), INTENT(OUT) :: Bess_In

	INTEGER (KIND = i9), INTENT(IN) :: n

	iexp = maxexponent(x)/2           ! initialisation parameter

    a1 = 1.0D0; a2 = 0.0D0; a3 = 2.0D0 / x
    Bess_In = 0.0D0

    j1 = 200

    m = 2 * ( (n + INT (SQRT (real (j1 * n)) )))

		 recursion: DO k = m, 1, -1   ! perform downward recursion

			 a4 = a2 + k * a1 * a3
			 a2 = a1
			 a1 = a4

				 IF (exponent(a1) > iexp) THEN

						 Bess_In = scale(Bess_In, -iexp)
						 a1      = scale(a1, -iexp)
						 a2      = scale(a2, -iexp)

				 END IF
						 IF (k == n) Bess_In = a2

		 END DO recursion

!   solve Bessel function I0(x) with subroutine call
!   and perform normalisation to I0(x)

    CALL Bessf_I0(x, Bess_I0)

    Bess_In = Bess_In * Bess_I0 / a1

    END Subroutine Bessf_In_gt2
!***
!******************************************************************************
!***
END MODULE Bessf_In_x