/*****************************************************************************
*                                                                            *
*   GAUSSFIT.SPL    Copyright (C) 2009 DSP Development Corporation           *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Fits y(x) = A * exp(-B*(x-C)^2)                            *
*                                                                            *
*   Revisions:     2 Apr 2009  RRR  Creation - from expfit.spl               *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_GAUSSFIT

    GAUSSFIT

    Purpose: Fits y(x) = A * exp(-B*(x-C)^2) using linearization.

    Syntax:  GAUSSFIT(s)

               s - input series or array


    Returns: A series, the fitted exponential curve.

             (fit, coef) = gaussfit(s)

             Returns both the fit and the coefficients as a series.

             (fit, A, B, C) = gaussfit(s)

             Returns the fit as a series and the coefficients as
             separate scalars.


    Example:
             W1: 10 * exp(-0.5 * ((-10..0.01..10) - 2)^2)

             W2: gaussfit(w1);overp(w1, lred)

             Overplots the original data with the calculated gaussian fit.


    Example:
             (fit, coef) = gaussfit(w1)

             fit is the same series as in W2

             coef == {10.0, 0.5, 2}

    Example:
             (fit, A, B, C) = gaussfit(w1)

             fit is the same series as in W2

             A == 10.0
             B == 0.5
             C == 2.0

    Remarks:
             GAUSSFIT fits an exponential curve of the form y = A*e^(B(x-C)^2).
             The fit is accomplished by fitting a polynomial to the following
             equation:

             ln(y) = -B*x^2 + 2*B*C*x - B*C^2 + ln(A)

             If y is not positive, a complex series will result.

    See Also:
             Expfit
             Polyfit
             Powfit
             Trend
#endif


/* fit y(x) = A * exp(B * (x - C)^2) */
ITERATE gaussfit(s)
{
        local coef, result, A, B, C;

        if (argc < 1) error("gaussfit - input series required");

        /* transform */
        s = ln(s);

        /* deal with infinities */
        s[find(isinf(s))] = 0;

        /* fit */
        if (isxy(s))
        {
                coef = polyfit(s, xvals(s), 2, -1);
        }
        else
        {
                coef = polyfit(s, 2, -1);
        }

        B = -coef[3];
        C = coef[2] / (2 * B);
        A = exp(coef[1] + B * C ^ 2);

        result = A * exp(-B * (xvals(s) - C) ^ 2);

        if (isxy(s))
        {
                result = xy(xvals(s), result);
        }

        if (outargc > 1)
        {
                if (outargc > 2)
                {
                        /* return fit and coefficients as scalars */
                        return(result, A, B, C);
                }

                /* return fit and coefficients as a series */
                return(result, {A, B, C});
        }

        /* return fit only */
        return(result);
}