%% Set seed, simulation parameters
rng(1)

tic % Start stopwatch timer to time the simulation

ScriptVar = containers.Map({'N_Particles', 'xLimit', 'yLimit', 'Sensing',...
    'number_of_cusps', 'init_final_StepLim', 'MobilityForEachParticleType',...
    'filenametrial_ID'},...
    {[10000], 1, 0.45, false, 0.1, [0.04 0.04], [1.0 1.0],...
    'SomeFileNameID'});

ScriptVar('StepLim_functional_form') = 'default';
% ScriptVar - stores simulation parameters
% N_Particles - Number of particles of each type; single-type used
% xLimit - Box length
% yLimit - Box width 
% init_final_StepLim - [a b] - a = StepLim at x=0, b = StepLim at x=xLimit
% StepLim denotes the max step size that can be sampled in any given step
% for a particle
% number_of_cusps - number of sharp discontinuities in going from StepLim = a 
%   to StepLim = b from x=0 to x=xLimit; value not relevant if
%   StepLim_functional_form is 'linear' and not 'stepwise'
% MobilityForEachParticleType - 1xk array where k = number of particle
%   types; [1 0.5 0.2] => basal StepLim for type 1 is determined from 
%   init_final_StepLim but basal StepLim for type 2 is halved and basal StepLim
%   for type 3 is one-fifth of init_final_StepLim value respectively. 
%   For all simulations conducted, [1.0 1.0] was used
% filenametrial_ID - Code used to identify this simulation run (for
% reference)
% StepLim_functional_form = 'default' if StepLim is same throughout the
%   box, 'linear' if it increases linearly from left to right, 'stepwise' if
%   increase in a discrete fashion from left to right. In the case of
%   stepwise, number_of_cusps >= 1

if ScriptVar('Sensing')
    ScriptVar('SensingStrength') = 4; % k x k matrix for k particle types 
    ScriptVar('rSense') = 0.02; % k x 1 matrix for k particle types
    % modified StepLim = initial StepLim * exp(-SensingStrength * number of
    % neighbours lying within rSense um away from particle in question
end
% Sensing is used to model interactions; see MakeStep.m for more info

%% Defining zones for particle initialisation
ScriptVar('xyStartFigure') = {[0 0 ScriptVar('xLimit') ScriptVar('yLimit')]}; 
% 1 x k cell where k denotes the number of particle types
% [0 0 xLimit yLimit] - 0, 0 specify the bottom-left corner of the box,
% xLimit denotes the length, yLimit denotes the breadth
if length(ScriptVar('xyStartFigure')) ~= length(ScriptVar('N_Particles'))
    error('Error. The number of arrays in xyStartFigure does not match the number of particle types. For each particle type, specify a separate initialisation zone as an array element in the xyStartFigure cell.')
end

%% Setting inhomogenous diffusivity - choose option A or option B (below)
% The variable viscous_granules is to be set as true if Option A is
% chosen, false if not.

%% Defining viscous zones - Option A - viscous granules
viscous_granules = false;
xLimit = ScriptVar('xLimit');
yLimit = ScriptVar('yLimit');
xyZoneGranules = {[0 0 xLimit yLimit]}; 
% Initialize all granules within xyZoneGranules{1}
if viscous_granules
    Radius_granules = 0.05; % Granule radius
    OccFracLimit = 0.4;     % Granule packing density limit
    for i = 1 % loop meant only for wrapping code
        xyZoneGranules = xyZoneGranules{1};
        xyHoney = zeros(1,3); 
        % After the while loop below, xyHoney would be a k x 3 matrix 
        % where k denotes the number of granules, columns 1-3 contain the 
        % x,y co-ordinates of the granule centers, radius
        xyHoney(:,3) = Radius_granules;
        xyHoney(1,1) = xyZoneGranules(1) + xyZoneGranules(3)*rand(1,1);
        xyHoney(1,2) = xyZoneGranules(2) + xyZoneGranules(4)*rand(1,1);
        [OccFrac, ~, ~] = AreaFraction(xyHoney, {xyZoneGranules}, false);
        % The code below runs (and generates new granules) until the 
        % packing density reaches OccFracLimit
        while OccFrac < OccFracLimit % While current packing density is 
            % less than desired packing density, keep adding granules
            try_x = xyZoneGranules(1) + xyZoneGranules(3)*rand(1,1);
            try_y = xyZoneGranules(2) + xyZoneGranules(4)*rand(1,1);
            dist_xy = sqrt((try_x-xyHoney(:,1)).^2 + (try_y - xyHoney(:,2)).^2);
            dist_xy = min(dist_xy);
            if dist_xy >= 2*Radius_granules && ...
                    MinDistToWall(try_x, try_y, {xyZoneGranules}) ...
                    > Radius_granules && max(dist_xy) <= 2.001*Radius_granules
                % Checks if the new granule is non-overlapping with all
                % others, if it lies entirely inside the box, and if it is
                % sufficiently close to the old granules
                % Third condition is necessary only for packing up to
                % OccFracLimit of 0.6; for packing density beyond 0.6,
                % refer to GenerateGranules.m
                xyHoney = cat(1, xyHoney, [try_x try_y Radius_granules]);
                [OccFrac, ~, ~] = AreaFraction(xyHoney, {xyZoneGranules}, false);
            end
        end
    end 
    xyHoney = mat2cell(xyHoney,ones(1,size(xyHoney,1)));
    [~,~,~] = AreaFraction(xyHoney, {xyZoneGranules}, true); % print area fractions
    clear dist_xy try_x try_y

    % Add parameters, granules to ScriptVar
    ScriptVar('xyZoneGranules') = xyZoneGranules;
    ScriptVar('Radius_granules') = Radius_granules;
    ScriptVar('OccFrac') = OccFrac;
    ScriptVar('OccFracLimit') = OccFracLimit;
    ScriptVar('xyHoney') = xyHoney;
end
clear OccFrac OccFracLimit Radius_granules xyHoney xyZoneGranules % As stored within ScriptVar
clear i xLimit yLimit viscous_granules

%% Defining viscous zones - Option B - other box, circ zones
% Comment the code below if granules are used instead

xyHoney = {[0 0 ScriptVar('xLimit')*0.2 ScriptVar('yLimit')]};
% Cell containing all the viscous zone(s)
ScriptVar('xyHoney') = xyHoney;
clear xyHoney
%%%%%%%%%%%%% IMPORTANT NOTE %%%%%%%%%%%%%
% If viscous zone(s) are not required, set N_FramesRegion = 0 (see section below)

%% Number of frames/timesteps/epochs in the run
ScriptVar('N_Frames') = 1000;            % Total number of frames
ScriptVar('N_FramesStartBoundaries') = 1; % For initialisation zone(s); Set >= 1 
ScriptVar('N_FramesRegion') = 1000;      % For viscous zone(s) 
ScriptVar('xyHoneyRatio') = sqrt(1/4);   % Ratio of StepLim inside, outside viscous zones
% Note that N_FramesStartBoundaries overlaps with N_FramesRegion and both 
% should be <= N_Frames

%% Initialise particles (initialise frame 1)
N_Particles = ScriptVar('N_Particles');
N_Frames = ScriptVar('N_Frames');
AllCoords = zeros(sum(N_Particles),2,N_Frames);
% AllCoords = Coords of all particles for all frames
% AllCoords, ScriptVar will be saved at the end of the run
xyStartFigure = ScriptVar('xyStartFigure');
AllCoordsInit = [];
for i = 1:length(N_Particles)
    AllCoordsInit = cat(1, AllCoordsInit,...
        InitialiseParticles(ScriptVar, N_Particles(i), xyStartFigure(i)));
end
clear N_Particles N_Frames xyStartFigure InitCoordsCell 
% Cleared as stored within ScriptVar
clear i 

%% Run simulation
N_Frames = ScriptVar('N_Frames');
% Get Frame 1
AllCoords(:,:,1) = NextSetCoords(ScriptVar, 1, AllCoordsInit); 
% Run (Frames 2 to N_Frames)
for frame_i = 2:N_Frames
    AllCoords(:,:,frame_i) = NextSetCoords(ScriptVar, frame_i,...
        AllCoords(:,:,frame_i-1));
    if ~mod(frame_i,100)
        fprintf("Frame %d done\n",frame_i) 
        % Print frame, elapsed time in case of long-running simulations
        toc
    end
end
toc % End stopwatch timer
clear N_Frames frame_i AllCoordsInit

%% Save simulation parameters, coords data
% Create a six-letter code to uniquely identify this run, besides 
% filenametrial_ID, to avoid overwriting while running the same script
% multiple times if required
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
rng(1) % Set seed back to 1

filenametrial_ID = ScriptVar('filenametrial_ID');
filename = strcat(date,'_',filenametrial_ID,'_reg-box_',filename_code,'_table-params.mat');
%%%%%%%%%%%%% IMPORTANT NOTE %%%%%%%%%%%%%
% In case of small AllCoords, use:
save(filename,'AllCoords','ScriptVar')
% In case of a potentially large AllCoords variable, use:
% save(filename,'AllCoords','ScriptVar', '-v7.3')

clear filename filename_code filenametrial_ID