function [filtarr,time] = squarefilt(tmin,tmax,tdif,RC,sigma,tzero,A,plotstat,plotnum)
%
% SQUAREFILT.M - calculate an RC-filtered Top Hat curve
%
%    The function GAUSSFILT.M calculates a square pulse curve
%    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.
%
%    SQUAREFILT.M requires a large number of inputs in the following order:
%
%    function [filtarr,t] = squarefilt(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 (corresponds to half-width
%    of square pulse). Default is 1.
%    tzero - approximate time at which peak of function occurs; corresponds
%            to centre of square pulse. Default is 0.
%    A - amplitude scaling factor for distribution peak (peak height depends
%        upon A, RC 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.
%

    time = [tmin:tdif:tmax] ;

    starttime = tzero - sigma ;  endtime = tzero + sigma ;

    for i = 1:length(time)

        if time(i) < starttime

            filtarr(i) = 0 ;
            
        elseif time(i) > starttime & time(i) < endtime

            filtarr(i) = ( A * ( 1 - exp( -( time(i) - starttime ) / RC ) ) ) ;
            
        elseif time(i) > endtime

            filtarr(i) = ( A * ( 1 - exp( -( endtime - starttime ) / RC ) ) * ( exp( -( time(i) - endtime ) / RC ) ) ) ;
            
        end

    end

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

    return