% XYFIT.M - finds an optimised X-Y distribution
%
%    The program XYFIT.M fits an optimised X-Y distribution to the
%    emittance data taken with the ion source.  The fit is optimised based
%    upon:  1) The assumption that the beam boundary should be an ellipse;
%           2) The integrated X and Y beam profiles
%
%    The X-Y profile is initially assumed to be an ellipse with a uniform
%    particle distribution within that ellipse.  The input matrix for the
%    fitting function is a single matrix that represents histogram data of
%    the particle distribution within the ion source beam.  This matrix is
%    then summed in both X and Y directions to give integrated profiles of
%    the X and Y planes, which are compared to the real data; the quality
%    of this fit is based upon a least-squares method.
%
%    XYFIT uses the Matlab routine FMINSEARCH to minimise a chi-squared fit
%    to the data.  This chi-squared is minimised by varying the input
%    matrix ie. varying the X-Y particle distribution within the simulated
%    beam.  FMINSEARCH in turn uses the function XYCHISQUARE to produce the
%    chi-squared value based upon the X and Y profiles.
%
%    The data is loaded from V:\NLCTA\FONT\phaselockdata\phaselockdata.mat;
%    the actual dataset to be fitted must be specified within the program
%    on the line 'select data to be fitted'.  SHORTFITTER.M uses the Matlab
%    routine FMINSEARCH to minimise a chi-squared fit to the data.  The type
%    of fit is specified both with the 'type' variable and is passed to the
%    function that FMINSEARCH calls to calculate the chi-squared: CHISQUAREFUN.
%    This function will either fit an RC-filtered Gaussian or Square pulse to the
%    data using the 4 specified variables (which are optimised by FMINSEARCH):
%    RC (RC time constant), sigma (standard deviation), tzero (position of
%    peak of fitted function) and A (amplitude scale factor).
%
%    SHORTFITTER.M loads all data into the workspace, selects the data to be
%    used for the fit, calls FMINSEARCH (and CHISQUAREFUN) to minimise the
%    chi-squared, displays and plots the data and does a quick integration
%    (if the fitted function is a square pulse) to check the charge summation.
%

    clear all ; close all ;

%    Load data into workspace

    load E:\df-mat\dfhistdat.mat ;

%    Calculate sums and set up initial fitting matrix

    histpart = 10000 ;

    custmap = jet(64) ;
    custmap(1,:) = [1 1 1] ;

    xsum = sum(xhistdat) ; xpsum = sum(xhistdat') ;
    ysum = sum(yhistdat) ; ypsum = sum(yhistdat') ;
    xyposhist = ysum'*xsum ;
    xprof = xsum./max(xsum) ; yprof = ysum./max(ysum) ;
    xyprof = xyposhist./max(max(xyposhist)) ;
    xsumm = sum(xyprof) ; ysumm = sum(xyprof') ;
    xscal = max(xsumm) ; yscal = max(ysumm) ;
    fitmat = xyprof ;
    sqmat = ( xyprof > 0 ).*1.0 ;
    zermat = zeros(size(sqmat)) ;

%    Call 'fminsearch' to minimise 'xychisquare' function

%    [fittedmat,chi2] = fminsearch(@xychisquare,fitmat,[],xprof,yprof,xscal,yscal) ;
    [fittedmat,chi2] = fminsearchbnd(@xychisquare,fitmat,zermat,sqmat,[],xprof,yprof,xscal,yscal) ;

    break

%    Set fitted variables and display them

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

    disp(' ') ;
    disp(['   Fit type = ' num2str(type)]) ;
    disp(' ') ;
    disp(['   chi-squared = ' num2str(chi2)]) ;
    disp(' ') ;
    disp(['   Fitted values: ']) ;
    disp(' ') ;
    disp(['   RC = ' num2str(RC) ' ; sigma = ' num2str(sigma)]) ;
    disp(['   tzero = ' num2str(tzero) ' ; A = ' num2str(A)]) ;
    disp(' ') ;

%    Plot BPM data and display fitted variables in title

    if type == 1
        titletype = char(['RC-Gaussian']) ;
    elseif type == 2
        titletype = char(['RC-Square']) ;
    end

    figure(100) ; hold on ;
    plot(time,bpm,'-r'), xlabel('Time (ns)'), ylabel('Output (V)'), ...
        title(['\bfShort Pulse X-band data (red) and ' titletype ' fit (black) \newline\chi^{2} = ' num2str(chi2) ...
            ';   \tau = ' num2str(RC) 'ns;   \sigma = ' num2str(sigma) 'ns']) ;

%    Calculate and plot fitted function using fitted variables

    tmin = time(1) ; tmax = time(length(time)) ; tdif = ( tmax - tmin ) / 1000 ;

    if type == 1

        [filtarr,t] = gaussfilt(tmin,tmax,tdif,RC,sigma,tzero,A,1,100) ;

%    If type=2 (i.e. square pulse) integrate function and display summed
%    voltage and predicted integrated voltage

    elseif type == 2

        [filtarr,t] = squarefilt(tmin,tmax,tdif,RC,sigma,tzero,A,1,100) ;

%    Calculate and plot integrated function

        [yinteg,xinteg] = integ(t,filtarr,0,50) ;
%        plot(xinteg,yinteg,'-b') ;

        vtot = yinteg(round( ( 45/50 ) * length(xinteg) ) ) ;
        predict = 2 * sigma * A ;

%    Display predicted and actual voltage values

        disp(['   Predicted integrated voltage = ' num2str(predict)]) ;
        disp(' ') ;
        disp(['   Actual integrated voltage = ' num2str(vtot)]) ;
        disp(' ') ;

    end

    hold off ;

    clear time bpm tmin tmax tdif fitvals type xinteg yinteg vtot predict ;