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


#if @HELP_SININT

    SININT

    Purpose: Evaluates Si(z), the sine integral.

    Syntax:  SININT(z)

              z - A scalar or series.

    Returns: A scalar or series, the value of Si(z), the integration of
             sin(t) / t from 0 to z.

    Example:
             sinint(1)

             returns 0.946083

    Example:
             sinint(1 + i)

             returns 1.104223 + 0.882454i

    Example:
             W1: -10..0.01..10
             W2: sinint(w1);xlabel("x");ylabel("Si(x)");label("Sine Integral");

             displays 2001 samples of Si(x) with integration limits from
             -10 to 10.

    Example:
             W1: 0.01..0.01..10
             W2: sinint(W1)
             W3: integ(sin(w1) / w1)

             W2 contains the sine integral and W3 computes an approximation
             by directly integrating sin(t) / t.

    Remarks:
             The sine integral, Si(z) is the integration of
             sin(t) / t from 0 to z. The input z may be complex.

             Si(z) = (Ei(i * z) - Ei(-i * z)) / (2i) - pi / 2

             for real(z) > 0 and Ei is the Cauchy principal value
             exponential integral implemented by EXPINTEI.

    See Also:
             Cosint
             Expint
             Expintei
             Frenselc
             Frensels
#endif


/* sine integral for real or complex values */
sinint(z)
{
        local y;

        y = iscomplex(z) ? sinint_complex(z) : sinint_real(z);

        return(y);
}


/* purely complex sinint */
sinint_complex(x)
{
        local y, z, neg;

        /* make a copy */
        z = x;

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

        /* complex formulation of x = abs(x) */
        z[neg] = -1.0 * conj(z[neg]);

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

        /* Si(-z) = -Si(z) */
        y[neg] = -1.0 * conj(y[neg]);

        /* Si(0) = 0 */
        y[find(z == 0)] = 0;

        return(y);
}


/* purely real sinint */
sinint_real(x)
{
        local y;

        /* real sine integral for x > 0 */
        y = imag(expint(i * abs(x))) + pi / 2;

        /* Si(-x) = -Si(x) */
        y[find(x < 0)] *= -1;

        /* Si(0) = 0 */
        y[find(x == 0)] = 0;

        return(y);
}