% THIS IS THE 2D SIMULATION
close all; clear all; clc;
N = 80; N_mid = N/2;
Theta = linspace(-pi/2, 3*pi/2, N);
R_out = 14; R_in = 10;

x_space = 10; y_space = 40;
L_size = N*x_space;

% Init Position
Ax = R_out*cos(Theta); Ay = R_out*sin(Theta);
Bx = R_in*cos(Theta); By = R_in*sin(Theta);
Axn = Ax; Ayn = Ay; Bxn = Bx; Byn = By; 
%{
Ax = linspace(0, L_size, N);
Bx = Ax;
Ay = zeros(1, N); By = Ay - y_space;
%}

% Init equilibrium points
Ao = sqrt((Ax(1:(end-1)) - Ax(2:end)).^2 + (Ay(1:(end-1)) - Ay(2:end)).^2); Ai = Ao;
Bo = sqrt((Bx(1:(end-1)) - Bx(2:end)).^2 + (By(1:(end-1)) - By(2:end)).^2);
Lo = sqrt((Ax - Bx).^2 + (Ay - By).^2);
Yo = polyarea(Bx, By);

V1x = Ax(1:(end-1)) - Ax(2:end); V1y = Ay(1:(end-1)) - Ay(2:end);
V2x = Ax(1:(end-1)) - Bx(2:end);V2y = Ay(1:(end-1)) - By(2:end); 
V3x = Ax(1:(end-1)) - Bx(1:(end-1));V3y = Ay(1:(end-1)) - By(1:(end-1));

Vo = abs(V1x.*V2y - V1y.*V2x)/2 + abs(V2x.*V3y - V2y.*V3x)/2;
% Init strengths
M = .6; S = 15; D = R_out^2*pi/4; M_o = .05;
KL_val = .005; 
KV_val = 0.01;
KB_val = .01; 
KY = .000001;
N_steps = 5000;

Con_size = 7;

X_tmp = linspace(-D, D, N-1);
KA = M*exp(-(X_tmp/S).^2) + M_o;% KA = M*ones(1, N-1);
KL = KL_val*ones(1, N); KV = KV_val*ones(1, N-1); KB = KB_val*ones(1, N-1);
%KB((N_mid-Con_size):(N_mid+Con_size)) = 0;

%close all;

Con_factor = 1.0;
Exp_factor = 1.0;


Ao((N_mid-Con_size):(N_mid+Con_size)) = Con_factor*Ao(1);
%Bo(2:(N_mid-Con_size-1)) = .8*Bo(1);%Bo((N_mid-Con_size-3):(N_mid-Con_size-1)) = Con_factor*Bo(1);
%Bo((N_mid+Con_size+1):end) = .8*Bo(1);
Bo = Bo;

%KB((N_mid-Con_size):(N_mid+Con_size)) = 0;

%[Gax, Gay, Gbx, Gby] = CALC_GRADIENT_NUM_2D(Ax, Ay, Bx, By, KV, Vo, KA, Ao, KL, Lo, .001);
[Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Ax, Ay, Bx, By, KV, Vo, KA, Ao, KB, Bo, KL, Lo);

%[Gax_yolk, Gay_yolk] = CALC_GRADIENT_YOLK_APICAL(Ax, Ay, KY, R_out);
[Gbx_yolk, Gby_yolk] = CALC_GRADIENT_YOLK(Bx, By, KY, Yo);

%X_max = max(Ax) + 1; X_min = min(Ax) - 1; Y_max = max(Ay) + 3; Y_min = min(By) - 5;
X_max = max(Ax) + 1; X_min = min(Ax) - 1; Y_max = X_max; Y_min = X_min;

figure, plot(Gax, 'r');hold on;plot(Gbx, 'b');title('Force along X');legend('Apical', 'Basal');
figure, plot(Gay, 'r');hold on;plot(Gby, 'b');title('Force along Y');legend('Apical', 'Basal');
figure('Position', [50 500 500 400]), hold on;
plot(Ax, Ay, 'r-*', 'linewidth', 2); plot(Bx, By, 'b-*', 'linewidth', 2);
for i=1:N, plot([Ax(i) Bx(i)], [Ay(i) By(i)], 'k', 'linewidth', 2);end;axis equal;
quiver(Ax, Ay, Gax, Gay);
quiver(Bx, By, Gbx, Gby);
title('Initial Position');
%%
tic,
%Ao = Constrict_Apical(Ao, Con_size, Con_factor); figure, plot(KA); hold on; plot(Ao/max(Ao)*max(KA));
E_old = CALC_ENERGY_2D(Ax, Ay, Bx, By, KV, Vo, KA, Ao, KB, Bo, KL, Lo);disp(['E_old = ', num2str(E_old)]);
close all;
h = figure('Position', [300 30 800 800]);
Axn = Ax; Ayn = Ay; Bxn = Bx; Byn = By; 
id_movie=1;
Con_change = .02; Exp_change = .025;
while (1),
    Con_factor = Con_factor - Con_change; Con_factor = round(Con_factor*1000)/1000;
    
    if Con_factor < -Con_change/2, break;end; 
    
    if Con_factor == .1, Con_change = Con_change/4; Exp_change = Exp_change/4;end
    
    
    if (Con_factor <= .8 && Con_factor > .2), 
        Exp_factor = Exp_factor + Exp_change;
    elseif (Con_factor <= .2),
        Exp_factor = Exp_factor - 3*Exp_change;
    end
    
    if Exp_factor < 1.0, Exp_factor = 1.0;end
    
    Ao((N_mid-Con_size):(N_mid+Con_size)) = Con_factor*Ao(1);
    
    Lo((N_mid-Con_size+1):(N_mid+Con_size)) = Exp_factor*Lo(1);
    Lo(N_mid-Con_size) = (Lo(1) + Lo(end/2))/2;
    Lo(N_mid+Con_size+1) = (Lo(1) + Lo(end/2))/2;
    
    Bo((N_mid-Con_size-3):(N_mid-Con_size-1)) = Con_factor*Bo(1);
    Bo((N_mid+Con_size+1):(N_mid+Con_size+3)) = Con_factor*Bo(1);
    %[Gbx, Gby] = CALC_GRADIENT_FULL_2D(Ax, Ay, Bx, By, KV, Vo, KB, Bo, KL, Lo);
    %[Gax, Gay] = CALC_GRADIENT_FULL_2D(Bx, By, Ax, Ay, KV, Vo, KA, Ao, KL, Lo);
    [Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Axn, Ayn, Bxn, Byn, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    [Gbx_yolk, Gby_yolk] = CALC_GRADIENT_YOLK(Bxn, Byn, KY, Yo);
    
    [Axn, Ayn, Bxn, Byn] = E_MIN_EVOLVE_2D(Axn, Ayn, Bxn, Byn, KV, Vo, KA, Ao, KB, Bo, KL, Lo, KY, Yo, N_steps, R_out, R_in);
    
    
    
    hold off; plot(Axn, Ayn, 'r-*', 'linewidth', 2); hold on;
    plot(Bxn, Byn, 'b-*', 'linewidth', 2); quiver(Axn, Ayn, Gax, Gay); quiver(Bxn, Byn, Gbx, Gby); quiver(Bxn, Byn, Gbx_yolk, Gby_yolk);
    for i=1:N, plot([Axn(i) Bxn(i)], [Ayn(i) Byn(i)], 'k', 'linewidth', 2);end;
    axis([X_min X_max Y_min Y_max]);%axis equal;
    title(['A_0 = ', num2str(Ao(1)*Con_factor), ' \sigma = ', num2str(S), ' KL = ', num2str(KL_val), ' KB = ', num2str(KB_val), ' KY = ', num2str(KY)], 'FontSize', 20);
    drawnow; pause(.1);
    M_tmp(id_movie) = getframe(h);id_movie = id_movie+1;clc;pause(.1);
end

%movie2avi(M_tmp, 'Movie_temp');
E_new = CALC_ENERGY_2D(Axn, Ayn, Bxn, Byn, KV, Vo, KA, Ao, KB, Bo, KL, Lo); disp(['E_new = ', num2str(E_new)]);
T = toc;

disp(['Elasped time was ', num2str(T), ' seconds']);
movie2avi(M_tmp, 'Movie_temp');
%% STart decrasing
%{
Con_factor = 0;
N_steps = 5000;

Ao((N_mid-Con_size):(N_mid+Con_size)) = Con_factor*Ao(1);
Bo((N_mid-Con_size-3):(N_mid-Con_size-1)) = Con_factor*Bo(1);
Bo((N_mid+Con_size+1):(N_mid+Con_size+3)) = Con_factor*Bo(1);
Basal_Con = 1.0;
while (1),
    
    Basal_Con = Basal_Con - .1; if Basal_Con < 0, break;end
    
    KB((N_mid-Con_size):(N_mid+Con_size)) = KB(1)*Basal_Con;
    
    [Axn, Ayn, Bxn, Byn] = E_MIN_EVOLVE_2D(Axn, Ayn, Bxn, Byn, KV, Vo, KA, Ao, KB, Bo, KL, Lo, KY, Yo, N_steps, R_out, R_in);
    
    
    [Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Axn, Ayn, Bxn, Byn, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    
    %figure
    hold off; plot(Axn, Ayn, 'r-*', 'linewidth', 2); hold on;
    plot(Bxn, Byn, 'b-*', 'linewidth', 2); quiver(Axn, Ayn, Gax, Gay); quiver(Bxn, Byn, Gbx, Gby); quiver(Bxn, Byn, Gbx_yolk, Gby_yolk);
    for i=1:N, plot([Axn(i) Bxn(i)], [Ayn(i) Byn(i)], 'k', 'linewidth', 2);end;
    axis([X_min X_max Y_min Y_max]);%axis equal;
    title(['A_0 = ', num2str(Ao(1)*Con_factor), ' \sigma = ', num2str(S), ' KL = ', num2str(KL_val), ' KB = ', num2str(KB_val), ' KY = ', num2str(KY)], 'FontSize', 20);
    drawnow; pause(.1);
    M_tmp(id_movie) = getframe(h);id_movie = id_movie+1;clc;pause(.1);
end
movie2avi(M_tmp, 'Movie_temp');
%% Plot Final
%%{
figure, plot(KA);
figure('Position', [600 500 500 400]), hold on;
plot(Ax, Ay+2*y_space, 'r-*', 'linewidth', 2); plot(Bx, By+2*y_space, 'b-*', 'linewidth', 2);
for i=1:N, plot([Ax(i) Bx(i)], [Ay(i)+2*y_space By(i)+2*y_space], 'k', 'linewidth', 2);end;
plot(Axn, Ayn, 'r-*', 'linewidth', 2); plot(Bxn, Byn, 'b-*', 'linewidth', 2);
for i=1:N, plot([Axn(i) Bxn(i)], [Ayn(i) Byn(i)], 'k', 'linewidth', 2);end;axis equal;
title('final Position');
%%
Af = sqrt((Axn(1:(end-1)) - Axn(2:end)).^2 + (Ayn(1:(end-1)) - Ayn(2:end)).^2);
Bf = sqrt((Bxn(1:(end-1)) - Bxn(2:end)).^2 + (Byn(1:(end-1)) - Byn(2:end)).^2);
Lf = sqrt((Axn - Bxn).^2 + (Ayn - Byn).^2);

V1x = Axn(1:(end-1)) - Axn(2:end); V1y = Ayn(1:(end-1)) - Ayn(2:end);
V2x = Axn(1:(end-1)) - Bxn(2:end);V2y = Ayn(1:(end-1)) - Byn(2:end); 
V3x = Axn(1:(end-1)) - Bxn(1:(end-1));V3y = Ayn(1:(end-1)) - Byn(1:(end-1));

Vf = abs(V1x.*V2y - V1y.*V2x)/2 + abs(V2x.*V3y - V2y.*V3x)/2;
figure('Position', [600 50 500 400]); 
subplot(3,1,1), plot(Ai, 'r', 'linewidth', 2); hold on; title('Apical')
plot(Af, 'r--', 'linewidth', 2);plot(Ao, 'k', 'linewidth', 2);
subplot(3,1,2), plot(Lo, 'b', 'linewidth', 2);hold on; title('Length');
plot(Lf, 'b--', 'linewidth', 2); 
subplot(3,1,3), plot(Vo, 'm', 'linewidth', 2);hold on; title('Volume');
plot(Vf, 'm--', 'linewidth', 2);
%legend('Apical', 'Lengths', 'Volume');% title('Initial and Final Parameters');
%}
%%
%{
close all;
Ri = sqrt(Axn.^2 + Ayn.^2);
Xci = R_out./Ri.*Axn; Yci = R_out./Ri.*Ayn;
L_x = Xci - Axn; L_y = Yci - Ayn;
A_x = Xci(2:end) - Xci(1:(end-1)); A_y = Yci(2:end) - Yci(1:(end-1));
B_x = Axn(2:end) - Axn(1:(end-1)); B_y = Ayn(2:end) - Ayn(1:(end-1));
[Gax_yolk, Gay_yolk] = CALC_GRADIENT_YOLK(Axn, Ayn, KY, R_out);
Yolk = .5*abs(L_x(1:(end-1)).*A_y - L_y(1:(end-1)).*A_x) + .5*abs(L_x(2:end).*B_x - L_y(2:end).*B_y);
figure, plot(Yolk);

figure, plot(Gax_yolk, 'b'), hold on, plot(Gay_yolk, 'g');

figure, plot(Axn, Ayn, 'r-*', 'linewidth', 2); hold on;
    plot(Bxn, Byn, 'b-*', 'linewidth', 2); 
    quiver(Axn, Ayn, Gax_yolk, Gay_yolk);
    for i=1:N, plot([Axn(i) Bxn(i)], [Ayn(i) Byn(i)], 'k', 'linewidth', 2);end;
%}


