close all; clear all; clc;
load('V_cell'); N = size(V_cell{1}, 2)/6; L = length(V_cell);

% Declare the initial apical and basal positions
%for j=1:L,
    j=4;
    Ax_o = V_cell{j}(1, 1:N); Ay_o = V_cell{j}(2, 1:N);
    Bx_o = V_cell{j}(1, (5*N+1):(6*N)); By_o = V_cell{j}(2, (5*N+1):(6*N));
    
%end

Ax = Ax_o; Ay = Ay_o; 
Bx = Bx_o; By = By_o;

% Declare initial parameters
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);
%Lo = mean(Lo)*ones(size(Lo));


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;
%Vo = mean(Vo)*ones(size(Vo));

% Declare strengths
M = 0;
KL_val = .5; 
KV_val = .01; 
KB_val = .1;
N_steps = 5000;

Con_size = 7;

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);

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

%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, hold on;title(['Original Frame j = ', num2str(j), ' with initial forces'], 'Fontsize', 20);
plot(Ax_o, Ay_o, 'r-o', 'linewidth', 3);
plot(Bx_o, By_o, 'b-o', 'linewidth', 3);
for i=1:N, plot([Ax_o(i), Bx_o(i)], [Ay_o(i), By_o(i)], 'k-', 'linewidth', 3);end
quiver(Ax_o, Ay_o, Gax, Gay);
quiver(Bx_o, By_o, Gbx, Gby);

%% Evaluate the basal side
close all;
h = figure('Position', [300 100 800 800]);
id_movie = 1;
X_max = max(Ax) + 1; X_min = min(Ax) - 1; Y_max = max(Ay) + 3; Y_min = min(Bx) - 1;

for j=4:26,

% Read the intial     
Ax_o = V_cell{j}(1, 1:N); Ay_o = V_cell{j}(2, 1:N);
Bx_m = V_cell{j}(1, (5*N+1):(6*N)); By_m = V_cell{j}(2, (5*N+1):(6*N));
Ax = Ax_o; Ay = Ay_o; 


% Redeclare the initial parameters
Bo = sqrt((Bx_m(1:(end-1)) - Bx_m(2:end)).^2 + (By_m(1:(end-1)) - By_m(2:end)).^2);
Lo = sqrt((Ax_o - Bx_m).^2 + (Ay_o - By_m).^2);

Ax_n = Ax; Ay_n = Ay;
Bx_n = Bx; By_n = By;

beta = .01; R_out = 14; R_in = 10;
c = .9;

for n = 1:N_steps,
    [Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Ax_n, Ay_n, Bx_n, By_n, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    
    Gax = 0*Gax; Gay = 0*Gay;
    
    EMIN_FUN = @(p) CALC_ENERGY_2D_CURVE(Ax_n + p*Gax, Ay_n + p*Gay, Bx_n + p*Gbx, By_n + p*Gby, KV, Vo, KA, Ao, KB, Bo, KL, Lo, 0, 0, R_out, R_in);

    E_OLD = EMIN_FUN(0);
    
    while (1),
        if EMIN_FUN(beta) >= E_OLD && (beta > .00001),
            beta = c*beta;
        else
            break;
        end
    end
        
    
    while (1),  
            E_NEW = EMIN_FUN(beta);
            
            if E_NEW < E_OLD,
                beta = beta + .01;
                E_OLD = E_NEW;
            else
                break;
            end
    end
    
    Ax_n = Ax_n + beta*Gax; Ay_n = Ay_n + beta*Gay;
    Bx_n = Bx_n + beta*Gbx; By_n = By_n + beta*Gby;
    
     disp(['E_new = ', num2str(EMIN_FUN(beta)), ' beta = ', num2str(beta), ' n = ', num2str(n)]);
     if beta < .00001, break; end
    
end
Bx = Bx_n; By = By_n;

[Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Ax_n, Ay_n, Bx_n, By_n, KV, Vo, KA, Ao, KB, Bo, KL, Lo);

%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');

hold off;plot(Ax_n, Ay_n, 'r-o', 'linewidth', 3);hold on;
plot(Bx_n, By_n, 'b-o', 'linewidth', 3);plot(Bx_m, By_m, 'k--', 'linewidth', 3);
for i=1:N, plot([Ax_n(i), Bx_n(i)], [Ay_n(i), By_n(i)], 'k-', 'linewidth', 3);end
for i=1:N, plot([Ax_n(i), Bx_m(i)], [Ay_n(i), By_m(i)], 'k--', 'linewidth', 1);end

quiver(Bx_n, By_n, Gbx, Gby);
title(['Final Frame j = ', num2str(j), ' with forces KL = ', num2str(KL_val)], 'Fontsize', 20);
axis([X_min X_max Y_min Y_max]);
drawnow; pause(.1);
M_tmp(id_movie) = getframe(h);id_movie = id_movie+1;pause(.1);
end
movie2avi(M_tmp, 'Movie_temp');

