
% ------------------------------------------------------------------------
% Script: pnstress_analysis
% ------------------------------------------------------------------------
% Aim:
% To solve the Peierls stress via analytical formula
% ------------------------------------------------------------------------
% Ref:
% 1/ Joos, B., and M. S. Duesbery. Physical Review Letters 78.2 (1997): 266.
% ------------------------------------------------------------------------

%   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);

% the ideal slip stress 

    % To fit the gamma surface
    
    if ~exist('gammafitresult','var') || isempty(gammafitresult)
        
        % % To fit the gamma surface
        
        fprintf('\n');
        fprintf('        >> To start fitting the gamma surface \n');

        [gammafunc,gof]=gammaFit(Ux,Uz,SFE);

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

        % % To save the figure "fitgamma.fig"

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

        fprintf('            >> The figure of SFE vs displacement ux and uz is ploted and saved as "fitgamma.fig" \n');
        
        % % The gamma fitting function, which has high efficiency

        temp_AA=coeffnames(gammafunc);
        temp_BB=coeffvalues(gammafunc);

        temp_gfit=formula(gammafunc);
        for ii=1:numcoeffs(gammafunc)
            temp_gfit=subs(temp_gfit,temp_AA{ii},temp_BB(ii));
        end
        temp_gfit=vpa(temp_gfit,4);

        % % % To define the function 'gammafitresult'
        eval(strcat('gammafitresult=@(Ux,Uz)',char(temp_gfit),';'));
        
    end
    
    % % (GPa)
    
    dsfe_x=-(diff(gammafitresult(0:0.0001:1, 0.00))/0.0001)*10/BurVect;
    dsfe_x_fliplr=fliplr(dsfe_x);
    
    str_is_positive=dsfe_x(1);
    str_is_negative=dsfe_x_fliplr(1);
    
    for ij=1:size(dsfe_x,2)
        if dsfe_x(ij+1) < dsfe_x(ij)
        	str_is_positive=dsfe_x(ij+1);
        else
            break;
        end
    end

    for ij=1:size(dsfe_x_fliplr,2)
        if dsfe_x_fliplr(ij+1) > dsfe_x_fliplr(ij)
        	str_is_negative=dsfe_x_fliplr(ij+1);
        else
            break;
        end
    end
    
    str_is_x=min(abs(str_is_positive),abs(str_is_negative));

% the dislocation halp width (Angstrom)

half_width=Kx*BurVect/(4*pi*str_is_x);

% For narrow dislocation (i.e. half_with/Inpas <<1) (MPa)

str_pn_narrow=1E3*3.00*sqrt(3)/8.00*str_is_x*Inpas/(pi*half_width);

% For wide dislocation (i.e. half_with/Inpas >>1) (MPa)

str_pn_wide=1E3*Kx*BurVect/Inpas*exp(-2*pi*half_width/Inpas);

% To output the result of Peierls stress

fprintf('\n');
fprintf('        >> The ideal slide stress is %10.4f (GPa) \n', str_is_x);
fprintf('        >> The dislocation half width is %10.4f (A) \n', half_width);
fprintf('        >> The dimensionless quantity half_width/a'' is %10.4f \n', half_width/Inpas);
fprintf('        >> The Peierls stress with narrow dislocaiton (half_width/a''<<1) is %10.4f (MPa) \n', str_pn_narrow);
fprintf('        >>                         wide dislocaiton (half_width/a''>>1) is %10.4f (MPa) \n', str_pn_wide);
    
