
% ------------------------------------------------------------------------
% Script: pnstress
% ------------------------------------------------------------------------
% Aim:
% To solve the Peierls stress via discrete dislocation energy approach
% ------------------------------------------------------------------------
% Ref:
% 1/ Joos, B., and M. S. Duesbery. Physical Review Letters 78.2 (1997): 266.
% ------------------------------------------------------------------------

clear u upn
clear Epn_temp Epn strpn 

% space of u

space=0.01;

%% >>>>>> enter into main loop <<<<<<

    ii=1;

    % define some matrix
    Epn_temp=zeros(1,floor(Inpas*2/space)+1);
    strpn=zeros(1,floor(Inpas*2/space));

    % solve the P-N Energy(unit: 1E-10 J/m)

    for u=-Inpas:space:Inpas

        Epn_temp(ii)=0;

        for jj=-2000:2000

        Epn_temp(ii)=Epn_temp(ii)+Gamma(jj*Inpas-u)*Inpas;

        end

        ii=ii+1;

    end

    Epn=(Epn_temp(2:(floor(Inpas*2/space)+1))+Epn_temp(1:floor(Inpas*2/space)))/2;
    
    % solve the P-N Stress(unit: MPa)

    for ii=2:(floor(Inpas*2/space)+1)

        strpn(ii-1)=1/BurVect*(Epn_temp(ii)-Epn_temp(ii-1))/space*1e4;

    end

    upn=-(Inpas-0.005):0.01:(Inpas-0.005);

%% To plot

    [hAx,hLine1,hLine2]=plotyy(upn,Epn,upn,strpn);
    title('Energy Barrier and stress vs. position of dislocation upn');
    xlabel('Position of Dislocation upn (A)');
    ylabel(hAx(1), 'Energy Barrier (1E-10 J/m)'); % left y-axis
    ylabel(hAx(2), 'Stress (MPa)'); % right y-axis
    hLine1.LineStyle='--';
    hLine2.LineStyle=':';

    % To save the figure

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

    % To get the Peierls stress
    
    % % The Peierls energy barrier (WP) is defined as max(W(u))-min(W(u))

    fprintf('\n');
    fprintf('        >> The figure of energy barrier and stress vs. position of Dislocation upn is ploted and saved as "pnstress.fig" \n');
    fprintf('        >> The Peierls barrier and stress are %10.4f (1E-10 J/m) and %10.4f (MPa)\n',max(Epn)-min(Epn),max(abs(min(strpn)),max(strpn)));
    