
% ------------------------------------------------------------------------
% Function: createFit(xdata, ydata, Lower, Upper, StartPoint)
% ------------------------------------------------------------------------
% Aim:
% To fit the restoring force using the fitting toolbox of MATLAB
% ------------------------------------------------------------------------

function [fitresult, gof] = createFit(xdata, ydata, Lower, Upper, StartPoint)

    global Nmis;
    global mis_1st;

    [xData, yData] = prepareCurveData( xdata, ydata );

    % % misfunc_1st: 0: without 1st; 1: with 1st

    if mis_1st==0

        if Nmis==1

            fitfunc='fitfuncABC( ux, [Xcoef_1], [0], [Xwid_1] )';

        elseif Nmis>=2

            Xcoef='Xcoef_1';
            Xdist='Xdist_1';
            Xdist_N='-Xdist_1';
            Xwid='Xwid_1';

            for mmnn=2:Nmis-1

                XdistI=strrep(', Xdist_mn','mn',num2str(mmnn));
                Xdist_NI=strrep('-Xdist_mn','mn',num2str(mmnn));
                Xdist=strcat(Xdist, XdistI);
                Xdist_N=strcat(Xdist_N, Xdist_NI);

            end

            Xdist=strcat(Xdist, ',', 32, Xdist_N);

            for mmnn=2:Nmis

                XcoefI=strrep(', Xcoef_mn','mn',num2str(mmnn));
                XwidI=strrep(', Xwid_mn','mn',num2str(mmnn));
                Xcoef=strcat(Xcoef, XcoefI);
                Xwid=strcat(Xwid, XwidI);

            end

            fitfunc=strcat('fitfuncABC( ux, [',Xcoef,'], [',Xdist,'], [',Xwid,'] )');

        end

    elseif mis_1st==1

        if Nmis==1

            fitfunc='fitfuncABC_1st( ux, [Xcoef_1], [0], [Xwid_1], [Xalpha_1] )';

        elseif Nmis>=2

            Xcoef='Xcoef_1';
            Xdist='Xdist_1';
            Xdist_N='-Xdist_1';
            Xwid='Xwid_1';
            Xalpha='Xalpha_1';

            for mmnn=2:Nmis-1

                XdistI=strrep(', Xdist_mn','mn',num2str(mmnn));
                Xdist_NI=strrep('-Xdist_mn','mn',num2str(mmnn));
                Xdist=strcat(Xdist, XdistI);
                Xdist_N=strcat(Xdist_N, Xdist_NI);

            end

            Xdist=strcat(Xdist, ',', 32, Xdist_N);

            for mmnn=2:Nmis

                XcoefI=strrep(', Xcoef_mn','mn',num2str(mmnn));
                XwidI=strrep(', Xwid_mn','mn',num2str(mmnn));
                XalphaI=strrep(', Xalpha_mn','mn',num2str(mmnn));
                Xcoef=strcat(Xcoef, XcoefI);
                Xwid=strcat(Xwid, XwidI);
                Xalpha=strcat(Xalpha, XalphaI);

            end

            fitfunc=strcat('fitfuncABC_1st( ux, [',Xcoef,'], [',Xdist,'], [',Xwid,'], [',Xalpha,'] )');

        end

    end

    %% Set up fittype and options.
    
    ft = fittype( fitfunc, 'independent', 'ux', 'dependent', 'y' );
    opts = fitoptions( ft );
    opts.Algorithm = 'Trust-Region';
    % opts.Algorithm = 'Levenberg-Marquardt';
    opts.DiffMinChange = 1e-30;
    opts.Robust = 'Bisquare';
    % opts.Robust = 'LAR';
    opts.TolFun = 1e-30;
    opts.TolX = 1e-30;
    opts.Lower = Lower;
    opts.Upper = Upper;
    opts.StartPoint = StartPoint;
    opts.Display = 'iter';
    opts.MaxFunEvals = 500*(3*Nmis-1);
    opts.MaxIter = 500;

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

    %% Plot fit with data.
    figure( 'Name', '1d PNt Fit' );
    h = plot( fitresult, xData, yData );
    legend( h, 'The input restoring force', 'The fitting elastic resistance', 'Location', 'NorthWest' );
    % Title and Label axes
    title('The result of fitting restoring force with elastic resistance');
    xlabel( 'x (A)' );
    ylabel( 'Elastic resistance and restoring force (GPa)' );
    grid on

end


