
% 13/10/2021
% Velocity exponentially decaying with force applied
clear all
close all
addpath Functionals:Interpolation_Geometry:Reparametrization:Examples_static_NEW:Matrices
global phi Invariant_aux C_1 dt ini_Val Material Constraints V CC  C_psi 


input_file ='Examples_Dynamic/OM_Support.mat';
load(input_file)

MaxFunEvals=3000;
MaxIter=50000;
tol_X=1e-3;

options = optimoptions('fmincon','Algorithm','interior-point' ,'Display','final', 'CheckGradient', false, ...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,'MaxFunctionEvaluations',MaxFunEvals);

options = optimoptions(options,'Algorithm','interior-point',"EnableFeasibilityMode",true,...
        'TolX',1e-7,'Display','final');


iopen=ini_Val.iopen;

% Do you want to impose the fixed volume with support? 0: No 1: Yes
Constraints.iVol_s = 0;
Constraints.Pr_S = -lambda.eqnonlin(1);
Material = rmfield(Material,'h');
Aeq = [Aeq; zeros(1,length(Aeq))];
Aeq(end,end) = 1;
beq(end+1)   = x(end);
%How many steps prior to final shape you wish to start
% S_back =1;
% x = x_hist(:,end-S_back+1);

%%

%C_phi=C_phi_hist(:,end-S_back);

%Constraints.C_s=C_s_eq(end-S_back);
dt = [];
time = [];
max_dau=10;

%% Parameters
%Membrane viscosity
Material.mu = 3.65e-3;
%Cortex viscosity (of a thin layer)
Material.eta_c = 5e3 / 1000 / 10;
%Membrane support friction
Material.eta_s = 2.4e-1 * 5*0;
%Parameters for z-cutoff of friction.
Material.tmp_s = 10;
Material.h_cut_off = 1;
% Force at edge of membrane
Constraints.f0 = 0;

% EFFECTIVE VISCOSITY/ELASTICITY OF THE REST OF THE SYSTEM
Material.eta_END = Material.eta_s * pi * x(end-1)^2 * log(100);
Material.xi_Cor_perp = 0;
Material.K_outer = 1 * 0;

% This is membrane BAR protein
ap = 3;
Material.Chi=0;
Material.RT_a0 = 1/20 / ap;
Material.phi_max=0.35;
Constraints.C0=-0.03  ;
Diff = 10; % this is 1 um^2/s
Diff = 1; % this is 0.1 um^2/s
Material.eta_a0 =  1/20 / Diff / ap;
Material.phi_bar = 0.05;
phi_bar = 0*Material.phi_bar;



% This is regulator
Material.D = 1 / 5000 ;
Material.k1= 1 / 5000;
Material.k2= 1 / 5000;
Material.enrichment_treshold = 2;
% Coefficients that transform psi into contractility as f1 - f2*psi
Material.f1 = 2;
Material.f2 = 0.7;



%Material.eta_c = Material.eta_c * 1e12;
%Constraints.C0 = 0;

%% To fix the height
% Aeq         = [Aeq; zeros(1,length(Aeq))];
% Aeq(4,2)    = 1;
% beq(4)      = x(2);

%% Times
t_f = 1000 * 30 * 3;
%t_f = 1000 * 15
time = [0];
x_hist = x;
dt = 0.001 * 100 ;
i_max = 3000;
dt_max = 10 * 20 * 3;
dt_max = 10 * 10;


%% Normal code

ini_Val.UU =ini_Val.U(4:end-3);

C_phi = phi_bar * ones(size(C_phi));
[phi, d_phi] = interpolate_func_2(Invariant,C_phi);
[phi_s, d_phi_s] = interpolate_func_2(Invariant_sample,C_phi);

C_psi = 0 * ones(size(C_phi));
[psi, d_psi] = interpolate_func_2(Invariant,C_psi);


[C,C1,C2] = interpolate(Invariant,x);
Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2,phi);


[C_s,C1_s,C2_s] = interpolate(Invariant_sample,x);
% Update the local area constraints
Constraints.S0_L = localArea_functional(Aux);
Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
[tS] = Area_functional(Aux);
[tV] = Volume_functional_support(Aux);

figure(2)
hold on
plot(C_s( :,1),C_s( :,2))
axis equal
hold on
plot(x(1:2:end),x(2:2:end),'ko-')

C_phi_hist = [C_phi];
phi_s_hist = [phi_s];
C_psi_hist = [C_psi];
C_s_hist = [C_s];
tot_ener_hist = [];
area_hist = [tS];
Chm_p_hist= [];
force_hist=[Constraints.f0];
force_hist_z=[];
Change_V0 = [];
sigma = [];

x_0=x;
x_1=x;
ct_rep=0;
dt_inc=0;
ini_Val.cc = 0.2;
ini_Val.ce = 1;
i=0;
epsl=1e-4;
ctr=0;
count=0;

%Cortex_Flat_MAB(psi);
V_hist = [V];

% velocity of growth of dendritic network
v0_dendrite = 9e-3;
F0          = 0.5; %0.308
decay_f = 15;
% Specifies whether the dendritic network is in contact with the membrane
contact = 0;
% Change in magnitude of height (just to initialize here)
Ht_Chng = 0;



%figure(2)
%quiver(Invariant.gaussp*Material.R0,Invariant.gaussp*Material.R0*0,v_c_c',v_c_c'*0)
%tmp_z  = 0.5*(1-tanh(Material.tmp_s*(C(:,2)-Material.h_cut_off)));
%v_c_m = v_c_m .* tmp_z;
%quiver(C(:,1),C(:,2),v_c_m,v_c_m*0)

%return

Proj_area_0 = pi * x(end-1)^2;
Areal_Strain = 1;
max_psi_hist = [];
max_phi_hist = [];
max_V_c_hist = [];
tip_hist       = [];
tip_force_hist = [];

Material.R_ini = x(end-1);


k = 0;
%%
while time(end) < t_f & i < i_max & k < 100000
    i = i+1;
    k = k+1;
    fprintf('Instant %10.5f towards %10.5f\n', time(end)+dt, t_f)
        
    %
    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);
    %     Constraints.S0_L = localArea_functional(Aux);
    %     Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
    
    
    max(abs(Aux.dau))
    
    if any(abs(Aux.dau) > max_dau)
        disp('**********************')
        disp('We reparametrize')
        [C_fit,phi_fit] = original_shape_to_fit(x, C_phi); % at points uniformly spaced in arc-length
        %[C_fit] = original_shape_to_fit(x); % at points uniformly spaced in arc-length
        ce=ini_Val.ce; cc=ini_Val.cc;
        %creparametrize(x,ce,cc);
        x_new = reparametrize(C_fit,x);
        x=x_new;
        % Update protein coverage
        C_phi_new = reparametrize_func(phi_fit,C_phi);
        C_phi=C_phi_new;
        [phi, d_phi, dd_phi] = interpolate_func_2(Invariant,C_phi);

        
        % Update the local area constraints
        [C,C1,C2] = interpolate(Invariant,x);
        Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);
        Constraints.S0_L = localArea_functional(Aux);
        Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
        ct_rep=ct_rep+1;
        if(ct_rep==2)
            max_dau=max_dau*1.15;
            if(max_dau>200)
                max_dau = 200;
            end
            ct_rep=0;
        end
    else
        ct_rep=0;
    end
    
    phi_1 = phi;
    C_phi_1 = C_phi;
   
    if 2==1
        % Map phi to the flat cortex
        [N ,EDGES ,BIN ] = histcounts(C(:,1),ini_Val.UU*Material.R0);
        [N_,EDGES_,BIN_] = histcounts(Invariant.gaussp,ini_Val.UU);
        phi_flat_ = zeros(size(Invariant.gaussp));
        for kk = 1:length(ini_Val.UU)-1
            if sum(BIN==kk)>0
                phi_flat(kk) = max(phi(BIN==kk));
            else
                phi_flat(kk) = phi_flat(kk-1);
            end
            phi_flat_(BIN_==kk) = phi_flat(kk);
        end
        
        psi_1 = psi;
        C_psi_1 = C_psi;
        
        Rhs_psi = vector_rhs_psi(psi_1,phi);
        C_psi = zeros(length(Rhs_psi),1);
        [psi, d_psi] = interpolate_func_2(Invariant,C_psi);
    end
    
    
    x_1 = x;
    [C_1,C1_1,C2_1] = interpolate(Invariant,x_1);
    [x, fval, exitflag, output, lambda] = fmincon(@Incremental_Mech_Rayleighian,x,[],[],Aeq,beq,[],[],@Nonlin_constr,options);
    
    if i>0
        if  max(abs(x-x_1)) > 0.5 || count > 5
            disp('*****************')
            disp('We reduce the timestep')
            dt = dt/2;
            if(contact)
                Ht_Chng  = Ht_Chng/2;
                beq(end)   = beq(end) - Ht_Chng;
            end
            pause(0.2)
            x = x_1;
            continue
        end
        if max(abs(x-x_1)) < 5e-2
            dt = dt * 1.5;
            if dt >= dt_max
                dt = dt_max;
                
            else
                disp('*****************')
                disp('We increase the timestep')
            end
        end
        
    end
    
    % To change the height depending on the force
    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2,phi);
    
    %Invariant_mem_to_cortex=precal_generic(C(:,1)/Material.R0);

    %[psi_m]=interpolate_func_2(Invariant_mem_to_cortex,C_psi);
    %[tD_c,dD_dCP_c] = Incremental_cortex_diss(Aux,C,C_1,C1);
     
    if(contact)
        disp(lambda.eqlin(end))
        Ht_Chng = dt*v0_dendrite*(exp(decay_f*lambda.eqlin(end)) - exp(-decay_f*F0));
        beq(end)  = beq(end) + Ht_Chng;
    end
    
    
    time = [time time(end)+dt];
    tip_hist = [tip_hist x(2)];
    
    if(contact)
        tip_force_hist = [tip_force_hist lambda.eqlin(4)];
    else
        tip_force_hist = [tip_force_hist 0];
    end
    
    C_psi_hist = [C_psi_hist C_psi];
    [phi_s] = interpolate_func_2(Invariant_sample,C_phi);
    phi_s_hist = [phi_s_hist phi_s];
    [C_s,C1_s,C2_s] = interpolate(Invariant_sample,x);
    C_s_hist = [C_s_hist C_s];
    V_hist = [V_hist V];



    area_hist = [area_hist Area_functional(Aux)];
    max_psi_hist  = [max_psi_hist  max(psi)];
    max_phi_hist  = [max_phi_hist  max(phi_s)/Material.phi_bar];

    
    
    if (Constraints.iSurf_loc ==1) %&& count==0
        if(Constraints.iVol_s==1)
            tension = interpolate_func(Invariant,lambda.eqnonlin(2:end),Invariant.lower);
        else
            tension = interpolate_func(Invariant,lambda.eqnonlin,Invariant.lower);
        end
       
    else
        tension = lambda.eqnonlin*ones(size(Invariant.gaussp));
        %count=1;
    end
    sigma = [sigma tension];
    
    %plot respect to archlength
    D= (Aux.au.* Invariant.gweight);
    S=zeros(size(D));
    for i=1:size(D)
        S(i)= sum(D(1:i));
    end
    %%%%%%%%%%%%%%%%%%%%%
    figure(3)
    hold off
    plot(C_s( :,1),C_s( :,2))
    hold on
    plot(x(1:2:end),x(2:2:end),'ko')
    axis equal
   
   
     if(v0_dendrite*(1-exp(-decay_f*F0))*time(end)>x(2) && contact==0)
        dt          = dt/5; 
        Aeq         = [Aeq; zeros(1,length(Aeq))];
        Aeq(end,2)    = 1;
        Ht_Chng     = dt*v0_dendrite*(1-exp(-decay_f*F0));
        beq(end+1)      = x(2) + Ht_Chng;
        contact     = 1;
        disp('We establish contact between growing dendritic network and membrane')
    end
    
    
    Proj_area_i = pi*x(end-1)^2;
    %Areal_Strain = [Areal_Strain 100.*(Proj_area_i - Proj_area_0)./Proj_area_0];
    Areal_Strain = [Areal_Strain (pi*Material.R0^2 - Proj_area_i )/(pi*Material.R0^2 - Proj_area_0 )];
    figure (7)
    title('Normalized excess area vs Time')
    plot(time,Areal_Strain,'--ro','Markersize',4,'LineWidth',2)
    xlabel('Time','FontSize',20,'FontWeight','bold','Color','k')
    ylabel('Normalized excess area','FontSize',20,'FontWeight','bold','Color','k')
    pause(0.1)
    
end
%C_phi_hist = [C_phi_hist_1 C_phi_hist];

figure(990)
plot(time(2:end),max_phi_hist)
xlabel('time')
ylabel('maxphi')
figure(991)
plot(time(2:end),max_psi_hist)
xlabel('time')
ylabel('maxpsi')
figure(992)
plot(time(2:end),max_V_c_hist)
xlabel('time')
ylabel('maxVc')


return

%%
osave = 1;
output_file = 'Examples_Dynamic/OM_Support.mat';

if osave == 1
    save(output_file);
end


%%


















Proj_area_hist = pi.*x_hist(end-1,:).^2;
Areal_Strain = [Areal_Strain 100.*(Proj_area_hist - Proj_area_hist(1))./Proj_area_hist(1)];
figure (7)
title('Areal Strain vs Time')
time = time*1.4;
plot(time,Areal_Strain,'--ro','Markersize',4,'LineWidth',2)
xlabel('Time (in ms)','FontSize',20,'FontWeight','bold','Color','k')
ylabel('Areal Strain','FontSize',20,'FontWeight','bold','Color','k')
set (gca,'fontsize',20);
ofile = strcat('tmp/','A_Strain_Drag_3_Flatten','.png');
print(ofile,'-dpng','-r300')
