% % Gas Particles - No Forces
% % This code just has a bunch of particles bounce around in
% % a 2d area for a while in order to see if their velocities
% % converge to something like a Maxwell Distribution.

% % If you want to just watch the particles bounce around, 
% % change the if statements below.

% % Henry Peterson 2016

clearvars;
close all;

%number of particles
nParts = 100;
%length of the axes (walls)
L = 1;
%small change in time
dt = 1e-2;
%buffer that will determine if particles collide
buff = 1e-2;

%setup initial position arrays
x = rand(1, nParts)*L;  %btwn 0 and L
y = rand(1, nParts)*L;  %btwn 0 and L

%setup initial velocity arrays
maxv = .25;
minv = -maxv;
vx = rand(1, nParts)*(maxv-minv) + minv;  %btwn min and max
vy = rand(1, nParts)*(maxv-minv) + minv;  %btwn min and max


%this part is just for the distribution check at the end
max = maxv * 2;            %endpoint on domain
bars = 20;                 %num of bars
dv = max/(bars-1);         %change between bar limits
dist = zeros(max/dv+1,1);  %number of bars in an interval
limits = 0:dv:max;         %the limits of intervals


%keep track of the frame number (if you want to save images and make a gif
frame = 1;

%initial distribution
%for each particle, put its velocity btwn the right limits
v = sqrt(vx.^2+vy.^2);
for j=1:nParts
    for k=1:length(limits)
        if k == length(limits)
            break;
        end
        if v(j) > limits(k) && v(j) < limits(k+1)
            dist(k) = dist(k) + 1;
        end
    end
end
figure()
plot(limits,dist,'-ro');
title('Initial Velocity Distribution of Simulated Gas Particles')
xlabel('Velocity')
ylabel('Frequency')
set(gca,'XTick',[])
set(gca,'YTick',[])
pause(dt)


itMax = 1e6; 
for it=1:itMax
    
    %change this to a different ratio if you want the code
    %to run without showing a plot
    if it/itMax > 0
        %scatter plot for different colors (to see collisions)
        p = scatter(x, y, 25, 1:nParts, 'filled');
        axis([0,L,0,L]);
        colormap hsv
        set(gca,'YTick',[],'XTick',[],'Color','black')
        set(gcf,'Color','white','InvertHardCopy','off')  %InvertHardCopy makes the gif also black
        title('Gas Particles - No Forces')

        pause(dt);
        
        %to save images
        print(gcf,['gas_',num2str(frame,'%05d'),'.png'],'-dpng','-r300')
        frame = frame + 1;
        if frame > 250
            break
        end
        
        %check if the user closed the plot 'p' and break if so
        if ishandle(p) == 0
            break;
        end
    end
    
    %move parts
    x = x + vx*dt;
    y = y + vy*dt;
    
    %catch collisions
    used = zeros(1,nParts);
    I = 1;
    for i=1:nParts-1
        for j=(i+1):nParts
            if find(used == j)
                continue;
            end
            if abs(x(i) - x(j))<buff && abs(y(i) - y(j))<buff
                used(I) = j;
                I = I + 1;
                
                %here we want to switch components along collision:                
                %the vector along collision
                r = [x(i)-x(j),y(i)-y(j)];  
                
                %compositions of each particle's velocity vector onto r
                compi = (vx(i)*r(1) + vy(i)*r(2))/norm(r);  
                compj = (vx(j)*r(1) + vy(j)*r(2))/norm(r);
                
                %turn compositions into vectors along r
                proji = compi*[r(1)/norm(r),r(2)/norm(r)];  
                projj = compj*[r(1)/norm(r),r(2)/norm(r)];
                
                %replace the projections (assuming same mass)
                vx(i) = vx(i) - proji(1) + projj(1);        
                vx(j) = vx(j) - projj(1) + proji(1);
                vy(i) = vy(i) - proji(2) + projj(2);
                vy(j) = vy(j) - projj(2) + proji(2);
                
                break;
            end
        end
    end
    
    
    %catch walls
    for i=1:nParts
        if x(i) < 0
            x(i) = 0;
            vx(i) = -vx(i);
        elseif x(i) > L
            x(i) = L;
            vx(i) = -vx(i);
        end
        
        if y(i) < 0
            y(i) = 0;
            vy(i) = -vy(i);
        elseif y(i) > L
            y(i) = L;
            vy(i) = -vy(i);
        end
    end
    
    % progress for user
    if mod(it,itMax/100)==0
        fprintf('percent done: %.0f %%\n',it/itMax*100);
    end
    
    % now we gather data for the distribution of velocities
    % ignore first 10% of time (let things settle)
    if it>itMax/10
        v = sqrt(vx.^2+vy.^2);
        %for each particle, put its velocity btwn the right limits
        for j=1:nParts
            for k=1:length(limits)
                if k == length(limits)
                    break;
                end
                if v(j) > limits(k) && v(j) < limits(k+1)
                    dist(k) = dist(k) + 1;
                end
            end
        end
    end
end

% after long time of running the code (preferably without the 
% visualization) it should look pretty nice! Plot the 
% velocity dist here.
figure()
plot(limits,dist,'-ro');
title('Velocity Distribution of Simulated Gas Particles')
xlabel('Velocity')
ylabel('Frequency')
set(gca,'XTick',[])
set(gca,'YTick',[])

%thanks for checking out my code! 