% MS&E 319 Project 2
% By Tyler Zellmer

% Clear the workspace
clc;
clear;
close all;
tic;

% What fraction of B particles should our space have?
b_fraction = .01;

% Create a 3-D space of particles of dimension s x s x s
s = 10;
numB = 0;

% In x-direction
for l=1:s
    % In y-direction
    for k=1:s
        % In z-direction
        for j=1:s
            % Generate a random number between 0 and 1
            random_number = rand(1,1);

            % If the random number is greater than our b_fraction, then our
            % current point is an A particle. Otherwise, it's a B particle.
            % This ensures that, overall, approximately b_fraction of our
            % space will be B particles
            if random_number > b_fraction
                space(j,k,l) = 0;
            else
                space(j,k,l) = 1;
                numB = numB + 1;
                % Let's also record the location of each B particle
                loc(:,numB) = [j k l]';
            end
        end
    end
end

% Now, recurse through each B point and find nearest neighboring B particle
for a = 1:length(loc)
    difference(a) = 10*s;
    % Recurse through all B-particles in the space and measure their
    % from the current B-particle
    for b = 1:length(loc)
        % Find this particle's distance from the particle of interest. If
        % it is closer than the previous closest, update difference to
        % reflect that.
        temp_length = norm(loc(:,b) - loc(:,a));
        if (temp_length < difference(a) && temp_length ~= 0)
            difference(a) = temp_length;
        end
    end
end

time = toc;
% Display useful information
disp('Number of B-particles: ');
disp(numB);
disp('Max difference (units): ');
disp(max(difference));
disp('Least difference (units): ');
disp(min(difference));
disp('Mean difference (units): ');
disp(mean(difference))
disp('Execution time (s): ');
disp(time);

% Plot location of B-particles
hold on;
for i=1:length(loc)
    plot3(loc(1,i),loc(2,i),loc(3,i), 'o', 'linewidth', 2);
    grid on;
end


% In the x-direction
%     for l=1:s
%         % In y-direction
%         for k=1:s
%             % In z-direction
%             for j=1:s
%                 % If the point at this location is a  B particle, measure
%                 % the distance from
%                 if space(j,k,l) == 1
%                     temp_length = norm([j k l] - loc(a));
%                     if temp_length < difference(a) && temp_length ~= 0
%                         difference(a) = temp_length;
%                     end
%                 end
%             end
%         end
%     end