function [config,NumNodesX] = NewSimulationActin(InitialCondType,VelBoundaryCond,InitialDomainShape)

global LevelSetFunction;
    
%% Physical parameters

% non-dimensional parameters
config.v0                   = 1;
config.mu                   = pi/2;
config.alpha                = 0.5;

% dimensional parameters
config.D                    = 0.2;
config.ksi                  = 1;
config.L                    = 1;
config.M0                   = 1;
config.sigma                = 1;

% actin dynamics parameters
config.zeta                 = 1;
config.epsilon              = 0.01;
config.koff                 = 1;
config.F0                   = 1;

config.M_tot                = config.M0 * config.mu * config.L^2;
config.eta                  = config.alpha * (config.ksi * config.L^2);

%% Boundary velocity (kinematics) parameters
% dimensional parameters
config.VelBoundaryCond      = VelBoundaryCond;
config.V0                   = config.v0 * config.D * config.M0 / config.L;
config.K                    = 0.3;

% other non-dimensional parameters
config.beta                 = (config.sigma*config.M0) / (config.ksi*config.D);
config.k                    = config.K * config.L^3/config.D;

%% Geometric and mesh parameters     
config.InitialRadius = config.L;

% computational domain (big box containing the cell during the simulation)
DomainSizeX                 = 6*config.InitialRadius;
config.Domain               = {[-DomainSizeX/3 2*DomainSizeX/3],[-DomainSizeX/2 DomainSizeX/2]};

% for use in main function
config.DomainSizeX = DomainSizeX;

% Number of Mesh points in x direction (uniform mesh)
% This can be a vector for multiple runs (sequential),
% e.g. NumNodesX = [66 77 88];
NumNodesX                   = 96;

%% Initial domain shape

config.InitDomainShape      = InitialDomainShape;

switch config.InitDomainShape
  
  case 'Circle'
    InitArea                = pi*config.InitialRadius^2;
    LevelSetFunction        = @(x,y) sqrt(x^2 + y^2)- config.InitialRadius;
  
  case 'ShapeFromSavedData'
    
    % load data file from which we extract the initial domain shape at approximate time
    % SelectedTime, the coordinates are changed so that the centroid of the cell is
    % always at (0,0)
    FileNameInitShape       = 'my_data.mat';
    SelectedTime            = 20;
    
    [InitFront,~,InitArea,~,~] = ExtractInitDomainShape(FileNameInitShape,SelectedTime,1,0,0);
    
    % create a global level-set function using 'InitFront'
    LevelSetFunction        = @(x,y) p_poly_dist_my(x,y,InitFront.Polygon(:,1), InitFront.Polygon(:,2));
    config.InitFront        = InitFront;
end
config.InitArea             = InitArea;
pause(2);

%% Initial conditions for actin velocity and myosin

% In case of initial (@ t=0) Myosin gradient, it's defined as
% M(x,y) = k*(1-GradScale*x/MaxDomainX)  OR  M(x,y) = k*(1-GradScale*y/MaxDomainY),
% with 'k' such that:    Integral[M(x,y) dx dy] = M_tot. For cells with
% transformed coordinate, (x_c,y_c) = (0,0), or for symmetric shapes, k = M_tot/InitArea
%
% For uniform Myosin, M(x,y) = M_tot/InitArea

config.InitGradDirection    = 'x';
config.InitGradScale        = 1;

% actin velocity init. values
InitCondition.Ux            = @(x,y) 0;
InitCondition.Uy            = @(x,y) 0;

% actin concentration init. value
InitCondition.F             = @(x,y) config.F0;
 
% myosin init. value
switch config.InitDomainShape
  
  case 'Circle'
    
    switch config.InitGradDirection
      case 'x'
        MaxDomainX = 1;
        InitCondition.M     = @(x,y) config.M_tot/InitArea*(1-config.InitGradScale*x/MaxDomainX);
      case 'y'
        MaxDomainY = 1;
        InitCondition.M     = @(x,y) config.M_tot/InitArea*(1-config.InitGradScale*y/MaxDomainY);
    end  
    
  case 'ShapeFromSavedData'
    
    % myosin init. value
    switch config.InitGradDirection
      case 'x'
        MaxDomainX = max(InitFront.Polygon(:,1));
        InitCondition.M     = @(x,y) config.M_tot/InitArea*(1-config.InitGradScale*x/MaxDomainX);
      case 'y'
        MaxDomainY = max(InitFront.Polygon(:,2));
        InitCondition.M     = @(x,y) config.M_tot/InitArea*(1-config.InitGradScale*y/MaxDomainY);
    end    
    
end
config.InitCondition        = InitCondition;
config.InitialCondType      = InitialCondType;    

end