function results=main_Dicty(p, ICs)

%{
Runs the hybrid PDE-AB-Model for Dictystelium dynamics

*Input*
p .. structure containing all parameters
ICs (optional) .. structure containing pre-set initial conditions)

*Output*
results .. structure containing the saved simulation output

*Notes*
Depending on the simulation settings in p, all simulation files, a video,
etc. are also saved (see below)

--

Variable name notes

*Cell characteristics*

xPos .. Nx1 vector x positions of cells
yPos .. Nx1 vector y positions of cells
cellDir .. Nx1 vector of cell orientation (in radians)
cellSpeed .. Nx1 vector of cell speeds

actConc .. Nx1 vector concentration of activator
inhConc .. Nx1 vector concentration of inhibitor
chemSens .. Nx1 vector chemotactic sensitivity

*External cAMP dynamics*

nX .. number of gridpoints in x
nY .. number of gridpoints in y

sigConc .. nX x nY matrix with signal (cAMP) concentraion

Written 12/2021 by Angelika Manhart
%}

%%%%%%%%%%
%% Prelude
%%%%%%%%%%

%% Set-up spatial grid 


if p.optBC==0 % periodoc BCs
    
    % calculate grid spacing
    p.dx = p.Lx/p.nX;
    p.dy = p.Ly/p.nY;
    
    p.xVals = linspace(-p.Lx/2,p.Lx/2-p.dx,p.nX);
    p.yVals = linspace(-p.Ly/2,p.Ly/2-p.dy,p.nY);
    
elseif p.optBC==1 % no flux BCs
  
    % calculate grid spacing
    p.dx = p.Lx/(p.nX-1);
    p.dy = p.Ly/(p.nY-1);    

    p.xVals = linspace(-p.Lx/2,p.Lx/2,p.nX);
    p.yVals = linspace(-p.Ly/2,p.Ly/2,p.nY);
    
end

%% Initial conditions

if nargin==1 % if ICs are not provided .. 
    
    % define initial conditions
    [sigConc, xPos, yPos, actConc, inhConc, chemSens, cellDir, cellSpeed]=get_IC(p);    
   
    
else % if the ICs are provided ..
        
    % .. read them in
    sigConc=ICs.sigConc;
    xPos=ICs.cellPos(:,1);
    yPos=ICs.cellPos(:,2);
    actConc=ICs.cellChar(:,1);
    inhConc=ICs.cellChar(:,2);
    chemSens=ICs.cellChar(:,3);
    cellSpeed=ICs.cellChar(:,4);
    cellDir=ICs.cellChar(:,5);
    
    % Activate these code lines only if you are creating a de-nuovo spiral
    
%     sigConc(1:p.nX/2,:)=0;
%     inhConc(yPos<0)=2.5;
    
    % .. and clear
    clear ICs
end

%% Plot ICs

% define signal gradient as zero for initial timestep
dxSigConc=zeros(p.nX, p.nY);
dySigConc=zeros(p.nX, p.nY);

make_plot(sigConc, xPos, yPos, actConc, inhConc, chemSens, cellDir, cellSpeed, dxSigConc, dySigConc, p,0);

pause(0.1)

%% Make Results folder and save files

% save initial timestep
results(1).sigConc=sigConc;
results(1).cellPos=[xPos yPos];
results(1).cellChar=[actConc inhConc chemSens cellSpeed cellDir];
results(1).time=0;

if p.saveEvery>0 || p.isVideo==1
    
    % make new folder for results
    if isfield(p, 'folderName') % if the folder name is provided use that ..
        folderName=p.folderName;
    else % .. otherwise create folder using current time stamp (avoids writing over previous simulations)
        folderName=char(datetime('now', 'Format', 'dMMMy_HH_mm_ss'));
    end
    
    mkdir(strcat('Results/',folderName));
    
    % save the program files in there (for reproducability)
    copyfile('Program_Files', strcat('Results/',folderName,'/anc_Program_Files'))
    
    % save the parameter files in there
    save(strcat('Results/',folderName,'/p.mat'), 'p')
    
    % save IC pic
    print(gcf, strcat('Results/',folderName,'/pic_', num2str(0)), '-dpng')
   
    % set-up video
    if p.isVideo==1
        aviobj = VideoWriter([strcat('Results/',folderName,'/Movie.avi')]);
        aviobj.Quality=40;
        open(aviobj)
        
        % add first frame
        M=getframe(gcf);
        writeVideo(aviobj,M);
    end
    
end 


%%%%%%%%%%%%
%% Time Loop
%%%%%%%%%%%%
time = 0;
n=0;
nn=1;

% loop 
while time<p.maxTime
    n=n+1;
    tic
    
    %% STEP 0 - change parameters dynamically in time (optional: could be changed to use other parameters)    
    
    if p.timeDepParas==1 % only do this if this parameter is set to 1
            
        initVal=1/6; % intitial value of parameter
        finalVal=1/15; % final value of parameter
        
        % linear function between them
        p.intTimeScale=finalVal+(time-p.maxTime)*(finalVal-initVal)/p.maxTime;
        
    end
    
    %% STEP 1 - Update cell internal signalling state variables
    %{
    Done by integrating of the interval [0 p.dt]. To make code faster, this
    could be replaced by a one-step solver
    %}
    
    % determine the external signal concentration each cell feels
    sigConcPerCell=get_sigConcPerCell(sigConc, xPos, yPos,p);
    
    % create initial conditions
    y0=[actConc' inhConc'];
    
    % set time span
    tspan=[0; p.dt];

    % solve using a FitzHugh-Naguma type model (see RHS_modFHN)
    [~,y] = ode45(@(t,y)RHS_modFHN(t,y,sigConcPerCell,p),tspan,y0);

    % extract new values
    actConc=y(end,1:p.N)';
    inhConc=y(end,p.N+1:end)';

    % Apply heaviside function to activator concentration
    actConcUse=actConc;
    actConcUse(actConc<=0)=0;
    actConcUse(actConc>0)=1;
    
    %% STEP 2 - Update signal concentrations
    
    % solves the signal concentration equation
    sigConc=get_sigConc_p1(sigConc, actConcUse, xPos, yPos, p);
        
    %% STEP 3 - Update cell positions
    
    % determine repulsion force in case there is direct cell-cell
    % interaction
    if p.repStrength==0 % if no cell-cell interaction ..
        
        repForce=zeros(p.N,2); % .. set to zero       
        
    else % .. otherwise, calculate repulsion effect  
        
        repForce=get_repForce(xPos, yPos, p);
    end
    
    
    if p.speedModelP1~=0 % if cells can move ..
        
        % .. update their position, direction, speed and chemotactic
        % sensitivity
        [xPos, yPos, cellDir, cellSpeed, chemSens]=...
            get_cellMove_p1(xPos, yPos, cellDir, cellSpeed, chemSens, sigConc, sigConcPerCell, repForce, p);
    end
    
    % map back into domain
    if p.optBC==0
        xPos=mod(xPos+p.Lx/2,p.Lx)-p.Lx/2;
        yPos=mod(yPos+p.Ly/2,p.Ly)-p.Ly/2;
    else
        xPos(xPos>p.Lx/2)=p.Lx/2; xPos(xPos<-p.Lx/2)=-p.Lx/2;
        yPos(yPos>p.Ly/2)=p.Ly/2; yPos(yPos<-p.Ly/2)=-p.Ly/2;
    end
    
    %% STEP  4 - check, plot & store
   
    % update time
    time = time+p.dt;
    
    % are there problems?
    if isreal(sigConc)==0
        disp('Complex Values!!')
%         break
    end
    
    if min(sigConc(:))<0
        disp('Negative signal conc!!')
%         break
    end
    
    % plot
    if mod(n, p.plotEvery)==0
        % plot
        make_plot(sigConc, xPos, yPos, actConc, inhConc, chemSens, cellDir, cellSpeed, dxSigConc, dySigConc, p,time);
        pause(0.001)
        
        % save video frame
        if p.isVideo==1
            M=getframe(gcf);
            writeVideo(aviobj,M);
        end

    end
    
    % Store    
    if mod(n, p.saveEvery)==0
        nn=nn+1;
        % find better solution for storing that doesn't need size changes
        % to results structure
        results(nn).sigConc=sigConc;
        results(nn).cellPos=[xPos yPos];
        results(nn).cellChar=[actConc inhConc chemSens cellSpeed cellDir];
        results(nn).time=time;
    end
  
    toc
   
end
% END OF TIME LOOP

%%%%%%%%%%%%%%%
%% Final things
%%%%%%%%%%%%%%%%

% close video file
if p.isVideo==1
    close(aviobj)
end

if p.saveEvery>0 % .. if several time steps have been saved
    
   % save all results. Note: This might be a problem if the results file is too large! 
   save(strcat('Results/',folderName,'/results.mat'), 'results','p', '-v7.3')

else % .. otherwise, just store the final state in results
   results(2).sigConc=sigConc;
   results(2).cellPos=[xPos yPos];
   results(2).cellChar=[actConc inhConc chemSens cellSpeed cellDir];
   results(2).time=0;

end

% save last picture
if p.saveEvery>0 || p.isVideo==1
    
    % make a pic of last timestep
    print(gcf, strcat('Results/',folderName,'/pic_', num2str(n)), '-dpng')
    
    % save final result
    save(strcat('Results/',folderName,'/results.mat'), 'results','p')
    
end