
% ------------------------------------------------------------------------
% Script: 1d_fitting_toolbox
% ------------------------------------------------------------------------
% Aim:
% To solve 1D P-N model via the fitting toolbox of MATLAB
% ------------------------------------------------------------------------

global Kx;
global Xcoef Xdist Xwid;

%   edge or screw dislocation

    % "theta" is the angle between the dislocation line and its Burgers 
    %         vector for a general mixed dislocation in an isotropic solid

    if mistype=='e'   % for edge dislocation
        theta=pi/2;
    end
    
    if mistype=='s'   % for screw dislocation
        theta=0;
    end
    
    Kx=shear_modulus*(sin(theta)^2/(1-poisson_ratio)+cos(theta)^2);
    
%% To fit the restoring force 

    if Nmis==1

        opts.Lower = [Xcoef_range(1)*ones(1,1), Xwid_range(1)*ones(1,1)];
        opts.Upper = [Xcoef_range(2)*ones(1,1), Xwid_range(2)*ones(1,1)];
        opts.StartPoint = [Xcoef_range(3)*ones(1,1), Xwid_range(3)*ones(1,1)];

    elseif Nmis>=2

        opts.Lower = [Xcoef_range(1)*ones(1,Nmis), Xdist_range(1)*ones(1,Nmis-1), Xwid_range(1)*ones(1,Nmis)];
        opts.Upper = [Xcoef_range(2)*ones(1,Nmis), Xdist_range(2)*ones(1,Nmis-1), Xwid_range(2)*ones(1,Nmis)];
        opts.StartPoint = [Xcoef_range(3)*ones(1,Nmis), Xdist_range(3)*[ones(1,ceil((Nmis-1)/2.0)),-ones(1,(Nmis-1)-ceil((Nmis-1)/2.0))], Xwid_range(3)*ones(1,Nmis)];

    end

    toolbox_result=evalc('[RestoringForcefitresult, gof] = createFit(Uxdata, RestoringForce, opts.Lower, opts.Upper, opts.StartPoint);');

    fprintf('%s',strrep(['        ',toolbox_result],char(10),[10,'        ']));

    fprintf('\n');
    if gof.rsquare >= 0.9800
        fprintf('        >> The rsquare of fitting the restoring force is %10.4f \n', gof.rsquare);
    else
        fprintf('        >> WARNING!!!  The rsquare of fitting the restoring force is < 0.9800 \n');
    end

    % To output the fit values yyfit
    
    RestoringForcefit=RestoringForcefitresult(Uxdata);
    
    % To save the figure

    figurefile=strcat(filepath,'/fitrestoringforce.fig');
    saveas(1,figurefile);
    close(figure(1));

    fprintf('\n');
    fprintf('    >> The results of fitting restoring force with elastic resistance is ploted and saved as "fitrestoringforce.fig \n');

%% To get the results of variational constants

    % % Define some matrix

    Xcoef=zeros(1,Nmis);
    Xdist=zeros(1,Nmis);
    Xwid=zeros(1,Nmis);

    % % for Xcoef and Xwid

    for ij=1:Nmis

        Xcoef(ij)=eval(strrep('RestoringForcefitresult.Xcoef_mn','mn',num2str(ij)));
        Xwid(ij)=eval(strrep('RestoringForcefitresult.Xwid_mn','mn',num2str(ij)));

    end

    Xcoef=Xcoef/sum(Xcoef);

    % % for Xdist

    if Nmis==1

        Xdist(1)=0;

    elseif Nmis>=2

        for ij=1:Nmis-1
            Xdist(ij)=eval(strrep('RestoringForcefitresult.Xdist_mn','mn',num2str(ij)));
        end

        Xdist(Nmis)=-sum(Xdist(1:Nmis-1));

    end

%% To output the results of variational constants in disregistry function

    fprintf('    >> The results of variational constants in disregistry function (Angstrom) \n');
    fprintf('\n');
    fprintf('        >>   NO.   Xcoef   Xdist   Xwid \n');

    for ij=1:Nmis
        fprintf('           %3i %10.4f %10.4f %10.4f \n',ij, Xcoef(ij), Xdist(ij), Xwid(ij));
    end
    fprintf('\n');
