% % Charged Wire Simulation
% % This code uses my collision function, a potential function
% % modified for Coulomb forces, and a ODE solver to simulate
% % many charged particles moving in a wire with an electric
% % field and resistance.

% % Henry Peterson 2016

clearvars;
close all;

% the E field, Length of wire, 
% and c (a constant for potential) will be needed 
% in functions (collisions2.m, ode2.m)
global E L c
% electric field (a constant vector field - start at 0 for my video)
E = [0, 0];   
% length of axes
L = 1;
% this constant plays a big role in the potential.m function
c = 0.1;

% number of particles
nParts = 100;
% number of particles which are protons (nuclei with positive net charge)
nProts = 40;
% naming some spots on the wire makes code easier to read:
left = L/8;
length = L/2;
right = left + length;
% small change in time
dt = 2e-4;

% % set up arrays
% charge of each individual particle stored in 'qs'
q = 10;
qs = -q*ones(1,nParts);   %electrons negative
qs(1:nProts) = q;         %protons positive
% mass of each particle in 'ms'
mE = 1;                   %mass of electron
mP = 1e4;                 %mass of proton
ms = mE*ones(1,nParts);
ms(1:nProts) = mP;
%%%%NOTE - since (on the atomic level) gravitational forces are much
%%%%less important than electric forces, I do not account for them in
%%%%my potential function. instead, mass is used to calculate acceleration
%%%%and conserve momentum in collisions
%initial x and y positions of each particle
x = zeros(1,nParts);
y = zeros(1,nParts);
%setup the prots in a lattice
I = 1;  %index to keep track
Nx = 5; %number of prots in x dir
Ny = nProts/Nx;          %in y dir
for i=1:Nx
    for j=1:Ny
        x(I) = 13*length/15 * (i-1)/(Nx-1) + left + length/15;
        y(I) = 0.9*L * (j-1)/(Ny-1) + mod(i,2)*0.05 + L/50;
        I = I + 1;
    end
end
%electrons in random positions in wire
for i=I:nParts
    x(i) = rand()*length + left;
    y(i) = rand()*L;
end
%initial x and y velocities of each particle set to zero
vx = zeros(1,nParts);
vy = zeros(1,nParts);

%to count frames for movie (if desired)
frame = 1;

while 1
    
    %collisions (just electrons)
    [vx,vy] = collisions2(x,y,vx,vy,ms,nProts);
    
    %ode2
    A = [x; y; vx; vy];
    F1 = ode2(A, ms, qs);
    F2 = ode2(A + F1*dt/2, ms, qs);
    F3 = ode2(A + F2*dt/2, ms, qs);
    F4 = ode2(A + F3*dt, ms, qs);
    A = A + dt/6*(F1 + 2*F2 + 2*F3 + F4);
    
    x = A(1,:);
    y = A(2,:);
    vx = A(3,:);
    vy = A(4,:);
    
    %catch walls
    for i=1:nParts
        if x(i) < left
            x(i) = left;
            vx(i) = -vx(i);
        elseif x(i) > right
            x(i) = right;
            vx(i) = -vx(i);
        end
        %periodic boundary conditions on the top and bottom!
        if y(i) < 0
            y(i) = L + y(i);   
        elseif y(i) > L
            y(i) = y(i)-L;   
        end
    end
    
    %calculate the avg drift velocity of the electrons 
    v_drift = sum(vy)/(nParts-nProts);
    
    %plot
    p = plot(x((nProts+1:end)),y((nProts+1):end),'.',...
        'MarkerSize',8,'MarkerFaceColor','y','MarkerEdgeColor','y');
    
    hold on
    plot(x(1:nProts),y(1:nProts),'.',...
        'MarkerSize',20,'MarkerFaceColor','r','MarkerEdgeColor','r');
    
    plot([left,left],[0,L],'b-',...
        [right,right],[0,L],'b-');
    hold off
    
    text(right + 0.05*L, 0.85*L,...
        strcat( 'E = ', num2str( -E(2) ) ), 'Color', 'yellow', 'FontSize', 14)
    text(right + 0.05*L, 0.7*L,...
        strcat( 'v_{drift} = ', num2str(v_drift,2) ), 'Color', 'white', 'FontSize', 14)
    
    set(gca,'YTick',[],'XTick',[],'Color','black')
    set(gcf,'Color','white','InvertHardCopy','off')  %InvertHardCopy makes the gif also black
    title('Charged Particles In A Conducting Wire')
    axis([0, L, 0, L]);
    legend({'Electrons (-)','Nuclei (+)'},'Location','southeast','Color','white')
    
    pause(dt);
    
    %if you want to make a movie, use this print statement to save images
    print(gcf,['charges_',num2str(frame,'%05d'),'.png'],'-dpng','-r300');
    frame = frame + 1;
    %the following are some automatic adjustments that I used for a nice
    %movie in my presentation. fiddle with them to see if you can find
    %better results!
    if frame == 50
        E = [0, -500];
    elseif frame==200
        E = [0, 1000];
    elseif frame>300
        break;
    end
    %give the user an idea of status
    if mod(frame,10)==0
        fprintf('frame: %d\n',frame);
    end
    
    %break if user closed the plot
    if ishandle(p) == 0
        break;
    end
end

%thanks for checking out my code!