/*****************************************************************************
*                                                                            *
*   COSINT.SPL Copyright (C) 2014 DSP Development Corporation                *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:    Randy Race                                                    *
*                                                                            *
*   Synopsis:  Computes the cosine integral                                  *
*                                                                            *
*   Revisions: 15 May 2014  RRR  Creation                                    *
*                                                                            *
*****************************************************************************/


#if @HELP_COSINT

    COSINT

    Purpose: Evaluates Ci(z), the cosine integral.

    Syntax:  COSINT(z)

              z - A scalar or series.

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

    Example:
             cosint(1)

             returns 0.337404

    Example:
             cosint(1 + 1i)

             returns 0.882172 + 0.287249i

    Example:
             W1: 0..0.01..10
             W2: cosint(w1);xlabel("x");ylabel("Ci(x)");label("Cosine Integral");

             displays 1001 samples of Ci(x) with integration limits from
             0 to 10.

    Example:
             W1: 01..0.01..10
             W2: cosint(W1)
             W3: gamm + ln(w1) + integ((cos(w1) - 1) / w1)

             W2 contains the cosine integral and W3 computes an approximation
             by directly integrating (cos(t) - 1) / t

    Remarks:
             The cosine integral, Ci(z) is the integration of
             cos(t) / t from z to inf. The input z may be complex.

             The function can also be computed as:

             gamm + ln(z) + integ((cos(z) - 1) / z)

             from 0 to z where gamm is the Euler constant.

             Ci(z) = -(E1(i * z) + E1(-i * z)) / 2

             for real(z) > 0 and E1 is the exponential integral implemented
             by EXPINT.

    See Also:
             Expint
             Expintei
             Frenselc
             Frensels
             Sinint
#endif


/* cosine integral for real or complex values */
cosint(z)
{
        local y;

        y = iscomplex(z) ? cosint_complex(z) : cosint_real(z);

        return(y);
}


/* purely complex cosint */
cosint_complex(z)
{
        local y, np, nn, zr;

        /* negative values */
        np = find(real(z) < 0 && imag(z) > 0);
        nn = find(real(z) < 0 && imag(z) < 0);

        /* zeros */
        zr = find(real(z) == 0);

        /* complex cosine integral for real(z) > 0 */
        y = -0.5 * (expint(i * z) + expint(-i * z));

        /* Ci(-z) = -Ci(z) - i * pi */
        y[np] += i * pi;
        y[nn] -= i * pi;
        y[zr] =  real(y[zr]) + i * phase(z[zr]);

        return(y);
}


/* purely real cosint */
cosint_real(x)
{
        local y, neg;

        /* negative values */
        neg = find(x < 0);

        /* real cosine integral for x > 0 */
        y = -1.0 * real(expint(i * abs(x)));

        /* Ci(-x) = -Ci(x) - i * pi */
        y[neg] += i * pi;

        return(y);
}