N = 20000;
M = str2num(get(handles.A_mag_value, 'String'));
S = str2num(get(handles.A_sig_value, 'String'));

W = N/5;
W_mag = str2num(get(handles.W_mag_value, 'String'));
M_b = str2num(get(handles.B_mag_value, 'String'));
X_cut = str2num(get(handles.B_radius_value, 'String'));

R_out = 14*ones(1, N);
R_in = 9*ones(1, N);
X = linspace(-5, 5, N);
Theta = linspace(-pi/2, pi/2, N);
B_id = R_in.*cos(Theta) > X_cut;
id = min(find(B_id > 0));
id2 = max(find(B_id > 0));





%R_in(id:(id+1250)) = R_in(id:(id+1250)) + M*linspace(0, 5, 1251);
%R_in((id2-1250):id2) = R_in((id2-1250):id2) + M*linspace(5, 0, 1251);
% NOTE THAT id < id2
R_in((id-W):id) = R_in((id-W):id) + W_mag*linspace(0, 5, W+1); % Left ear
R_in(id2:(id2+W)) = R_in(id2:(id2+W)) + W_mag*linspace(5, 0, W+1); % Rigth ear
%R_in((id+1):(id2-1)) = R_in((id+1):(id2-1)) + abs(W_mag*linspace(-sqrt(5), sqrt(5), id2-id-1).^2);
K = 6; R_in((id+1):(id2-1)) = R_in((id+1):(id2-1)) + abs(W_mag*linspace(-(5)^(1/K), (5)^(1/K), id2-id-1).^K);
X_in = R_in.*sin(Theta); Y_in = R_in.*cos(Theta);




Y_tmp = 1/2*(X).^2 + (1/8*(X).^4 + 1/16*(X).^6 + 1/32*(X).^8)*.5 + abs(X);


Y_tmp = -M_b*Y_tmp;%./max(Y_tmp);
Y_tmp = Y_tmp - Y_tmp(id+1);

Y_in((id+W):(id2-W)) = Y_in(id);
Y_in = Y_in - Y_tmp.*B_id;
%R_in(B_id) = R_in(id);
%R_in = R_in - Y_tmp.*B_id;





L_tmp = length(handles.V_cell{handles.V_cell_id}(1, :));
V_RAW_X = handles.V_cell{handles.V_cell_id}(1, :); V_X_tmp = reshape(V_RAW_X, L_tmp/6, 6);
V_RAW_Y = handles.V_cell{handles.V_cell_id}(2, :); V_Y_tmp = reshape(V_RAW_Y, L_tmp/6, 6);

V_RAW_X = V_X_tmp*cos(handles.rotate) + V_Y_tmp*sin(handles.rotate);
V_RAW_Y = -V_X_tmp*sin(handles.rotate) + V_Y_tmp*cos(handles.rotate);
save('TEST', 'V_RAW_X', 'V_RAW_Y');
%[M, S] = Apical_fit(V_RAW_X(:, 1), V_RAW_Y(:, 1));

R_out = R_out - M*exp(-(X/S).^2);
X_out = R_out.*sin(Theta); Y_out = R_out.*cos(Theta);

Apical_Dist = handles.Apical_Dist;
Basal_Dist = ones(length(Apical_Dist), 1)*.8*8/12;
%V_apical = Populate(X_out, Y_out, Apical_Dist);
V_apical =  [V_RAW_X(:, 1), V_RAW_Y(:, 1)];
V_basal =  [V_RAW_X(:, 6), V_RAW_Y(:, 6)];
%V_basal = Populate(X_in, Y_in, Basal_Dist);
Volumes = zeros(size(V_apical, 1)-1, 1);

%V_basal = Calc_Basal_Constant_Volume(V_apical, X_in, Y_in, 6.65);

set(handles.msg_disp, 'String', 'Calculating...');pause(.1);
%V_basal = Calc_Basal_E_MIN(V_apical, V_basal, 6.65);
set(handles.msg_disp, 'String', 'DONE!!!');
V_grad = Calc_Gradient(V_apical, V_basal, 1);


for i=1:(size(V_apical, 1)-1),
   A1 = V_apical(i, :) - V_basal(i, :);
   A2 = V_apical(i, :) - V_basal(i+1, :);
   A3 = V_apical(i, :) - V_apical(i+1, :);
    
   Volumes(i) = abs(A1(1)*A2(2) - A1(2)*A2(1)) + ...
       abs(A2(1)*A3(2) - A2(2)*A3(1));
   
end

plot(handles.axes3, Volumes, 'r-*', 'linewidth', 3);




set(handles.axes1, 'NextPlot', 'replace');
plot(handles.axes1, V_RAW_X(:, 1), V_RAW_Y(:, 1), 'm-*', 'linewidth', 3);
set(handles.axes1, 'NextPlot', 'add');
plot(handles.axes1, X_out, Y_out, 'b', 'linewidth', 3);
plot(handles.axes1, V_RAW_X(:, 6), V_RAW_Y(:, 6), 'm-*', 'linewidth', 3);
plot(handles.axes1, X_in, Y_in, 'g', 'linewidth', 3);
plot(handles.axes1, V_apical(:, 1), V_apical(:, 2), 'r-*');
plot(handles.axes1, V_basal(:, 1), V_basal(:, 2), 'r-*');

for i = 1:size(V_apical, 1),
  plot(handles.axes1, [V_apical(i, 1) V_basal(i, 1)], [V_apical(i, 2) V_basal(i, 2)], 'k', 'linewidth', 3); 
end

quiver(handles.axes1, V_basal(:, 1), V_basal(:, 2), V_grad(:, 1), V_grad(:, 2));

%axis equal;