%% Master file for alternative Myosin transport supplementary results in:
% A free-boundary model of a motile cell explains turning behavior, by
% Masoud Nickaeen et al. 2017
%
% Tested and developed using MATLAB R2016b, the codes run in both
% windows and linux environments (64-bit).
% 
% Set appropriate parameters and choose methods from the list of available
% options. The results will be saved to the 'myosin' folder in the current
% working DIR.
%
% If you want to "load solution from a saved file" AND at the same time
% "change some parameters", make sure to put new parameters in this
% file by changing 'config' data structure. For example, if you
% want to load data and change V0 and M_total, you need to add these lines
% to the '%% Final changes to model parameters' block:
%     config.v0                   = 7.5; 
%     config.mu                   = 0.25*pi; 
%     config.M_tot                = config.M0 * config.mu * config.L^2;
%     config.V0                   = config.v0 * config.D * config.M0 / config.L;
%
% Note the Matlab code is based on the dimensional form of the equations.
% Therefore, varying any dimensionless parameter should be accompanied by
% further correcting for the corresponding dimensional parameters.

clear all;
clc;
format long;
close all;
 

%% List of available options for running the code

VelBoundaryCond             = {'ZeroDirichlet';'ZeroStress'};              % 'ZeroDirichlet' is the ZV model, and 'ZeroFlux' is the ZS model
MyosinBoundaryCond          = {'ZeroFlux','Flux'};                         % 'Flux' sets the free myosin diffusion flux equal to the bound myosin advection flux at the boundaries
FrontPropagationMethod      = {'first_order';'fourth_order'};
FrontRedistributionMethod   = {'DEFAULT';'EXPANSION_REDIST';'NO_REDIST'};
TimeDiscrSchemes            = {'Stokes'};
DiffusionCoefApproxMethod   = {'HarmonicAverage';'ArithmeticAverage'};
InitialDomainShape          = {'Circle','ShapeFromSavedData'};
InitialConditionsType       = {'NewSimulation','SolutionFromSavedData'};


%% Initial Condition type Selection

InitialCondType             = InitialConditionsType{1};

switch InitialCondType
  
  case 'SolutionFromSavedData'
    FileName = 'my_data.mat';
    SelectedTime       = 20;
    [config,NumNodesX] = LoadDataBFMyosin(InitialCondType,FileName,SelectedTime);
  case 'NewSimulation'
    [config,NumNodesX] = NewSimulationBFMyosin(InitialCondType,VelBoundaryCond{2},InitialDomainShape{1});
    SelectedTime       = 0;
    FileName           = '';
end
config.FileName             = FileName;
config.SelectedTime         = SelectedTime;
config.MyosinBoundaryCond   = MyosinBoundaryCond{2};


%% Final changes to model parameters
% 
% config.v0                   = 0.1; 
% config.mu                   = 0.9572; 
% config.M_tot                = config.M0 * config.mu * config.L^2;
% config.V0                   = config.v0 * config.D * config.M0 / config.L;

%% Time discretization parameters  dt = TimeStepScale * MeshSize/(10*MaxVelocity), MaxVelocity=1

config.T_max                = 5;
config.TimeDiscrScheme      = TimeDiscrSchemes{1};
config.TimeStepScale        = 1;   

%% Segregated solver control parameters

config.FixedPointIterationsMax  = 35;
config.FixedPointIterationsTol  = 1e-10;
config.FixedPointStopCrite      = 1;                                       % 1 Absolute solution error
config.DisplayFixedPointError   = false;
config.SolveMyosinsInOneStep    = true;                                    % 'true' solves bound and free myosins in one step at each segregated solver iteration


%% avoiding singularity solution

% by advection and diffusion terms in myosin equation:
%     u = u*(1-M_b/MaxMyosinAdvection), D = D*(1-M_u/MaxMyosinDiffusion)
%
% Note that for comparison against the one-version transport equation, i.e.
% main text results, the thresholds for anti-crowding terms need to be
% modified by reaction kinetics.
config.AvoidSingbyAdvec     = true;
config.AvoidSingbyDiffu     = true;
config.MaxMyosinAdvection   = 15  * config.k_on/(config.k_off + config.k_on);
config.MaxMyosinDiffusion   = 125 * config.k_off/(config.k_off + config.k_on);

if config.AvoidSingbyDiffu
  config.DiffCoefApproxMeth = DiffusionCoefApproxMethod{1};
end


%% avoiding collapsing cells [ Vp = V0/M0*(A0/A) - K.(A-A0(A0/A)^n)   OR  V0/(M0+M)*(A0/A) - K.(A-A0(A0/A)^n) ]

config.A0                   = pi;
config.AreaExponent         = 2;

% fixed boundary problem, i.e. v|_f = 0
config.FixedBoundary        = false;


%% Other settings

config.returnAllData        = true;
config.NumFrames            = 300;
config.Visualization        = false;
config.FrontFormula         = [];
config.smoothing            = false;
config.bilinearInterp       = true;
config.FrontRefinementScale = 1;
config.FrontPropagate       = FrontPropagationMethod{1};
config.FrontRedistMethod    = FrontRedistributionMethod{1};
config.FrontRedistFrequency = 5;
config.EdgeLength           = sqrt(pi);
config.SaveDuringSolve      = false;
config.CreateMovie          = true;
config.CreateAreaCentPlot   = true;
config.MaxMinMyosin         = true;


%% Membrane velocity (implemented in an m-file that is passed to FronTier)

switch config.VelBoundaryCond
    case 'ZeroStress'
      if config.bilinearInterp
        copyfile('FT_velo_func_NeumannBC_bilinear.m', 'FT_velo_func.m', 'f');
      else
        copyfile('FT_velo_func_NeumannBC.m', 'FT_velo_func.m', 'f');
      end
    case 'ZeroDirichlet'
        copyfile('FT_velo_func_DirichletBC.m', 'FT_velo_func.m', 'f');
end 

if config.FixedBoundary
    copyfile('FT_velo_func_FixedBoundary.m', 'FT_velo_func.m', 'f');
end

pause(2);


%% Get the solution

for i = 1:length(NumNodesX)
  
  config.NumNodesX         = NumNodesX(i);
  config.MeshSize          = config.DomainSizeX/config.NumNodesX;
  
  switch config.TimeDiscrScheme
    case 'Stokes'
      disp ('Stokes-type Velocity Equation - Semi-Implicit Time Discretization')
      [Ux_All, Uy_All, M_All_b, M_All_u, Time_All, PosX, PosY, Front_All, ...
        FrontVelocity_All, Global_Area_All,Area_All,InitialMass,str] = ...
        GetSolutionStokesBFMyosin(config,SelectedTime); 
  end
  
end  

%% Post-processing

% plotting area, centroid of the cell
if config.CreateAreaCentPlot
  figure(1);
  [x_c,y_c,r_c,Area_Green] = GetAreaCentroid(Front_All);
  subplot(2,2,1)
  plot(Time_All',x_c,'-bo');legend('x_c')
  hold on;
  subplot(2,2,2)
  plot(Time_All',y_c,'-ro');legend('y_c')
  hold on;
  subplot(2,2,3)
  plot(Time_All',r_c,'-ko');legend('r_c')
  hold on;
  subplot(2,2,4)
  plot(Time_All',Area_Green,'-go');
  hold on;
  plot(Time_All,Area_All,'-b.');
  legend('Area_{Green}','Area_{Voronoi}') 
end

% plotting minimum and maximum values of myosin
if config.MaxMinMyosin
  
  FontSize = 14;
  figure(2);
  M_All = M_All_b+M_All_u;
  max_Myosin = max(M_All,[],1);
  plot(Time_All,max_Myosin,'ko')
  hold all
  min_Myosin = 0*max_Myosin;
  for i=1:size(M_All+M_All,2)
    min_Myosin(i) = min(M_All(M_All(:,i)~=0,i));
  end
  plot(Time_All,min_Myosin,'ro')
  legend('Max. Myosin','Min. Myosin')
  axis square
  xlabel('Time (t)','FontSize', FontSize)
  ylabel('Myosin Concentration','FontSize', FontSize)
  set(gca,'FontSize', FontSize)

end

% create movie
if config.CreateMovie
    
    SkipTimeFrames = floor(config.NumFrames/20);
    XYRatio = diff(config.Domain{1})/diff(config.Domain{2});
    FigurePosistion = [100 100 round(XYRatio*1000) 1000];
    XLIM = config.Domain{1};
    YLIM = config.Domain{2};
    XTicks = round(XLIM(1):4:XLIM(2));
    YTicks = round(YLIM(1):4:YLIM(2));
    FontSize = 18;
    FontName = 'Times new roman';
    QuiverSpec = {'color','w','LineWidth',1};
    VelocityScale = config.k_on/(config.k_on+config.k_off);
    CreateMovie(M_All_b+M_All_u, VelocityScale*Ux_All, VelocityScale*Uy_All, PosX, PosY, Time_All, Front_All, str, config.MeshSize,...
      config.v0,config.mu, config.alpha, SkipTimeFrames, FigurePosistion,XLIM,YLIM,FontSize,FontName,XTicks,YTicks,QuiverSpec);

end


