/*****************************************************************************
*                                                                            *
*   AIRY.SPL     Copyright (C) 2009 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Ai(z) and Bi(z) Airy functions                              *
*                                                                            *
*   Revisions:   25 Nov 2009  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_AIRY

    AIRY

    Purpose: Evaluates the Airy function the first or second kind.

    Syntax:  AIRY(k, z)

                 k - Optional, an integer. The kind of Airy function to
                     evaluate.

                       0: compute Ai(x) (default)
                       1: compute the derivative of Ai(x)
                       2: compute Bi(x)
                       3: compute the derivative of Bi(x)

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

    Returns:
             A scalar or series, Ai(x), Bi(x) or the derivatives.

    Example:
             airy(1.5)

             returns 0.0717490, the value of Ai(1.5).

    Example:
             airy(2, 1.5)

             returns 1.878942, the value of Bi(1.5).

    Example:
             W1: airy(0, -10..0.01..4)
             W2: airy(1, -10..0.01..4)
             W3: airy(2, -10..0.01..4)
             W4: airy(3, -10..0.01..4)

             For z between -10 and 4, W1 displays the values of Ai(z),
             W2 displays the derivative of Ai(z), W3 displays Bi(z) and
             W4 displays the derivative of Bi(z).

    Remarks:
             Airy functions are solutions to the differential equation:

                2
             d y  
             --- + z y = 0
               2
             dz 

             Ai(z) as a solution of the first kind. Bi(z) is a linearly 
             independent solution of the second kind.

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

   See Also:
             Besselh
             Besseli
             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


/* airy function of the first or second kind */
ITERATE airy(k, z, opt)
{
        local y, scalarval, type, id;

        (k, z, opt) = airy_parse_args(k, z, opt);

        if (not(isscalar(k)))
        {
                if (length(k) == 1)
                {
                        k = castreal(k);
                }
                else
                {
                        error("airy - Unknown Type, k must be an integer 0, 1, 2 or 3");
                }
        }

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

        switch (k)
        {
                case 0:
                        type = "a";
                        id   =  0;
                        break;

                case 1:
                        type = "a";
                        id   =  1;
                        break;

                case 2:
                        type = "b";
                        id   =  0;
                        break;

                case 3:
                        type = "b";
                        id   =  1;
                        break;

                default:
                        error(sprintf("airy - Unknown Type:%d, k must be 0, 1, 2 or 3", k));
        }

        /* use built-in */
        y = besselx(type, {z}, {id}, opt);

        if (scalarval) y = y[1];

        return(y);
}


airy_parse_args(k, z, opt)
{
        if (argc < 3)
        {
                opt = 0;
                
                if (argc < 2)
                {
                        z = k;
                        k = 0;
                        
                        if (argc < 1) error("airy - input required");
                }
        }
        
        return(k, z, opt);
}