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

#include 

#if @HELP_BESSELK

    BESSELK

    Purpose: Evaluates the modified Bessel function of the second kind, Kv(z)

    Syntax:  BESSELK(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(z)

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

    Example:
             besselk(0, 3)

             returns 0.034740, the value of K0(3).

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

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

             {inf, 4.775973, 2.184354, 1.302835, 0.861782, 0.601907}


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

             Evaluates the modified bessel function of the second 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

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

             If z is real and z < 0, the result is generally complex.

             BESSELK 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 besselk(v, z, opt)
{
        local y, scalarval;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("besselk - 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("k", {z}, {v}, opt);

        if (scalarval) y = y[1];

        return(y);
}