/*****************************************************************************
*                                                                            *
*   GAMMAINC.SPL Copyright (C) 2008 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    incomplete gammma function                                  *
*                                                                            *
*   Revisions:   31 Mar 2008  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_GAMMAINC

    GAMMAINC

    Purpose: Evaluates the incomplete gamma function

    Syntax:  GAMMAINC(x, a, tail)

                 x - A series or scalar, the integration limit

                 a - Optional, a scalar. the exponential factor.
                     Defaults to 1.0.

              tail - Optional, a string. Specifies whether the upper
                     or lower function is evaluated.

                       'lower': lower incomplete gamma (default)
                       'upper': upper incomplete gamma

    Returns: A series or scalar.

    Example:
             gammainc(10, 2)

             returns 0.999501.

    Example:
             gammainc(3+2i, 2)

             returns 0.992332 + 0.222522i

    Example:
             gammainc(0..0.01..10, 2)

             plots the incomplete gamma function with x ranging from
             0 to 10 and a == 2.

    Remarks:
             The incomplete gamma function is defined as gammainc(x, a) =

              (1/gamma(a)) * integral from 0 to x of t^2 * (a-1) exp(-t) dt

             For a real and a >= 0, gammainc(x, a) approaches 1 as x
             approaches infinity.

             For a and x real and both a and x small, gammainc(x, a)
             approximates x^a, thus gammainc(0, 0) == 1.

             Both x and a may be complex.

    See Also:
             Erf
             Erfc
             Gamma
             Gammap
             Gammaq
#endif


/* incomplete gamma function */
ITERATE gammainc(x, a, tail)
{
        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("gammainc - input required");
                        
                        a = 1.0;
                }
                
                tail = 'lower';
        }
        
        if (tail == 'upper')
        {
                return(gammaq(a, x));
        }
        else
        {
                return(gammap(a, x));
        }
}