function [filtarr,t] = gaussfilt(tmin,tmax,tdif,RC,sigma,tzero,A,plotstat,plotnum)
%
% GAUSSFILT.M - calculate an RC-filtered Gaussian curve
%
%    The function GAUSSFILT.M calculates a Gaussian distribution
%    that is skewed to the right - the amount of skew depends upon
%    a low pass-style filter with a time constant RC. The idea is
%    replicate the output of the X-band BPM in response to a ns-length
%    beam - the skew corresponds to the response of the BPM. All times
%    are assumed to be in nanoseconds.
%
%    GAUSSFILT.M requires a large number of inputs in the following order:
%
%    function [filtarr,t] = gaussfilt(tmin,tmax,tdif,RC,sigma,tzero,A,plotstat,plotnum)
%
%    tmin, tmax & tdif - time range, in nanoseconds (minimum, maximum and
%                        point-to-point difference) over which the function
%                        is to be evaluated.
%    RC - RC time constant (a.k.a. tau) for the response of the filter. Default
%         is 2.
%    sigma - standard deviation of Gaussian distribution. Default is 1.
%    tzero - approximate time at which peak of function occurs; corresponds
%            to mean position of Gaussian distribution. Default is 0.
%    A - amplitude scaling factor for distribution peak (peak height depends
%        upon A and sigma).
%    plotstat - 1 if the data is to be plotted, 0 if it is not.
%    plotnum - the figure number on which the data is to be plotted.
%

    t = [tmin:tdif:tmax] ;

    for i = 1:length(t)

        filtarr(i) = ( A * ( sqrt(pi/2) ) * sigma * exp( -0.5*( 2*t(i)*RC - 2*tzero*RC - sigma^2 ) / RC^2 ) ) * ...
            ( ( erf( (sqrt(2) * ( t(i)*RC - tzero*RC - sigma^2 ) ) / (2 * sigma * RC) ) ) + ...
            ( erf( (sqrt(2) * ( 30*RC + tzero*RC + sigma^2 ) ) / (2 * sigma * RC) ) ) ) ;

    end

    if plotstat == 1
        figure(plotnum) ; hold on ; plot(t,filtarr,'-k') ;
    end

    return