
% ------------------------------------------------------------------------
% Script: 1d_lsq_curvefit_1st
% ------------------------------------------------------------------------
% Aim:
% 1/ To solve 1D P-N model via the lsqcurvefit function of MATLAB with considering 
%    the first-order approximation in the trial function of disregistry vector
% ------------------------------------------------------------------------

global Kx;
global Xcoef Xdist Xwid Xalpha;

%   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), Xalpha_range(1)*ones(1,1)];
        opts.Upper = [Xcoef_range(2)*ones(1,1), Xwid_range(2)*ones(1,1), Xalpha_range(2)*ones(1,1)];
        coef = [Xcoef_range(3)*ones(1,1), Xwid_range(3)*ones(1,1), Xalpha_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), Xalpha_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), Xalpha_range(2)*ones(1,Nmis)];
        coef = [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), Xalpha_range(3)*ones(1,Nmis)];

    end

    % options=optimoptions('lsqcurvefit','MaxIterations',500,'Display','iter','FunctionTolerance',1E-15);
    % options=optimoptions(options,'MaxFunctionEvaluations',500*(3*Nmis-1),'StepTolerance',1E-8);
    % options=optimoptions(options,'PlotFcn', {[]},'UseParallel',false);

    options=optimoptions('lsqcurvefit','MaxIter',500,'Display','iter','TolFun',1E-15);
    options=optimoptions(options,'MaxFunEvals',500*(3*Nmis-1)); 

    lsq_result=evalc('[coef,resnorm]=lsqcurvefit(@fitfuncABCx_1st, coef, Uxdata, RestoringForce, opts.Lower, opts.Upper, options);');
    fprintf('%s',strrep(['        ',lsq_result],char(10),[10,'        ']));

    % R-square

    rsquare=1-resnorm/sum((RestoringForce-mean(RestoringForce)).^2);

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

    % To get the figure of fitrestoringforce.fig

    RestoringForcefit=fitfuncABCx_1st(coef,Uxdata);
    plotxy = plot(Uxdata,RestoringForce,Uxdata,RestoringForcefit);
    set(plotxy(1),'MarkerSize',8,'Marker','o','LineStyle','none','Color','b');
    set(plotxy(2),'Color','r');

    legend( '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

    % To save the figure

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

    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);
    Xalpha=zeros(1,Nmis);

    % % for Xcoef, Xwid and Xalpha

    for ij=1:Nmis

        Xcoef(ij)=coef(ij);
        Xwid(ij)=coef(2*Nmis+ij-1);
        Xalpha(ij)=coef(3*Nmis+ij-1);

    end

    Xcoef=Xcoef/sum(Xcoef);

    % % for Xdist

    if Nmis==1

        Xdist(1)=0;

    elseif Nmis>=2

        for ij=1:Nmis-1
            Xdist(ij)=coef(Nmis+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   Xalpha \n');

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