function [Axn, Ayn, Bxn, Byn] =  BC_CONDITIONS(Ax, Ay, Bx, By, R_out, R_in)

Axn = Ax; Ayn = Ay; Bxn = Bx; Byn = By;

R_tmp = sqrt(Ax.^2 + Ay.^2);
Axn(R_tmp > R_out) = Axn(R_tmp > R_out).*R_out./R_tmp(R_tmp > R_out);
Ayn(R_tmp > R_out) = Ayn(R_tmp > R_out).*R_out./R_tmp(R_tmp > R_out);

R_tmp = sqrt(Bx.^2 + By.^2);
Bxn(R_tmp > R_in) = Bxn(R_tmp > R_in).*R_in./R_tmp(R_tmp > R_in);
Byn(R_tmp > R_in) = Byn(R_tmp > R_in).*R_in./R_tmp(R_tmp > R_in);

%{
D_r = [Bx(2:(end)) - Ax(1:(end-1)); By(2:(end)) - Ay(1:(end-1))];
M_D_r = sqrt(sum(D_r.^2, 1));
D_l = [Bx(1:(end-1)) - Ax(2:end); By(1:(end-1)) - Ay(2:end)];
M_D_l = sqrt(sum(D_l.^2, 1));

L = [Bx - Ax; By - Ay];

L_par_r = D_r./repmat(M_D_r.^2.*sum(D_r.*L(:, 1:(end-1)), 1), 2, 1);
L_perp_r = L(:, 1:(end-1)) - L_par_r;
L_par_rx = Ax(1:(end-1)) + L_par_r(1, :); L_par_ry = Ay(1:(end-1)) + L_par_r(2, :);

SIN_r = D_r(1, :).*L_perp_r(2, :) - D_r(2, :).*L_perp_r(1, :);

Bx_tmp = Bxn(1:(end-1)); By_tmp = Byn(1:(end-1));
Bx_tmp(SIN_r < 0) = L_par_rx(SIN_r < 0); By_tmp(SIN_r < 0) = L_par_ry(SIN_r < 0);
Bxn(1:(end-1)) = Bx_tmp; Byn(1:(end-1)) = By_tmp;

L_par_l = D_l./repmat(M_D_l.^2.*sum(D_l.*L(:, 2:end), 1), 2, 1);
L_perp_l = L(:, 2:end) - L_par_l;
L_par_lx = Ax(2:end) + L_perp_l(1, :); L_par_ly = Ay(2:end) + L_perp_l(2, :);

SIN_l = D_l(1, :).*L_perp_l(2, :) - D_l(2, :).*L_perp_l(1, :);

Bx_tmp = Bxn(2:end); By_tmp = Byn(2:end);
Bx_tmp(SIN_l > 0) = L_par_lx(SIN_l > 0); By_tmp(SIN_l > 0) = L_par_ly(SIN_l > 0);
Bxn(2:end) = Bx_tmp; Byn(2:end) = By_tmp;
%}


end