xMin = 0; xMax = 100; xWall = 50;
yMin = 0; yMax = 100; holeMin = 45; holeMax = 55;
particles = 5;
minVelocity = -3;
maxVelocity = 3;


% Create wall.
wallLength = yMax - (holeMax - holeMin);
wallX = zeros(wallLength, 1);
for i = 1:size(wallX)
    wallX(i) = xWall;
end
wallY = zeros(wallLength,1);
for i = 1:holeMin
    wallY(i) = i;
end
for i = holeMin + 1:wallLength
    wallY(i) = i + (yMax - wallLength);
end

% Create particle starting positions.
positionsX = zeros(particles,1); positionsY = zeros(particles,1);
for i = 1:particles
    positionsX(i) = xWall*rand();
    positionsY(i) = yMax*rand();
end

% Create particle initial velocities.
velocitiesX = zeros(particles,1); velocitiesY = zeros(particles,1);
for i = 1:particles
    velocitiesX(i) = (maxVelocity - minVelocity) * rand() + minVelocity;
    velocitiesY(i) = (maxVelocity - minVelocity) * rand() + minVelocity;
end

plot(positionsX, positionsY, 'o');

for t = 1:300
    for i = 1:particles
        
        for j = 1:particles
            if i ~= j % if not comparing position of same particle
                if abs(positionsX(i) - positionsX(j)) <= 1 && abs(positionsY(i) - positionsY(j)) <= 1
                    disp "collision"
                    tempX = velocitiesX(i);
                    tempY = velocitiesY(i);
                    velocitiesX(i) = velocitiesX(j);
                    velocitiesY(i) = velocitiesY(j);
                    velocitiesX(j) = tempX;
                    velocitiesY(j) = tempY;
                    break
                end
            end
        end
        
        newXposition = positionsX(i) + velocitiesX(i);
        newYposition = positionsY(i) + velocitiesY(i);
%         positionsX(i) = positionsX(i) + velocitiesX(i);
%         positionsY(i) = positionsY(i) + velocitiesY(i);
        
        if newXposition < 0
            velocitiesX(i) = velocitiesX(i) * -1;
            positionsX(i) = newXposition * -1;
        elseif newXposition > xMax
            velocitiesX(i) = velocitiesX(i) * -1;
            positionsX(i) = xMax - (newXposition - xMax);
        elseif newXposition > xWall && positionsX(i) <= xWall
            if newYposition >= holeMax || newYposition <= holeMin
                %disp 'wall bounce';
                velocitiesX(i) = velocitiesX(i) * -1;
                positionsX(i) = xWall - (newXposition - xWall);
            else
                positionsX(i) = newXposition;
            end
        elseif newXposition < xWall && positionsX(i) >= xWall
            if newYposition >= holeMax || newYposition <= holeMin
                %disp 'wall bounce';
                velocitiesX(i) = velocitiesX(i) * -1;
                positionsX(i) = xWall + (xWall - newXposition);
            else
                positionsX(i) = newXposition;
            end
        else
            positionsX(i) = newXposition;
        end
        if newYposition < 0
            velocitiesY(i) = velocitiesY(i) * -1;
            positionsY(i) = newYposition * -1;
        elseif newYposition > yMax
            velocitiesY(i) = velocitiesY(i) * -1;
            positionsY(i) = yMax - (newYposition - yMax);
        else
            positionsY(i) = newYposition;
        end
        
        plot(positionsX, positionsY, 'o');
        hold on
        plot(wallX, wallY, 'x')
        xlim([xMin xMax]); ylim([yMin yMax]);
        pause(0.0001)
        hold off
    end
    saveas(gcf, ['gas_',sprintf('%03d',t),'.png'])
end