clf;
clear;
clc;

global nElectrons q m k nProtons n X Y Z
n=5;  %number of rows of electrons and protons for face-centered lattice
nElectrons=(n^3)+(n-1)^3;
nProtons=1;

X=15;  %bounding for proton locations
Y=15;
Z=12;

k=8.99e9;  %epsilon naught constant
q=-1.602176565e-19; %charge on an electron;
m=9.10938e-31;

dt=.0015;
nn=1000;
t=0;
f=0;

for i=1:nProtons   %set proton positions (randomly, but clustered to some extent in center)
    xp(i)=i*X/nProtons;
    yp(i)=Y/3+(rand()*Y/3);
    zp(i)=Z/3+(rand()*Z/3);
end

yProton=[xp yp zp];

for i=0:nProtons-1  %put into format of x(i) y(i) z(i)
    rp(3*i+1:3*i+3)=yProton(i+1:nProtons:nProtons*3);
end
yProton=rp; %set location of protons.

yElectron=initial_conditions_electrons(n);

s=zeros(3*nElectrons,1);

yElectron=[yElectron; s];  %put electron positions and velocities in vector of [... x(i) y(i) z(i) ... vx(i) vy(i) vz(i)...]
yElectron=initial_velocities(yElectron);

for it=1:nn
    f1=ode(yElectron,yProton);   % Runge Kutta method for new electron positions, velocities
    f2=ode(yElectron+(dt/2).*f1,yProton);
    f3=ode(yElectron+(dt/2).*f2,yProton);
    f4=ode(yElectron+dt.*f3,yProton);
    
    yElectron=yElectron+dt/6.*(f1+f2+f3+f4);
        
    yElectron=reflective_boundary_conditions(yElectron);  %check to see whether any electrons are outside of boundary of box
    
    xxe = [yElectron(1:3:3*nElectrons)];
    yye = [yElectron(2:3:3*nElectrons)];
    zze = [yElectron(3:3:3*nElectrons)];
    
    xxp=[yProton(1:3:3*nProtons)];
    yyp=[yProton(2:3:3*nProtons)];
    zzp=[yProton(3:3:3*nProtons)];
    
    
    
    if mod(it+1,2) == 0
       f=f+1;
      scatter3(xxe,yye,zze,'MarkerEdgeColor','m');   %use to plot current location of all electrons. protons not included. previous locations of electrons not included
     set(gca,'plotboxaspectratio',[2 2 2])
    axis([0 15 0 15 0 15])
    title('Electron motion')
    saveas(gcf, ['emotion',sprintf('%03d',f),'.png'])
    pause(0.1)
    %t=t+dt;
    end
    %fprintf('%f\n', t)
    
    %if mod(it+1,2) == 0
       % f=f+1;
       % scatter3(xxe,yye,zze,'m')  %plot both protons, electrons with trajectory of movement
       % hold on
       % scatter3(xxp,yyp,zzp,'MarkerEdgeColor','k','MarkerFaceColor','k');
       % hold off
       % set(gca,'plotboxaspectratio',[2 2 2]);
        %axis([0 X 0 Y 0 Z])
       % title('Mapping electron motion with protons');
      %  legend('electrons','protons')
     %   saveas(gcf, ['emotiontrajectory',sprintf('%03d',f),'.png'])
    %    pause(0.1)
    %end
end





