% ------------------------------------------------------------------------
% Script: StartCal
% ------------------------------------------------------------------------
% Aim:
% Start to run the code PNADIS
% ------------------------------------------------------------------------

% To turn off the warning

warning('off','all')

% clear all variates before starting

clear
clear global

clc

% close all figures before starting

close all

fprintf('%% %% ---------------------------------------------------------------------- \n');
fprintf('%% %% INPUTTING THE PARAMETER VALUE \n');
fprintf('%% %% ---------------------------------------------------------------------- \n');
fprintf('\n');

% set global variables

global BurVect
global Nmis fitcut;
global Ux Uz SFE SFEfit;
global gammafitresult;
global filepath misdim mistype pnmode;
global shear_modulus poisson_ratio;

global mis_1st;
global fitmethod minimethod;
global dax Inpas;
global StressEnergy GradientEnergy NonlocalEnergy;
% global str_coef grad_coef nonloc_coef;

global pressfld_latxz pressfld_latyy;
global dis_component_xz;
global eb_input es_input;
global solute_concentration LatC_spacing coefV;
global wc_range;

% To input the parameters

infile_path=strcat(pwd,'/infile.m');
run(infile_path);

% For the gradient, nonlocal and stress energy, which are not included in V_1.0.0

    % % whether consider the gradient energy or NOT and the gradient
    %   coefficient only for SVPN model (i.e. pnmode=2)
    
    GradientEnergy=0;
    % grad_coef=[];
    
    % % whether consider the nonlocal energy or NOT and the nonlocal
    %   coefficient only for SVPN model (i.e. pnmode=2) 
    
    NonlocalEnergy=0;
    % nonloc_coef=[];
    
    % % whether consider the stress energy or NOT and the applied stress
    %   only for SVPN model (i.e. pnmode=2)
    
    StressEnergy=0;
    % str_coef=[];

% To get the default value of parameters

run Default.m

if defaultend==1
    fprintf('    >> The inputing of the parameter value is finished \n');
    fprintf('\n');
else
    fprintf('\n');
    fprintf('    >> Error happens during inputing the parameter value \n');
    fprintf('\n');
    return
end

% To creat a new outfile 'outfile.txt'

    % % outfile=strcat(filepath,'/outfile.txt');
    % % FID=fopen(outfile, 'wt');
    % % % FID=fopen(outfile, 'at');
    % % fclose(FID);

fprintf('%% %% ---------------------------------------------------------------------- \n');
fprintf('%% %% STARTING THE PNADIS CODE \n');
fprintf('%% %% ---------------------------------------------------------------------- \n');

%% ------------------------------------------------------------------------
% To calculate the dislocation core structure
% ------------------------------------------------------------------------

if strcmp(Dislocation_Core_Structure,'TRUE')

    fprintf('\n');
    fprintf('%% %% INITIATION OF CALCULATING THE DISLOCATION CORE STRUCTURE \n');
    fprintf('\n');
    
    % ------------------------------------------------------------
    % The 1D Peierls-Nabarro model
    % ------------------------------------------------------------

    if misdim==1

        fprintf('    >> The 1D Peierls-Nabarro Model is choden \n');

        % % To fit the gamma surface

        fprintf('    >> To start fitting the gamma surface \n');

        [Uxdata,RestoringForce]=fitgamma(Ux,SFE);

        % % The First Mothed "1d_fitting_toolbox"

        if strcmp(fitmethod,'1d_fitting_toolbox')

            fprintf('    >> Curve Fitting Toolbox was chosen to fit the restoring force \n');

            if mis_1st==1

                fprintf('    >> First approximation was used in disregistry fucntion \n');

                run fitting_toolbox_1d_1st.m
            elseif mis_1st==0

                fprintf('    >> NO First approximation was used in disregistry fucntion \n');

                run fitting_toolbox_1d.m
            end

        % % The Second Mothed "1d_lsq_curvefit"

        elseif strcmp(fitmethod,'1d_lsq_curvefit')

            fprintf('    >> "lsqcurvefit" function was chosen to fit the restoring force \n');

            if mis_1st==1

                fprintf('    >> First approximation was used in disregistry fucntion \n');

                run lsq_curvefit_1d_1st.m
            elseif mis_1st==0

                fprintf('    >> NO First approximation was used in disregistry fucntion \n');

                run lsq_curvefit_1d.m
            end

        end

        % % To plot

        % % % To get the figure of disregistry and misfit density vs x/b

        xmis=-20*BurVect:0.1:20*BurVect;
        misfit=misfit1d(xmis);
        misdens=misdens1d(xmis);
        xmis=xmis'/BurVect;
        misfit=misfit';
        misdens=misdens';
        misfit=misfit/BurVect;

        [hAx,hLine1,hLine2]=plotyy(xmis,misfit,xmis,misdens);

        title('Disregistry and misfit density vs. x/b');
        xlabel('x/b');
        ylabel(hAx(1), 'Disregistry/b'); % left y-axis
        ylabel(hAx(2), 'Misfit density'); % right y-axis
        hLine1.LineStyle='--';
        hLine2.LineStyle=':';

        % % % To save the figure

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

        fprintf('    >> The figure of disregistry and misfit density vs. x/b is ploted and saved as "misdens.fig" \n');
    
    % ------------------------------------------------------------
    % The 2D Peierls-Nabarro model
    % ------------------------------------------------------------

    elseif misdim==2

        fprintf('    >> The 2D Peierls-Nabarro Model is choden \n');

        % % To fit the gamma surface

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

        % % To plot the gammafiterror.fig
        
            % % % To output the fit values SFEfit
    
            SFEfit=gammafunc(Ux,Uz);
            gammafiterror=SFEfit-SFE;
            
            % % % To plot

            surf(Ux,Uz,0*gammafiterror,gammafiterror);            
            
            % % % % Title and Label axes
            title('The SFE fitting error (J/m2)');
            xlabel( 'Displacement Ux/b' );
            ylabel( 'Displacement Uz/b' );
            zlabel( 'SFE fitting error (J/m2)' );
            grid on
            view(0,90);
            set(gca,'DataAspectRatio',[sqrt(3) 1 1]);
            
            % % % % Create colorbar
            colorbar;
            colormap('jet');

            % % % % To save the figure "gammafiterror.fig"
            figurefile=strcat(filepath,'/gammafiterror.fig');
            saveas(1,figurefile);
            close(figure(1));

            fprintf('        >> The figure of the difference between the fitting and input SFEs is ploted and saved as "gammafiterror.fig" \n');
            fprintf('\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),';'));

        % % The Classical P-N model or Semidiscrete Variational P-N model
        
        if pnmode==1
            
            fprintf('    >> The Classical Peierls-Nabarro model was chosen \n');
            
        elseif pnmode==2
            
            fprintf('    >> The Semidiscrete Variational Peierls-Nabarro model was chosen \n');
            fprintf('\n');

            % if GradientEnergy==0
            %     fprintf('        >> The gradient energy was not considered \n');
            % elseif GradientEnergy==1
            %     fprintf('        >> The gradient energy was considered with grad_coef=%8.4f \n', grad_coef);
            % end
            % 
            % if NonlocalEnergy==0
            %     fprintf('        >> The nonlocal energy was not considered \n');
            % elseif NonlocalEnergy==1
            %     fprintf('        >> The nonlocal energy was considered with nonloc_coef=%8.4f \n', nonloc_coef);
            % end
            % 
            % if StressEnergy==0
            %     fprintf('        >> The stress energy due to the applied stress was not considered \n');
            % elseif StressEnergy==1
            %     fprintf('        >> The stress energy due to the applied stress was considered with str_coef=%8.4f \n', str_coef);
            % end

            % fprintf('\n');          
        end
            
        % % The First Mothed "2d_particle_swarm_optimization"

        if strcmp(minimethod,'2d_PSO')

            fprintf('    >> Particle_Swarm_Optimization (PSO) was chosen to minimize the total dislocation energy \n');
            fprintf('\n');

            run particle_swarm_optimization_2d.m

        % % The Second Mothed "2d_genetic_algorithm"

        elseif strcmp(minimethod,'2d_GA')

            fprintf('    >> Genetic_Algorithm (GA) was chosen to minimize the total dislocation energy \n');
            fprintf('\n');            

            run genetic_algorithm_2d.m

        end

        % % To plot

        % % % To get the figure of disregistry vs x/b

        xmis=-20*BurVect:0.1:20*BurVect;
        misfit_x_axis=misfit_x(xmis);
        misfit_z_axis=misfit_z(xmis);

        xmis=xmis/BurVect;
        misfit_x_axis=misfit_x_axis/BurVect;
        misfit_z_axis=misfit_z_axis/BurVect;
        plot(xmis,misfit_x_axis,xmis,misfit_z_axis);

        title('Disregistry ux and uz vs. x/b');
        xlabel('x/b');
        ylabel('Disregistr/b');

        % % % To save the figure

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

        fprintf('    >> The figure of disregistry vs. x/b is ploted and saved as "misfit.fig" \n');

        % % % To get the figure of misfit density vs x/b

        xmis=-20*BurVect:0.1:20*BurVect;
        misdens_x_axis=misdens_x(xmis);
        misdens_z_axis=misdens_z(xmis);

        xmis=xmis/BurVect;
        plot(xmis,misdens_x_axis,xmis,misdens_z_axis);

        title('Misfit density dux/dx and duz/dx vs. x/b');
        xlabel('x/b');
        ylabel('Misfit density');

        % % % To save the figure

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

        fprintf('    >> The figure of misfit density vs. x/b is ploted and saved as "misdens.fig" \n');
        
    end

    fprintf('\n');
    fprintf('%% %% THE END \n');

elseif strcmp(Dislocation_Core_Structure,'FALSE')

    fprintf('\n');
    fprintf('%% %% NO dislocation core structure was calculated \n'); 
    fprintf('\n');
    
    if (strcmp(Peierls_Stress,'TRUE') && pnstrmethod==2) || strcmp(Pressure_Field,'TRUE') || ( misdim==2 && strcmp(Solution_Strengthening,'TRUE'))
    
        if exist('pnadis.mat','file')
            clear gammafitresult;
            save('pnadis.mat','*','-append');
            load('pnadis.mat');

            fprintf('    >> Load the data file "pnadis.mat" \n');

        else

            fprintf('\n');
            fprintf('%% %% ERROR!!! NO pnadis.mat is obtained for Dislocation_Core_Structure=''FALSE'' \n');
            fprintf('    Please set Dislocation_Core_Structure=''TRUE'' OR prepare the file pnadis.mat \n');
            fprintf('\n');

            return

        end
      
    else  
        
        fprintf('    >> The data file "pnadis.mat" is not needed \n');
        
    end
    
end
    
%% ------------------------------------------------------------------------
% To calculate the Peierls stress
% -------------------------------------------------------------------------
    
if strcmp(Peierls_Stress,'TRUE')

    fprintf('\n');
    fprintf('%% %% INITIATION OF CALCULATING THE PEIERLS STRESS \n');
    fprintf('\n');
        
    if misdim==1
        
        if pnstrmethod==1
            
            fprintf('    >> The Peierls stress was calculated via "Analytical formula" \n');
            run pnstress1d_analysis.m
            
        elseif pnstrmethod==2
            
            fprintf('    >> The Peierls stress was calculated via "Discrete dislocation energy approach" \n');
            run pnstress1d.m
            
        end
    
    elseif misdim==2
        
            if pnstrmethod==1

                fprintf('    >> The Peierls stress was calculated via "Analytical formula" \n');
                run pnstress_analysis.m
                
            elseif pnstrmethod==2

                fprintf('    >> The Peierls stress was calculated via "Discrete dislocation energy approach" \n');
                run pnstress.m
                
            % elseif pnstrmethod==3

                % fprintf('    >> The Peierls stress was calculated via the method "xxx" \n');
                % run pnstress_str.m
                
            end
        
    end
           
    fprintf('\n');
    fprintf('%% %% THE END \n');
    
elseif strcmp(Peierls_Stress,'FALSE')
        
    fprintf('\n');
    fprintf('%% %% NO Peierls stress was calculated \n');
    fprintf('\n');
        
end
    
%% -----------------------------------------------------------------------
% To calculate the pressure field
% ------------------------------------------------------------------------
    
if strcmp(Pressure_Field,'TRUE')

    fprintf('\n');
    fprintf('%% %% INITIATION OF CALCULATING THE PRESSURE FIELD \n');
    fprintf('\n');

    fprintf('    >> The pressure field around dislocation core was calculated \n');
    
    if misdim==1
        [pressfld_latxz_mis, pressfld_latyy_mis, pressfld] = pressfield1d (pressfld_latxz, pressfld_latyy);
    elseif misdim==2
        [pressfld_latxz_mis, pressfld_latyy_mis, pressfld] = pressfield (pressfld_latxz, pressfld_latyy); 
    end
    
    fprintf('\n');
    fprintf('%% %% THE END \n');

elseif strcmp(Pressure_Field,'FALSE')

    fprintf('\n');
    fprintf('%% %% NO pressure field was calculated \n'); 
    fprintf('\n');

end
    
%% ------------------------------------------------------------------------
% To calculate the solid solute strengthening only for 2D P-N model
% ------------------------------------------------------------------------

if misdim==2

    if strcmp(Solution_Strengthening,'TRUE')

        fprintf('\n');
        fprintf('%% %% INITIATION OF CALCULATING THE SOLID SOLUTION STRENGTHENING \n');
        fprintf('\n');

        fprintf('    >> The solid solution strengthening was calculated \n');

        run Solstrengthen.m

        fprintf('\n');
        fprintf('%% %% THE END \n');

    elseif strcmp(Solution_Strengthening,'FALSE')

        fprintf('\n');
        fprintf('%% %% NO solid solution strengthening was calculated \n'); 
        fprintf('\n');

    end
  
end

%% ------------------------------------------------------------------------
% To save the workspace variables to file
% ------------------------------------------------------------------------

clear i j ii jj ij II JJ;

% % Used to define the variable "gammafitresult"
clear temp_AA temp_BB temp_gfit;
clear gammafunc;
% % the angle between the dislocation line and its Burgers vector
clear theta;

% % The temporary variable of fitting gamma surface, fitting restoring stress
%   and minimizing the dislocation energy
clear nvars;
clear RestoringForcefitresult;
clear toolbox_result lsq_result pso_result ga_result;
clear variable_ga variable_pso fval;
clear gof opts options rsquare resnorm;

% % The plot with two y-axes
clear hAx hLine1 hLine2;
% % In the file of Default.m
clear defaultend;
% % Used to save the figure file xxx.fig
clear figurefile;
% % In the files of pnstress1d.m and pnstress.m
clear u Epn_temp space;
clear dsfe dsfe_fliplr;
clear dsfe_x dsfe_x_fliplr;
clear str_is_positive str_is_negative;

% % The temporary variable of solid solution strengthening
clear Eint_plot_temp;

% For the gradient, nonlocal and stress energy, which are not included in V_1.0.0
clear GradientEnergy NonlocalEnergy StressEnergy;

save('pnadis.mat');

fprintf('\n');
fprintf('%% %% The workspace variables was saved in file "pnadis.mat" \n'); 
fprintf('\n');
