%Final Project 
%Rustin Domingos 
%EPS 109 - Fall 2015

clear

%Data for copper
global bulk m r L kb N eps sigma
m = 1.05520675*(10^-25); %mass of copper atom
r = 354.5706371; %scaled from 1.28A
kb = 105.6551221; %scaled
eps = 714804.5122; %scaled
sigma = 627.9565133; %scaled
bulk = .0506086; %scaled
L = 2500; %scale everything so 1000 = 3.61A

N = 96; %number of particles
f = 0;


%initialize vectors
y = initialconditions();

rv = zeros(3,N); %position vector
count = 1;
for i = 1:6:N*6-5
    rv(:,count) = [y(i); y(i+1); y(i+2)];
    count = count + 1;
end

vv = zeros(3,N); %velocity vector
count = 1;
for i = 4:6:N*6-2
    vv(:,count) = [y(i); y(i+1); y(i+2)];
    count = count + 1;
end

Y = zeros(6,N); %main vector
for i = 1:N
    Y(1:3,i) = rv(:,i);
    Y(4:6,i) = vv(:,i);
end

%Scale Velocities and Calculate Temp
Y(4:6,:) = Y(4:6,:);
K = kinetic_energy(Y);
%T = K*2/3/N/kb;


%Perform Simulation
tMax = 3*10^-12;
t = 0;
dt = 10^-15;
YN = Y;
greencount = .45;

while t < tMax
    for i = 1:10
        for particle = 1:N
            F1 = ode(t,Y(:,particle),Y,particle);
            F2 = ode(t+dt/2,Y(:,particle) + dt/2.*F1,Y,particle);
            F3 = ode(t+dt/2,Y(:,particle) + dt/2.*F2,Y,particle);
            F4 = ode(t+dt,Y(:,particle) + dt.*F3,Y,particle);
            YN(:,particle) = Y(:,particle) + dt/6.*(F1+2.*F2+2.*F3+F4);
        end
         Y = YN;
        t = t + dt;
    end
if t <= 2*10^-12
    Y(4:6,:) = Y(4:6,:)*1.01158;
    greencount = greencount -.001;
end


    scatter3(Y(1,:),Y(2,:),Y(3,:),6000,'MarkerEdgeColor','k','MarkerFaceColor',[.75 greencount 0]);
    axis([-500,3000,-500,3000,-500,3000])
    view([-65,23])
    set(gca,'Zdir','reverse')
    %title(['T = ',num2str(T),' K'],'Color','r')
    axis off
    saveas(gcf, ['copper_',sprintf('%03d',f),'.png'])
    f = f + 1;
    clf
end
    




