function [Gbx_y, Gby_y] = CALC_GRADIENT_YOLK(Bx, By, KY, Yo)

S = length(Bx);
Gbx_y = zeros(size(Bx));
Gby_y = zeros(size(By));
dx = .0001;

Bx_right = Bx; Bx_left = Bx;
By_up = By; By_down = By;

for i = 20:(S-19),
   Bx_right(i) = Bx_right(i) + dx; Bx_left(i) = Bx_left(i) - dx;
   By_up(i) = By_up(i) + dx; By_down(i) = By_down(i) - dx;
   
   %E_right = KY*(Yo - polyarea(Bx_right, By))^2;
   %E_left = KY*(Yo - polyarea(Bx_left, By))^2;
   
   %E_up = KY*(Yo - polyarea(Bx, By_up))^2;
   %E_down = KY*(Yo - polyarea(Bx, By_down))^2;
   
   E_right = KY*(Yo - f_closed_area([Bx_right', By']))^2;
   E_left = KY*(Yo - f_closed_area([Bx_left', By']))^2;
   
   E_up = KY*(Yo - f_closed_area([Bx', By_up']))^2;
   E_down = KY*(Yo - f_closed_area([Bx', By_down']))^2;
   
   Gbx_y(i) = -(E_right - E_left)/(2*dx);
   Gby_y(i) = -(E_up - E_down)/(2*dx);
   
   Bx_right(i) = Bx(i); Bx_left(i) = Bx(i);
   By_up(i) = By(i); By_down(i) = By(i);
end



end