/*****************************************************************************
*                                                                            *
*   BESSELI.SPL  Copyright (C) 2009 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Iv(z) bessel function of fractional order                   *
*                                                                            *
*   Revisions:    8 Nov 2009  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_BESSELI

    BESSELI

    Purpose: Evaluates the modified Bessel function of the first kind, Iv(z)

    Syntax:  BESSELI(v, z, opt)

                 v - A real or real series, the order. The order must
                     be real but need not be an integer.

                 z - Any scalar or series. The input value, can be complex.

               opt - Optional. An integer, the scaling method.

                       0: no scaling (default)
                       1: scale by exp(-mag(real(z)))

    Returns:
             A scalar or series, the value of Iv(z), where v is the order
             and z is the input.

    Example:
             besseli(0, 3)

             returns 4.880793, the value of I0(3).

    Example:
             W1: 0..0.2..1
             W2: besseli(1, w1)

             Evaluates I1(z) for z bewteen 0 and 1. W2 contains the series:

             {0.0, 0.100501, 0.204027, 0.313704, 0.432865, 0.565159}

    Example:
             besseli(3..9, 0..0.2..10)

             Evaluates the modified bessel function of the first kind
             for orders 3 through 9 with inputs from 0 to 10.  Each
             column of the result contains the output for the specified
             order.

    Remarks:
             Modified Bessel functions are solutions to the differential 
             equation:

                 2
              2 d y     dy     2   2
             z  --- + z -- - (z + v ) y = 0
                  2
                dz      dz

             Iv(z) as a solution of the first kind. Kv(z) is a linearly 
             independent solution of the second kind.

             If the order is not an integer and z is real where z < 0, 
             the result is generally complex.

             BESSELI is based on a FORTRAN library written by Donald E. Amos.

   See Also:
             Besselh
             Besselj
             Besselk
             Bessely
             Jn
             Yn

 References:
             [1] Abramowitz and Stegun
                 Handbook of Mathematical Functions (9th printing 1970)
                 US Gov. Printing Office
                 Section 9.1.1, 9.1.89, 9.12

             [2] Amos, D.E.
                 A Subroutine Package for Bessel Functions of a Complex 
                 Argument and Nonnegative Order
                 Sandia National Laboratory Report
                 SAND85-1018, May, 1985. 

             [3] Amos, D.E. 
                 A Portable Package for Bessel Functions of a Complex
                 Argument and Nonnegative Order
                 Trans. Math. Software, 1986

#endif


/* modified bessel function of the second kind */
ITERATE besseli(v, z, opt)
{
        local y, scalarval;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("besseli - input required");
                        
                        z = v;
                        v = 0;
                }
                
                opt = 0;
        }

        /* is input a single value */
        scalarval = not(isarray(v)) && not(isarray(z));

        /* use built-in */
        y = besselx("i", {z}, {v}, opt);

        if (scalarval) y = y[1];

        return(y);
}