%% Output of code
% A .mat file comprising xyHoney - a k x 3 array where k denotes the number
% of granules, 3 cols denotes x-coord of center, y-coord of center, radius
% To reduce time taken to run, reduce iterates to 1e8 (however, local
% ordering will be more)

% Run the entire script file directly on HPC; reduce iterations to run on
% home PC

%% Generate ordered arr
% Generate ordered (cubic close packing in 2D) array of granules
filenametrial_ID = 'SomeTrialID';

OccFracLimit = 0.7;
Radius_granules = 0.005;
xLimit = 1;
yLimit = 0.45;
gap_x = 1e-6;
gap_y = 1e-6;
xyZoneGranules = [0 0 xLimit yLimit];
xyHoney = [];
N_iandj = floor(sqrt(0.7854/(pi*Radius_granules*Radius_granules)));
% 0.7854 is the packing limit in 2D

for i = 1:N_iandj
    for j = 1:N_iandj
        [OccFrac, ~, ~] = AreaFraction(xyHoney, {xyZoneGranules}, false);
        if OccFrac < OccFracLimit
            k_x = (j-1)*2*Radius_granules + (j-1)*gap_x + Radius_granules;
            k_y = (i-1)*2*Radius_granules + (i-1)*gap_y + Radius_granules;
            if k_x >= Radius_granules && k_x <= xLimit-Radius_granules ...
                    && k_y >= Radius_granules && k_y <= yLimit-Radius_granules
                xyHoney = cat(1,xyHoney,[k_x k_y Radius_granules]);
            end
        end
    end
end
[~, ~, ~] = AreaFraction(xyHoney, {xyZoneGranules}, true);

fprintf('\n********\nBegin spreading of granules\n********\n')

%% Spread granules apart
% Incrementally move granules apart
tic
N = size(xyHoney,1);
for i = 1:1e9
    rand_granule_index = randi(N);
    rand_granule = xyHoney(rand_granule_index,:);
    try_x = rand_granule(1) + 0.01*(2*rand(1,1)-1);
    try_y = rand_granule(2) + 0.01*(2*rand(1,1)-1);
    xyHoney_without_rand_granule = xyHoney;
    xyHoney_without_rand_granule(rand_granule_index,:) = [];
    dist_xy = sqrt((try_x-xyHoney_without_rand_granule(:,1)).^2 +...
        (try_y - xyHoney_without_rand_granule(:,2)).^2);
    dist_xy = min(dist_xy);
    if dist_xy >= 2*Radius_granules && ...
            MinDistToWall(try_x, try_y, {xyZoneGranules}) ...
            > Radius_granules 
        xyHoney(rand_granule_index,:) = [try_x try_y Radius_granules];
    end
end
toc

%% Save xyHoney
rng('shuffle') % Unset seed = 1, shuffle seed
s='A':'Z';
filename_code = [s(randi(26)) s(randi(26)) s(randi(26))...
    s(randi(26)) s(randi(26)) s(randi(26))];
disp(filename_code)
clear s
filename = strcat(date,'_',filenametrial_ID,...
    '_','OccFracLimit_',num2str(OccFracLimit),...
    '_','Radius_granules_',num2str(Radius_granules),...
    '_',filename_code,'_granules-config_.mat');
save(filename,'xyHoney')

