
% ------------------------------------------------------------------------
% Script: Solstrengthen
% ------------------------------------------------------------------------
% Aim: 
% To start the calculation of solid solution strengthening
% ------------------------------------------------------------------------
% Note:
% 1/ The input parameters eb_input and es_input are in unit of 1 
% ------------------------------------------------------------------------

%% the vector of X and Y that have added the misfit

    global pressfld_latxz_mis pressfld_latyy_mis;

    pressfld_latxz_mis=pressfld_latxz;
    pressfld_latyy_mis=pressfld_latyy;
    
    if mistype=='e'
    
        for ii=1:size(pressfld_latxz,1)
            for jj=1:size(pressfld_latxz,2)
                if (pressfld_latyy(ii,jj) > 0.0000)
                    pressfld_latxz_mis(ii,jj)=pressfld_latxz(ii,jj)+misfit_x(pressfld_latxz(ii,jj));
                end
            end
        end    
    
    elseif mistype=='s'
    
        for ii=1:size(pressfld_latxz,1)
            for jj=1:size(pressfld_latxz,2)
                if (pressfld_latyy(ii,jj) > 0.0000)
                    pressfld_latxz_mis(ii,jj)=pressfld_latxz(ii,jj)+misfit_z(pressfld_latxz(ii,jj));
                end
            end
        end 
    
    end

%% The position-dependent dislocation-solute interaction energy

    if strcmp(Einteraction_Plot,'TRUE')

        Einact_plot=zeros(size(pressfld_latxz_mis,1),size(pressfld_latxz_mis,2),size(eb_input,1),size(eb_input,2));
        
        for II=1:size(eb_input,1)
            for JJ=1:size(eb_input,2)

            for ii=1:size(pressfld_latxz_mis,1)
                for jj=1:size(pressfld_latxz_mis,2)

                    Einact_plot(ii,jj,II,JJ) = Einact( pressfld_latxz_mis(ii,jj), pressfld_latyy_mis(ii,jj), eb_input(II,JJ), es_input(II,JJ) );
                    % Unit from 10-20 J to eV
                    Einact_plot(ii,jj,II,JJ)=Einact_plot(ii,jj,II,JJ)*1E-1/1.6000;

                end
            end

            % To plot
            
            Einact_plot_temp=Einact_plot(:,:,II,JJ);
            scatter(pressfld_latxz_mis(:),pressfld_latyy_mis(:),100,Einact_plot_temp(:));

            title('The position-dependent solute/dislocation interaction energy (eV)');
            xlabel('Pressfld_Latxz (A)');
            ylabel('Pressfld_Latyy (A)');

            set(gca,'LineWidth',2)
            axis equal

            % % % Create colorbar

            colorbar;
            colormap('jet');

            % % % To save the figure

            figurefile=strcat(filepath,'/Einteraction_',num2str(eb_input(II,JJ)),'_',num2str(es_input(II,JJ)),'.fig');
            saveas(1,figurefile);
            close(figure(1));

            end
        end

        fprintf('    >> The figure of position-dependent solute/dislocation interaction energy is ploted and saved as "Einteraction.fig" \n');

    end    
    
%% To calculate the solid solution strengthening

    Ebarrier=zeros(size(eb_input,1),size(eb_input,2));
    crss0=zeros(size(eb_input,1),size(eb_input,2));
    wc=zeros(size(eb_input,1),size(eb_input,2));

    fprintf('    >> The results of solid solution strengthening according to the Nabarro-Labusch-Leyson model \n');
    fprintf('\n');
    fprintf('        >> <volume misfit>  <slip misfit>  <bow-out distance (A)>  <Etot/L_dislocation (eV/A)>  <energy barrier at 0K (eV)>  <CRSS at 0K (MPa)> \n');

    for ii=1:size(eb_input,1)
        for jj=1:size(eb_input,2)

            % % zero-temperature energy barrier( ev ) and corresponding CRSS ( MPa )
            [crss0(ii,jj),Ebarrier(ii,jj), wc(ii,jj)]=CRSS(eb_input(ii,jj),es_input(ii,jj));                            % % Etotal(wc(ii,jj))/1600: 1E-10J --> eV/A
            fprintf('           %10.4f %15.4f %18.4f %25.4f %28.4f %26.4f \n', eb_input(ii,jj), es_input(ii,jj), wc(ii,jj), Etotal(wc(ii,jj))/1600, Ebarrier(ii,jj), crss0(ii,jj));

        end
    end
