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

#include 

#if @HELP_FRESNELC

    FRESNELC

    Purpose: Evaluates the Fresnel Cosine Integral

    Syntax:  FRESNELC(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 C(x), the integration of
             cos(pi/2 * t^2) from 0 to x.

             (fc, fg) = FRESNELC(x, xc) returns both the integral and
             derivative.

    Example:
             fresnelc(1)

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

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


             returns {0.099998, 0.199921, 0.779893, 0.488253}
             the values of the integral with limits 0.1, 0.2, 1.0 and 2.0.

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

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

    Remarks:
             FRESNELC approximates C(x), the integration of
             cos(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 FRESNELS to evaluate the sine form of the Fresnel
             integral.

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

  See Also:
             ERF
             FRESNELS
             GAMMA
             GAMMALN

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

                 Section 7.3 p300

#endif


/* approximate cosine Fresnel integral */
fresnelc(x, xc)
{
        local sgnx, f1, f2, s, ss, is, fc, ff, fg, sx, cx, 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];

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

        // Approximations for x < 1.6 accurate to eps
        fc[f1] = x[f1] * (1 + (-.1 + (.46296296296296296296e-2 + (-.10683760683760683761e-3 + (.14589169000933706816e-5 + (-.13122532963802805073e-7 + (.83507027951472395917e-10 + (-.39554295164585257634e-12 + (.14483264643598137265e-14 + (-.42214072888070882330e-17 + (.10025164934907719167e-19 + (-.19770647538779051748e-22 + (.32892603491757517328e-25 + (-.46784835155184857737e-28 + .57541916439821717722e-31 * 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);

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

        fc = fc * sgnx;

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

        return(fc);
}