% Copyright: Wing Kam Liu, Adrian Kopacz, Eduard Karpov,Harold Park,
%            David Farrell, Mechanical Engineering,
%            Northwestern University, 2004 


function [hx, hy, vhx, vhy] = md2d_sol(x,y,vx,vy,sL,tS,nIt,bC) 

% input
% initial positions (x,y)
% initial velocities (vx,vy)
% space dimension (sL)
% time step (tS)
% number of iterations (nIt)
% boundry condition (1 = fixed) (bC)

% output
% position of particles for each iteration (hx, hy)
% velocity of particles for each iteration (hvx, hvy)
   
% initialize positions and velocities to zero
[nP,n]=size(x);

hx=zeros(nP,nIt);
hy=zeros(nP,nIt);
vhx=zeros(nP,nIt);
vhy=zeros(nP,nIt);

%hx(:,1)=x-vx*tS;
%hy(:,1)=y-vy*tS;

% previous position
xp=x-vx*tS;
yp=y-vy*tS;

% solving Newton's equation (Verlet Algorithm)  
% Enter with current (x,y) and previous (xp, yp) positions 
for it=1:nIt
    
    % initialize forces to zero
    fx=zeros(nP,nP);
    fy=zeros(nP,nP);  
    
    for i=1:nP-1               
        for j=i+1:nP
            dx=x(i)-x(j); 
            dy=y(i)-y(j);
            
            % check for boundary conditions 
            if(bC == 0)
                if(dx>sL/2) 
                    dx=dx-sL; 
                end 
                if(dx<-sL/2) 
                    dx=dx+sL; 
                end 
                if(dy>sL/2) 
                    dy=dy-sL; 
                end 
                if(dy<-sL/2) 
                    dy=dy+sL; 
                end 
            end
            r=sqrt(dx^2+dy^2); 
            
            % Lennard-Jones model
            f=24*(2/r^13-1/r^7); 
            fx(i,j)=-f*dx/r;
            fy(i,j)=-f*dy/r; 
        end 
    end 
    fx=fx-fx'; 
    fy=fy-fy'; 
    
    % update the position   
    fxSum=(sum(fx))';
    fySum=(sum(fy))';
    
    % save positions
    hx(:,it)=x;
    hy(:,it)=y;
    
    x=2*x-xp+fxSum*tS^2; 
    y=2*y-yp+fySum*tS^2; 
    
    % save velocities
    vhx(:,it)=(x-xp)/(2*tS); 
    vhy(:,it)=(y-yp)/(2*tS); 
    
    % check for boundary condition 
    if(bC == 0)
        obC=find(x<0); 
        x(obC)=x(obC)+sL; 
        hx(obC,it)=hx(obC,it)+sL; 
        
        obC=find(x>sL); 
        x(obC)=x(obC)-sL; 
        hx(obC,it)=hx(obC,it)-sL; 
        
        obC=find(y<0); 
        y(obC)=y(obC)+sL; 
        hy(obC,it)=hy(obC,it)+sL; 
        
        obC=find(y>sL); 
        y(obC)=y(obC)-sL; 
        hy(obC,it)=hy(obC,it)-sL; 
    end;
    
    % fixed boundary condition
    if(bC == 1)
        obC=find(x<0); 
        vhx(obC,it)=-1*vhx(obC,it);
        x(obC)=x(obC)*-1; 
        hx(obC,it)=hx(obC,it)*-1; 
        
        obC=find(x>sL); 
        vhx(obC,it)=-1*vhx(obC,it);
        x(obC)=((x(obC)-sL)*-1)+sL; 
        hx(obC,it)=((hx(obC,it)-sL)*-1)+sL; 
        
        obC=find(y<0); 
        vhy(obC,it)=-1*vhy(obC,it);
        y(obC)=y(obC)*-1; 
        hy(obC,it)=hy(obC,it)*-1; 
        
        obC=find(y>sL); 
        vhy(obC,it)=-1*vhy(obC,it);
        y(obC)=((y(obC)-sL)*-1)+sL; 
        hy(obC,it)=((hy(obC,it)-sL)*-1)+sL; 
    end;
    
    xp = hx(:,it); 
    yp = hy(:,it); 
end;
