function chi2 = chisquarefun(x,type,time,bpm)
%
% CHISQUAREFUN.M - calculates a chi-squared
%
%    The function CHISQUAREFUN.M calculates a chi-squared fit for the
%    short pulse data taken with the X-band BPM, stored in
%    V:\NLCTA\FONT\phaselockdata (using the data file in the .MAT
%    file 'phaselockdata.mat' - particularly the array 'ydat_12').
%
%    The function to which the data is fitted is either a Gaussian distribution
%    skewed by a low pass filter i.e. with an added tail, or a square pulse.
%    For an example, see the function GAUSSFILT.M.
%
%    CHISQUAREFUN.M requires 4 inputs in the following order:
%
%    function [chi2] = chisquare(RC,sigma,tzero,A)
%
%    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 25.
%    A - amplitude scaling factor for distribution peak (peak height depends
%        upon A and sigma). Default is 0.45.
%

%    Set the type of fit required: 1 = Gaussian; 2 = Top Hat

%    type = 1 ;
%    type = 2 ;

    RC = x(1); sigma = x(2); tzero = x(3); A = x(4) ;

%    load V:\NLCTA\FONT\phaselockdata\phaselockdata.mat ;

%    time = ydat_14(:,1) ; bpm = ydat_14(:,2) ;
    
    if type == 1

        for i = 1:length(time)
	
            filtarr(i) = ( A * ( sqrt(pi/2) ) * sigma * exp( -0.5*( 2*time(i)*RC - 2*tzero*RC - sigma^2 ) / RC^2 ) ) * ...
                ( ( erf( (sqrt(2) * ( time(i)*RC - tzero*RC - sigma^2 ) ) / (2 * sigma * RC) ) ) + ...
                ( erf( (sqrt(2) * ( 30*RC + tzero*RC + sigma^2 ) ) / (2 * sigma * RC) ) ) ) ;
	
        end

    elseif type == 2

        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
        
    end

    chi2 = sum( ( bpm(250:375) - filtarr(250:375)' ).^2 ) ;

    return