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

#include 

#if @HELP_BESSELY

    BESSELY

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

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

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

    Example:
             bessely(0, 3)

             returns 0.376850, the value of Y0(3).

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

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

             {-inf, -3.323825, -1.780872, -1.260391, -0.978144, -0.781213}


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

             Evaluates the 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:

             Bessel functions are solutions to the differential equation:

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

             Yv(z) is a solution of the second kind. For for non-integer v:

                         Jv(z) cos(v pi) - J-v(z)
                Yv(z) =  ------------------------
                                sin(v pi)

             See YN for a faster evaluation of the bessel function of the
             second kind for strictly integer orders.

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

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

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




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

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

        if (scalarval) y = y[1];

        return(y);
}