
% ------------------------------------------------------------------------
% Function: output_Edis
% ------------------------------------------------------------------------
% Aim:
% 1/ To determine the elasitc and misfit energies as the total
%    dislocation energy is minimum
% ------------------------------------------------------------------------

function output_Edis

    global pnmode;
    global Ax Dx Wx dx Az Dz Wz dz;
    global BurVect;
    global Kx Kz;
    global Nmis;
    
    global shear_modulus;
    global dax Inpas;
    global StressEnergy GradientEnergy NonlocalEnergy;
    global str_coef grad_coef nonloc_coef; 
    
    if pnmode==1

        % % outer cut-off radius

        Rcut=5000*BurVect;
        
        % % elastic energy (1e-11*J/m)

        if Nmis>=1

            % % To define some Matrix

            Xcoef=zeros(1,Nmis-1);
            Xdist=zeros(1,Nmis-1);
            Xwid=zeros(1,Nmis);

            Zcoef=zeros(1,Nmis-1);
            Zdist=zeros(1,Nmis-1);
            Zwid=zeros(1,Nmis);

            % % x direction

            for j=1:(Nmis-1)
                Xcoef(j)=Ax(j);
            end

            if Nmis==1
                Xcoef(Nmis)=1;
            else 
                Xcoef(Nmis)=1-sum(Xcoef(1:(Nmis-1)));
            end

            for j=1:(Nmis-1)
                Xdist(j)=Dx(j);
            end

            if Nmis==1
                Xdist(Nmis)=0;
            else 
                Xdist(Nmis)=0-sum(Xdist(1:(Nmis-1)));
            end

            for j=1:Nmis
                Xwid(j)=Wx(j);
            end

            % % z direction

            for j=1:(Nmis-1)
                Zcoef(j)=Az(j);
            end

            if Nmis==1
                Zcoef(Nmis)=0;
            else 
                Zcoef(Nmis)=0-sum(Zcoef(1:(Nmis-1)));
            end

            for j=1:Nmis-1
                Zdist(j)=Dz(j);
            end

            if Nmis==1
                Zdist(Nmis)=0;
            else 
                Zdist(Nmis)=0-sum(Zdist(1:(Nmis-1)));
            end

            for j=1:Nmis
                Zwid(j)=Wz(j);
            end

            % % x energy

            [Xcoefm,Xcoefn]=meshgrid(Xcoef,Xcoef);
            [Xdistm,Xdistn]=meshgrid(Xdist,Xdist);
            [Xwidm,Xwidn]=meshgrid(Xwid,Xwid);

            Xtemp=Xcoefm.*Xcoefn.*log(1+(Xdistm-Xdistn).^2./(Xwidm+Xwidn).^2);
            Xtempr=Xcoefm.*Xcoefn.*log(Rcut./abs(Xwidm+Xwidn));

            % % z energy

            [Zcoefm,Zcoefn]=meshgrid(Zcoef,Zcoef);
            [Zdistm,Zdistn]=meshgrid(Zdist,Zdist);
            [Zwidm,Zwidn]=meshgrid(Zwid,Zwid);

            Ztemp=Zcoefm.*Zcoefn.*log(1+(Zdistm-Zdistn).^2./(Zwidm+Zwidn).^2);
            Ztempr=Zcoefm.*Zcoefn.*log(Rcut./abs(Zwidm+Zwidn));

            Eel=1/(4*pi)*Kx*BurVect^2*(sum(sum(Xtempr))-1/2*sum(sum(Xtemp)))...
               +1/(4*pi)*Kz*BurVect^2*(sum(sum(Ztempr))-1/2*sum(sum(Ztemp)));      

        elseif Nmis==0

            Eel=1/(4*pi)*Kx*(BurVect^2*log(abs(1/(2*Wx)))-Ax*(1-Ax)*BurVect^2*log(1+Dx^2/(2*Wx)^2))...
               -1/(4*pi)*Kz*(Az*(-Az)*BurVect^2*log(1+Dz^2/(2*Wz)^2));

        elseif Nmis==-1

            Xcoef=[1/2*(1/2-Ax);Ax;1/2*(1/2-Ax);1/2*(1/2-Ax);Ax;1/2*(1/2-Ax)];
            Xdist=[Dx/2+dx;Dx/2;Dx/2-dx;-(Dx/2+dx);-Dx/2;-(Dx/2-dx)];
            Zcoef=[1/2*(sqrt(3)/6-Az);Az;1/2*(sqrt(3)/6-Az);-1/2*(sqrt(3)/6-Az);-Az;-1/2*(sqrt(3)/6-Az)];
            Zdist=[Dz/2+dz;Dz/2;Dz/2-dz;-(Dz/2+dz);-Dz/2;-(Dz/2-dz)];

            [Xcoefm,Xcoefn]=meshgrid(Xcoef,Xcoef);
            [Xdistm,Xdistn]=meshgrid(Xdist,Xdist);
            [Zcoefm,Zcoefn]=meshgrid(Zcoef,Zcoef);
            [Zdistm,Zdistn]=meshgrid(Zdist,Zdist);

            Xtemp=Xcoefm.*Xcoefn.*log(1+(Xdistm-Xdistn).^2/(2*Wx)^2);
            Ztemp=Zcoefm.*Zcoefn.*log(1+(Zdistm-Zdistn).^2/(2*Wz)^2);        

            Eel=1/(4*pi)*Kx*BurVect^2*(log(abs(Rcut/(2*Wx)))-1/2*sum(sum(Xtemp)))-1/(4*pi)*Kz*BurVect^2*1/2*sum(sum(Ztemp));

        end

        % %   misfit energy (1e-10*J/m)

        EA=quadgk(@Gamma,-Rcut,Rcut,'AbsTol',1E-15,'RelTol',1E-10,'MaxIntervalCount',2000);

        % %   total energy(unit: 1e-10*J/m) of a dislocation

        MIN=1E-1*Eel + EA;
        
        fprintf('\n');
        fprintf('    >> The results of the various energies of dislocation (1e-10*J/m)\n');
        fprintf('\n');
        fprintf('       Total energy (Etot) : %10.4f \n', MIN);
        fprintf('       Elastic energy (Eel): %10.4f \n', 1E-1*Eel);
        fprintf('       Misfit energy (EA)  : %10.4f \n', EA);  
        
    elseif pnmode==2

        % % outer cut-off radius

        Rcut=50;           
        
        % % elastic energy (1e-11*J/m)

        [n,nn]=meshgrid(((-Rcut+dax):(Rcut+dax)),((-Rcut+dax):(Rcut+dax)));

        Chi=3/2*Inpas^2+1/2*((n-1)*Inpas-(nn-1)*Inpas).^2.*log(abs((n-1)*Inpas-(nn-1)*Inpas))...
                      +1/2*(n*Inpas-nn*Inpas).^2.*log(abs(n*Inpas-nn*Inpas))...
                      -1/2*(n*Inpas-(nn-1)*Inpas).^2.*log(abs(n*Inpas-(nn-1)*Inpas))...
                      -1/2*((n-1)*Inpas-nn*Inpas).^2.*log(abs((n-1)*Inpas-nn*Inpas));

        Chi(1,1)=3/2*Inpas^2-Inpas^2*log(Inpas);
        for m=2:(2*Rcut+1)
            Chi(m,m)=3/2*Inpas^2-Inpas^2*log(Inpas);
            Chi(m,m-1)=3/2*Inpas^2+Inpas^2*log(Inpas)-1/2*(2*Inpas)^2*log(2*Inpas);
            Chi(m-1,m)=Chi(m,m-1);
        end

        % Eel=1/(4*pi)*Kx*BurVect^2*log(Rcut*BurVect)+sum(sum(Chi.*(1/(4*pi)*Kx*(misfit_x(n*Inpas)-misfit_x((n-1)*Inpas))/Inpas.*(misfit_x(nn*Inpas)-misfit_x((nn-1)*Inpas))/Inpas...
        %                                                          +1/(4*pi)*Kz*(misfit_z(n*Inpas)-misfit_z((n-1)*Inpas))/Inpas.*(misfit_z(nn*Inpas)-misfit_z((nn-1)*Inpas))/Inpas)));

        Eel=sum(sum(Chi.*(1/(4*pi)*Kx*(misfit_x(n*Inpas)-misfit_x((n-1)*Inpas))/Inpas.*(misfit_x(nn*Inpas)-misfit_x((nn-1)*Inpas))/Inpas...
                         +1/(4*pi)*Kz*(misfit_z(n*Inpas)-misfit_z((n-1)*Inpas))/Inpas.*(misfit_z(nn*Inpas)-misfit_z((nn-1)*Inpas))/Inpas)));        
        
        % % misfit energy (1e-10*J/m)

        mm=-(Rcut*10+dax):(Rcut*10+dax);
        EA=sum((Gamma(mm*Inpas)*Inpas));

        % % stress interaction term (1e-10*J/m and unit of str_coef: MPa)
    
        if StressEnergy==0
    
            Estr=0.0000;
            
        elseif StressEnergy==1
            
            % Estr=-1E-4*str_coef*sum(((mm*Inpas).^2-((mm-1)*Inpas).^2)/2.*(misfit_x(mm*Inpas)-misfit_x((mm-1)*Inpas))/Inpas);
            Estr=-1E-4*str_coef*sum((misfit_x(mm*Inpas)+misfit_x((mm+1)*Inpas))*Inpas)/2;
        
        end
        
        % % gradient energy (1e-10*J/m and unit of grad_coef: shear_modulus*Burvect)
    
        if GradientEnergy==0
            
            Egradient=0.0000;
            
        elseif GradientEnergy==1
            
            Egradient=1E-1*grad_coef*shear_modulus*BurVect*Inpas*(sum(((misfit_x((mm+1)*Inpas)-misfit_x((mm-1)*Inpas))/(2*Inpas)).^2)...
                                                                 +sum(((misfit_z((mm+1)*Inpas)-misfit_z((mm-1)*Inpas))/(2*Inpas)).^2));
                                  
        end
                              
        % % simplified nonlocal energy (1e-10*J/m and unit of nonlocal_coef: shear_modulus/Burvect)
        
        if NonlocalEnergy==0
            
            Enonlocal=0.0000;
        
        elseif NonlocalEnergy==1
            
            Enonlocal=1E-1*nonloc_coef*shear_modulus/BurVect*(sum(misfit_x(mm*Inpas).*(misfit_x(mm*Inpas)-(misfit_x((mm+1)*Inpas)+misfit_x((mm-1)*Inpas))/2)*Inpas)...
                                                             +sum(misfit_z(mm*Inpas).*(misfit_z(mm*Inpas)-(misfit_z((mm+1)*Inpas)+misfit_z((mm-1)*Inpas))/2)*Inpas));
        
        end
                          
        % % total energy(unit: 1e-10*J/m) of a dislocation

        MIN=1E-1*Eel + EA + Estr + Egradient + Enonlocal;

        fprintf('\n');
        fprintf('    >> The results of the various energies of dislocation (1e-10*J/m)\n');
        fprintf('\n');
        fprintf('       Total energy (Etot)         : %10.4f \n', MIN);
        fprintf('       Elastic energy (Eel)        : %10.4f \n', 1E-1*Eel);
        fprintf('       Misfit energy (EA)          : %10.4f \n', EA); 
        % fprintf('       Gradient energy (Egradient) : %10.4f \n', Egradient); 
        % fprintf('       Nonlocal energy (Enonlocal) : %10.4f \n', Enonlocal); 
        % fprintf('       Stress Energy (Estr)        : %10.4f \n', Estr); 
    
    end
    
end

