% 2-D Molecular Dynamics: Periodic Boundary Conditions

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

% Main code, calls all others
close all;
clear all;

sL= 21;              % space dimension of primary image
npu= 14;            % number of particles per each axis
nP= npu^2; % number of particles
mV=50;             % magnitude of initial velocities
tS=0.00000005;           % time step 0.0000001
nIt=2000;           % number of iterations
bC=0;              % boundry condition (0 = Periodic, 1 = Fixed, 2 = None)
mI=0;              % number of translation images (0 = primary image)
tR=0;              % trace option (0 = no trace) (1 = trace)
fN='MD-CLuster_PeriodicBoundary_dense_random';     % name of video file
fps = 40;          % movie output (frames per second)

% initialize positions and velocities to zero
x=zeros(nP,1);
y=zeros(nP,1);
vx=zeros(nP,1);
vy=zeros(nP,1);

% assign positions and velocities
for j=1:npu
    for i=1:npu
        index=(j-1)*npu+i;
        x(index)=((i-0.5)/npu)*sL;
        y(index)=((j-0.5)/npu)*sL;
    end
end

% Random initial velocities
for j=1:npu
    for i=1:npu
        index=(j-1)*npu+i;
        vx(index)=mV*10000*rand(1); % mV instead of sL
        vy(index)=mV*10000*rand(1);
    end
end

% solving Newton's equation
[hx, hy, vhx, vhy] = md2d_sol(x,y,vx,vy,sL,tS,nIt,bC);

% simulating results
md2d_sim(hx,hy,sL,nIt,mI,tR,fN,fps);
