/*****************************************************************************
*                                                                            *
*   EXPINTEI.SPL Copyright (C) 2014 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Computes the Cauchy Principal Value Exponential Integral    *
*                                                                            *
*   Revisions:   28 May 2014  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_EXPINTEI

    EXPINTEI

    Purpose: Evaluates Ei(z), the Cauchy principal value exponential integral.

    Syntax:  EXPINTEI(z)

              z - A scalar or series.

    Returns: A scalar or series, the value of Ei(z), the integration of
             -exp(t) / t from -inf to z.

    Example:
             expintei(1)

             returns 1.895118

    Example:
             expintei(0)

             returns -inf

    Example:
             expintei(1 + 1i)

             returns 1.764626 + 2.387770i

    Example:
             W1: -5..0.01..5
             W2: expintei(W1);xlabel("x");ylabel("Ei(x)");label("Ei Integral");

             displays 1001 samples of Ei(x) with integration limits from
             -5 to 5.

    Example:
             W1: 01..0.01..10
             W2: expintei(W1)
             W3: integ(exp(w1) / w1)

             W2 contains the Cauchy principal value exponential integral and W3
             computes an approximation by directly integrating exp(t) / t.

    Remarks:
             The Cauchy principal value exponential integral, Ei(z), is the
             integration of exp(t) / t from -inf to z.

             The input z may be complex. The function is multi-valued with a
             branch cut along the negative real axis.

             For real positive z:

             E1(z) = -Ei(-z)

             Where E1(z) is the generalized exponential integral implemented
             by EXPINT.

             Ei(z) is real for all real z.

             Ei(0) == -inf.

    See Also:
             Cosint
             Expint
             Frenselc
             Frensels
             Sinint
#endif


/* Cauchy principal value integral Ei(z) */
expintei(z)
{
        local y, pp, pn, np, nn, ipi;

        /* base value */
        y = -expint(-z);

        if (iscomplex(z))
        {
                /* branch points */
                pp = find(real(z) >  0 && imag(z) >  0);
                pn = find(real(z) >= 0 && imag(z) <= 0);
                np = find(real(z) <= 0 && imag(z) >  0);
                nn = find(real(z) <  0 && imag(z) <  0);

                ipi = i * pi;

                y[pp] += ipi;
                y[pn] -= ipi;
                y[np] += ipi;
                y[nn] -= ipi;
        }
        else
        {
                /* always real for real z */
                y = real(y);
        }

        return(y);
}