/*****************************************************************************
*                                                                            *
*   FRESNELS.SPL Copyright (C) 2008 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Evaluates the sin form of the Fresnel Integral S(x)         *
*                                                                            *
*   Revisions:   25 Mar 2008  RRR  Creation, from J. N. McElwaine            *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_FRESNELS

    FRESNELS

    Purpose: Evaluates the Fresnel Sine Integral

    Syntax:  FRESNELS(x, xc)

                x - A real or series, the integration limit.

               xc - Optional real, the cutoff limit that determines
                    the integration method. Defaults to 1.8.

    Returns:
             A real or series, the value of S(x), the integration of
             sin(pi/2 * t^2) from 0 to x.

    Example:
             fresnels(1)

             returns 0.438259, the value of the integration of
             sin(pi/2 * t^2) for t = 0.0 to 1.0.

    Example:
             fresnels({0.1, 0.2, 1, 2})

             returns {0.000524, 0.004188, 0.438259, 0.343416},
             the values of the integral with limits 0.1, 0.2, 1.0 and 2.0.

    Example:
             fresnels(-0.5..0.01..0.5)

             returns 1001 samples of the Fresnel Sine Integral from -5 to 5.

    Remarks:
             FRESNELS approximates S(x), the integration of
             sin(pi/2 * t^2) for t = 0 to x.

             For abs(x) < xc, a power series about x = 0 is used
             yielding an accuracy better than 5e-16.

             For abs(x) > xc, a minimax rational approximation based on
             auxilliary functions described in [1] is used yielding an 
             accuracy better than 1e-9.

             If the integral limit x is a series, the result is a series.

             See FRESNELC to evaluate the cosine form of the Fresnel
             integral.

             FRESNELS was developed from an algorithm by J. N. McElwaine.

  See Also:
             ERF
             FRESNELC
             GAMMA
             GAMMALN

  References:
             [1] Abramowitz and Stegun
                 Handbook of Mathematical Functions (9th printing 1970)
                 US Gov. Printing Office
                 Section 7.3 p300

#endif


/* approximate Fresnel sine integral */
fresnels(x, xc)
{
        local sgnx, f1, f2, s, ss, is, fs, ff, fg, cx, sx, scalar = FALSE;

        // Cut point for switching between power series and rational approximation
        if (argc < 2)
        {
                if (argc < 1)
                {
                        x = 1.0;
                }
                
                xc = 1.8;
        }

        if (isscalar(x))
        {
                scalar = TRUE;
                x      = {x};
        }

        sgnx = sign(x);
        x    = abs(x);

        f1 = find(x <= xc);
        f2 = find(x >  xc);

        s = pi / 2 * x[f1] ^ 2;

        ss = s ^ 2;
        is = sqrt(2 / pi) / x[f2];


        fs = zeros(length(x), 1);
        ff = fs;
        fg = fs;

        // Approximations for x < 1.6 accurate to eps

        fs[f1] = x[f1] * s * (1 / 3 + (-.23809523809523809524e-1 + (.75757575757575757576e-3 + (-.13227513227513227513e-4 + (.14503852223150468764e-6 + (-.10892221037148573380e-8 + (.59477940136376350368e-11 + (-.24668270102644569277e-13 + (.80327350124157736091e-16 + (-.21078551914421358249e-18 + (.45518467589282002862e-21 + (-.82301492992142213568e-24 + (.12641078988989163522e-26 + (-.16697617934173720270e-29 + .19169428621097825308e-32 * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss) * ss);
        ff[f2] = sqrt(2 / pi) * is * (0.36634901025764670680 + (0.84695666222194589182 + (-3.7301273487349839902 + (10.520237918558456228 + (-10.472617402497801126 + 4.4169634834911107719 * is) * is) * is) * is) * is) / (.73269802661202980832 + (1.6939102288852613619 + (-7.4599994789665215344 + (21.032436583848862358 + (-20.269535575486282590 + 8.9995024877628789836 * is) * is) * is) * is) * is);
        fg[f2] = sqrt(2 / pi) * is ^ 3 * (.10935692320079194659 + (.72025533055541994934 + (-2.5630322339412334317 + (7.2404268469720856874 + (-7.0473933910697823445 + 2.9315207450903789503 * is) * is) * is) * is) * is) / (.43742772185339366619 + (2.8810049959848098312 + (-10.250672312139277005 + (28.912791022417600679 + (-25.740131167525284201 + 15.145134363709932380 * is) * is) * is) * is) * is);

        s  = pi / 2 * x[f2] ^ 2;
        cx = cos(s);
        sx = sin(s);

        fs[f2] = 1 / 2 - ff[f2] * cx - fg[f2] * sx;

        fs = fs * sgnx;

        if (scalar)
        {
                fs = castreal(fs);
        }
        else
        {
                setdeltax(fs, deltax(x));
                setxoffset(fs, xoffset(x));
        }

        return(fs);
}