function [Ax_n, Ay_n, Bx_n, By_n] = E_MIN_EVOLVE_2D(Ax, Ay, Bx, By, KV, Vo, KA, Ao, KB, Bo, KL, Lo, KY, Yo, N_steps, R_out, R_in)

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

%dx = .001;
beta = .01;
c = .9;
Gbx_yolk = zeros(size(Ax_n)); Gby_yolk = zeros(size(Ay_n));
n_yolk_calc = 25;
for n_steps=1:N_steps,
   
    % Calculate the gradient
    
    %[Gax, Gay, Gbx, Gby] = CALC_GRADIENT_NUM_2D(Ax_n, Ay_n, Bx_n, By_n, KV, Vo, KA, Ao, KL, Lo, .001);
    [Gax, Gay, Gbx, Gby] = CALC_GRADIENT_ANYL_2D(Ax_n, Ay_n, Bx_n, By_n, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    if mod(n_steps, n_yolk_calc) == 0, [Gbx_yolk, Gby_yolk] = CALC_GRADIENT_YOLK(Bx_n, By_n, KY, Yo);end
    
    Gbx = Gbx + Gbx_yolk; Gby = Gby + Gby_yolk;
    %figure, plot(Gax, 'b'); hold on; plot(Gay, 'r');title('Ga');
    %figure, plot(Gbx, 'b'); hold on; plot(Gby, 'r');title('Gb');
    Gax(1) = 0; Gbx(1) = 0; Gax(end) = 0; Gbx(end) = 0;
    %Gay(1) = 0; Gby(1) = 0;  Gay(end) = 0;  Gby(end) = 0;
    % Move Until energy minimum
    %EMIN_FUN = @(p) CALC_ENERGY_2D(Ax_n + p*Gax, Ay_n + p*Gay, Bx_n + p*Gbx, By_n + p*Gby, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    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, KY, Yo, R_out, R_in);
    %EMIN_FUN = @(p) CALC_ENERGY_2D(Ax_n - p*Gax, Ay_n - p*Gay, Bx_n - p*Gbx, By_n - p*Gby, KV, Vo, KA, Ao, KB, Bo, KL, Lo);
    %beta = fminsearch(EMIN_FUN, 0);
    %beta = .001;
    
    %%{
    E_OLD = EMIN_FUN(0);
    %Scale the gradient to size dx
    
    
    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;
    Ay_n(Ay_n > 0) = 0;
    Bx_n = Bx_n + beta*Gbx; By_n = By_n + beta*Gby;
    By_n(By_n > 0) = 0;
    %}
    
    
    Ax_n = Ax_n + beta*Gax; Ay_n = Ay_n + beta*Gay;
    Bx_n = Bx_n + beta*Gbx; By_n = By_n + beta*Gby;
    
   
    R_tmp = sqrt(Ax_n.^2 + Ay_n.^2);
    Ax_n(R_tmp > R_out) = Ax_n(R_tmp > R_out).*R_out./R_tmp(R_tmp > R_out);
    Ay_n(R_tmp > R_out) = Ay_n(R_tmp > R_out).*R_out./R_tmp(R_tmp > R_out);
    %{
    Ri = sqrt(Ax_n.^2 + Ay_n.^2);
    Ax_n(1:25) = R_out./Ri(1:25).*Ax_n(1:25); Ay_n(1:25) = R_out./Ri(1:25).*Ay_n(1:25);
    Ax_n((end-25):end) = R_out./Ri((end-25):end).*Ax_n((end-25):end); Ay_n((end-25):end) = R_out./Ri((end-25):end).*Ay_n((end-25):end);
    %}
    R_tmp = sqrt(Bx_n.^2 + By_n.^2);
    Bx_n(R_tmp > R_in) = Bx_n(R_tmp > R_in).*R_in./R_tmp(R_tmp > R_in);
    By_n(R_tmp > R_in) = By_n(R_tmp > R_in).*R_in./R_tmp(R_tmp > R_in);
    
    Ax_left = Ax_n(1:end/2); Ax_right = Ax_n((end/2+1):end);
    Ax_left(Ax_left < 0) = 0; Ax_right(Ax_right > 0) = 0;
    Ax_n(1:end/2) = Ax_left; Ax_n((end/2+1):end) = Ax_right;
    
     %{
    D_r = [Bx_n(12) - Ax_n(12), By_n(12) - Ay_n(12)];
    L_r = [Bx_n(13) - Ax_n(12), By_n(13) - Ay_n(12)];
    
    L_par_r = D_r./norm(D_r)^2*dot(D_r, L_r);
    L_perp_r = L_r - L_par_r;
    
    SIN_r = D_r(1)*L_r(2) - D_r(2)*L_r(1);
    
    if (SIN_r > 0),
        Bx_n(13) = Ax_n(12) + L_par_r(1); By_n(13) = Ay_n(12) + L_par_r(2);
    end
    
    D_l = [Bx_n(28) - Ax_n(28), By_n(28) - Ay_n(28)];
    L_l = [Bx_n(27) - Ax_n(28), By_n(27) - Ay_n(28)];
    
    L_par_l = D_l./norm(D_l)^2*dot(D_l, L_l);
    L_perp_l = L_l - L_par_l;
    
    SIN_l = D_l(1)*L_l(2) - D_l(2)*L_l(1);
    
    if (SIN_l < 0),
        Bx_n(27) = Ax_n(28) + L_par_l(1); By_n(27) = Ay_n(28) + L_par_l(2);
    end
    %}
    
    %[Ax_n, Ay_n, Bx_n, By_n] = BC_CONDITIONS(Ax_n, Ay_n, Bx_n, By_n, R_out, R_in);
    
    % if mod(n_steps, round(N_steps/10)) == 0, disp([num2str(round(n_steps/N_steps*100)), '% done and E_new = ', num2str(CALC_ENERGY_2D(Ax_n, Ay_n, Bx_n, By_n, KV, Vo, KA, Ao, KB, Bo, KL, Lo))]);end
    
   if mod(n_steps, 100) == 0, disp(['E_new = ', num2str(EMIN_FUN(beta)), ' beta = ', num2str(beta), ' n = ', num2str(n_steps)]);end
    
    
    if beta < .00001, break; end
end
disp(['Completed in ', num2str(n_steps)]);
end