
% ------------------------------------------------------------------------
% Function: minfunc ( Input )
% ------------------------------------------------------------------------
% Aim:
% The total energy of a dislocation according the CCPN model
% ------------------------------------------------------------------------
% Ref:
% 1/ Schoeck, Gunther. Materials Science and Engineering: A 400 (2005): 7-17.
% ------------------------------------------------------------------------

function MIN = minfunc ( Input )

    %   defining global variables

    global Dx Wx Dz Wz Ax Az;
    global BurVect;
    global Kx Kz;
    global dx dz;
    global Nmis;
    
    %   outer cut-off radius

    Rcut=5000*BurVect;
    
    %%   >>>>>> enter into main loop <<<<<<

    nx=size(Input,1);
    
    %   To define some Matrix
    
    MIN=zeros(nx,1);
    
    for i=1:nx
        
        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)
                Ax(j)=Input(i,j);
                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)
                Dx(j)=Input(i,Nmis-1+j);
                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
                Wx(j)=Input(i,2*Nmis-2+j);
                Xwid(j)=Wx(j);
            end
            
            % % z direction
    
            for j=1:(Nmis-1)
                Az(j)=Input(i,3*Nmis-2+j);
                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
                Dz(j)=Input(i,4*Nmis-3+j);
                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
                Wz(j)=Input(i,5*Nmis-4+j);
                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));
       
            % % elastic energy (1e-11*J/m)
            
            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

            %   some parameters

            Ax=1/2;
            Az=sqrt(3)/6;

            Dx=Input(i,1);
            Wx=Input(i,2);
            Dz=Input(i,3);
            Wz=Input(i,4);

            % % elastic energy (1e-11*J/m)
            
            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

            Ax=Input(i,1);
            Dx=Input(i,2);
            Wx=Input(i,3);
            dx=Input(i,4);
            
            Az=Input(i,5); 
            Dz=Input(i,6);
            Wz=Input(i,7);
            dz=Input(i,8);

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

            % % elastic energy (1e-11*J/m)            
            
            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);
        % EA=integral(@Gamma,-Rcut,Rcut,'AbsTol',1E-15,'RelTol',1E-10);
        
        %%   total energy(unit: 1e-10*J/m) of a dislocation

        MIN(i,1)=1E-1*Eel + EA;
        
    end
    
end
