/*****************************************************************************
*                                                                            *
*   CUMTRAPZ.SPL  Copyright (C) 2009-2015 DSP Development Corporation        *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Integration using the trapezoidal rule                      *
*                                                                            *
*   Revisions:   17 Nov 2009  RRR  Creation - from TRAPZ.SPL                 *
*                29 Sep 2015  RRR  difference equation with initial conds    *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_CUMTRAPZ

    CUMTRAPZ

    Purpose: Calculates the integral of a series using the trapezoidal rule.

    Syntax:  CUMTRAPZ(series)

              series - An interval or XY series, the input data.

    Alternative Syntax:  

              CUMTRAPZ(xseries, yseries)

              xseries - A series, the X values.

              yseries - A series, the Y values.

    Returns: A series

    Example:
             W1: {0, 2, 4, 6, 8, 10}
             W2: integ(w1)
             W3: cumtrapz(w1)

             W2 == {0.0, 1.3, 4.0, 9.3, 16.0, 25.3}

             W3 == {0.0, 1.0, 4.0, 9.0, 16.0, 25.0}

             W2 contains the integrated series using Simpson's rule and
             W3 contains the integrated series using the Trapezoidal rule.

    Example:
             y1 = gnorm(1000, .001);
             a1 = cumtrapz(y1);
             a2 = trapz(y1);

             Series a1 contains the cumulative area of y1 and scalar a2
             contains the value of the total area of y1.  Note that
             a1[end] == a2, the last point of the cumulative area
             equals the total area.

    Remarks:
             For series S, the trapezoidal rule calulates the running sum
             of:
                        deltax * (S[i+1] + S[i]) / 2.

    See Also:
             Area
             Integ
             Trapz
#endif


/* trapezoidal rule */
ITERATE cumtrapz(x, y)
{
        local t, dx, vu, b, have_x;

        /*
         *  calculate:
         *
         *  dx * sum( S[i] + (S[i] - S[i-1])/2 )
         *
         *  by using the moving average without end padding for XY series
         *  or a difference equation for interval series;
         */

        if (argc < 2)
        {
                if (argc < 1) error("cumtrapz - input series required");
                
                s      = refseries(x);
                have_x = FALSE;
        }
        else
        {
                s      = refseries(y);
                have_x = TRUE;
        }

        /* vertical units */
        vu = getvunits(integ(extract(s, 1, 2)));

        if (isxy(s) || have_x)
        {
                t = movavg(s, 2, 0);

                /* zero first point */
                t[1] = 0.0;

                /* XY series */
                if (not(have_x)) x = xvals(s);

                /* delta x's */
                dx = extract(x - delay(x, 1), 1, length(s));

                /* sum */
                t = partsum(dx * t);

                /* get correct length */
                t = extract(t, 1, length(s));

                /* units */
                setvunits(t, vu);

                /* return XY series */
                t = xy(x, t);
        }
        else
        {
                /* difference equation with proper initial condition */
                b = 0.5 * deltax(s);
                t = filteq({b, b}, {1, -1}, s, {-s[1]*b});

                /* units */
                setvunits(t, vu);
        }

        return(t);
}