
% ------------------------------------------------------------------------
% Function: gammaFit
% ------------------------------------------------------------------------
% Aim:
% To fit the GSFE data
% ------------------------------------------------------------------------

function [fitresult, gof] = gammaFit(X, Y, Z)

    global fitcut;

    if fitcut>=1

        R_0='';
        fitfunc='';
        mmnn=1;

        for m=0:10
            for n=[0,-1,1,-2,2,-3,3,-4,4,-5,5,-6,6,-7,7,-8,8,-9,9,-10,10]

                if ( sqrt(m^2+n^2) <= fitcut )

                    if ( m==0 && n<=0 )
                        continue;
                    end

                    cosfit=subs(subs(strrep('R_mn*cos(mm*2*pi*Ux+nn*2*pi*Uz)','mn',num2str(mmnn)),'mm',m),'nn',n);
                    sinfit=subs(subs(strrep('I_mn*sin(mm*2*pi*Ux+nn*2*pi*Uz)','mn',num2str(mmnn)),'mm',m),'nn',n);

                    R_0=strcat(R_0,strrep('+R_mn','mn',num2str(mmnn)));
                    fitfunc=strcat(fitfunc,'+',char(cosfit),'-',char(sinfit));

                    mmnn=mmnn+1;

                end

            end
        end

        R_0=strcat('-(',R_0,')');
        fitfunc=strcat(R_0,fitfunc);

    elseif fitcut==0

       fitfunc='-3*(R_1+R_2+R_3+2*R_4)+R_1*(cos(2*pi*Ux-2*pi*Uz)+cos(2*pi*Ux+2*pi*Uz)+cos(2*2*pi*Uz))+R_2*(cos(2*pi*Ux-3*2*pi*Uz)+cos(2*pi*Ux+3*2*pi*Uz)+cos(2*2*pi*Ux))+R_3*(cos(2*2*pi*Ux-2*2*pi*Uz)+cos(2*2*pi*Ux+2*2*pi*Uz)+cos(4*2*pi*Uz))+I_1*(sin(2*pi*Ux-2*pi*Uz)-sin(2*pi*Ux+2*pi*Uz)+sin(2*2*pi*Uz))+I_2*(sin(2*2*pi*Ux-2*2*pi*Uz)-sin(2*2*pi*Ux+2*2*pi*Uz)+sin(4*2*pi*Uz))+R_4*(cos(3*2*pi*Ux+2*pi*Uz)+cos(3*2*pi*Ux-2*pi*Uz)+cos(2*2*pi*Ux+4*2*pi*Uz)+cos(2*2*pi*Ux-4*2*pi*Uz)+cos(2*pi*Ux+5*2*pi*Uz)+cos(2*pi*Ux-5*2*pi*Uz))'; 

    elseif fitcut==-1

       fitfunc='-3*(R_1+R_2+R_3)+R_1*(cos(2*pi*Ux-2*pi*Uz)+cos(2*pi*Ux+2*pi*Uz)+cos(2*2*pi*Uz))+R_2*(cos(2*pi*Ux-3*2*pi*Uz)+cos(2*pi*Ux+3*2*pi*Uz)+cos(2*2*pi*Ux))+R_3*(cos(2*2*pi*Ux-2*2*pi*Uz)+cos(2*2*pi*Ux+2*2*pi*Uz)+cos(4*2*pi*Uz))+I_1*(sin(2*pi*Ux-2*pi*Uz)-sin(2*pi*Ux+2*pi*Uz)+sin(2*2*pi*Uz))+I_2*(sin(2*2*pi*Ux-2*2*pi*Uz)-sin(2*2*pi*Ux+2*2*pi*Uz)+sin(4*2*pi*Uz))'; 

    end

    %% Fit: 'Peierls-Nabarro model'.
    [xData, yData, zData] = prepareSurfaceData( X, Y, Z );

    % Set up fittype and options.
    ft = fittype( fitfunc, 'independent', {'Ux', 'Uz'}, 'dependent', 'SFE' );
    opts = fitoptions( ft );
    opts.Algorithm = 'Levenberg-Marquardt';
    opts.DiffMaxChange = 0.01;
    opts.DiffMinChange = 1e-20;
    opts.Display = 'Off';
    opts.Robust = 'LAR';
    opts.TolFun = 1e-12;
    opts.TolX = 1e-12;

    if fitcut>=1
        opts.StartPoint = zeros(1, 2*(mmnn-1));
    elseif fitcut==0
        opts.StartPoint = zeros(1, 6);
    elseif fitcut==-1
        opts.StartPoint = zeros(1, 5);
    end

    % Fit model to data.
    [fitresult, gof] = fit( [xData, yData], zData, ft, opts );

    %% Plot fit with data.
    figure( 'Name', 'The resutl of fitting gamma surface' );
    h = plot( fitresult, [xData, yData], zData );
    legend( h, 'SFE vs. Ux, Uz', 'Input SFE data', 'Location', 'NorthEast' );
    % Title and Label axes
    title('The resutl of fitting gamma surface');
    xlabel( 'Displacement Ux/b' );
    ylabel( 'Displacement Uz/b' );
    zlabel( 'SFE (J/m2)' );
    grid on

end
