function [ip1,ip2,status] = NE_int_E2_2circles(circle1,circle2)
% NE_int_E2_2circles - find R^2 intersection points of 2 circles
% On input:
%     circle1 (1x3 vector): center and radius of circle 1
%     circle2 (1x3 vector): center and radius of circle 2
% On output:
%     ip1 (1x2 vector): intersection point 1
%     ip2 (1x2 vector): intersection point 2
%     status (int): 0 if no intersection
%                   1 if 1 point
%                   2 if 2 points;
%                   3 if same circle
% Call:
%     [p1,p2,st] = NE_int_E2_2circles([0,0,2],[1,0,1]);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.001;

ip1 = [];
ip2 = [];
status = 0;

C1 = circle1(1:2);
r1 = circle1(3);
C2 = circle2(1:2);
r2 = circle2(3);
d12 = norm(C2-C1);

% no intersection
if (d12+r1<r2)|(d12+r2<r1)|(r1+r2<d12)
    return
end

% complete intersection
if d12<ZERO_THRESH&abs(r1-r2)<ZERO_THRESH
    status = 3;
    ip1 = circle1(4:5);
    ip2 = circle2(4:5);
    return
end

% 1 pt intersection: just 1 point intersection
if abs(r1+r2-d12)<ZERO_THRESH
    status = 1;
    dir = C2 - C1;
    dir = dir/norm(dir);
    ip1 = C1 + r1*dir;
    ip2 = ip1;
    return
end

% 1 pt intersection: C1 inside circle2
if abs(d12+r1-r2)<ZERO_THRESH
    status = 1;
    dir = C1 - C2;
    dir = dir/norm(dir);
    ip1 = C1 + r1*dir;
    ip2 = ip1;
    return
end

% 1 pt intersection: C2 inside circle1
if abs(d12+r2-r1)<ZERO_THRESH
    status = 1;
    dir = C2 - C1;
    dir = dir/norm(dir);
    ip1 = C2 + r2*dir;
    ip2 = ip1;
    return
end

% 2 pt intersection
status = 2;
C1p = [0,0];
C2p = C2 - C1;
theta = posori(atan2(C2p(2),C2p(1)));
R = [cos(theta),sin(theta); -sin(theta),cos(theta)];
C2r = R*C2p';
x = (r1^2-r2^2+C2r(1)^2)/(2*C2r(1));
y = sqrt(r1^2-x^2);
ip1 = (R'*[x;y]+C1')';
ip2 = (R'*[x;-y]+C1')';

tch = 0;
