function [ip1,ip2,status] = NE_int_line_circle(P,Q,circle)
% NE_int_line_circle - find intersection points of line and circle
% On input:
%     P (1x2 vector): point 1 defining line
%     Q (1x2 vector): point 2 defining line
%     circle (1x5 vector): center, radius and circle pt
% On output:
%     ip1 (1x2 vector): intersection point 1
%     ip2 (1x2 vector): intersection point 2
%     status (int): 0 if no intersection; 1 if single point; 2 if 2 points
% Call:
%     [p1,p2,s] = NE_int_line_circle([-0.5,0.6],[0.5,0.6],[0,0,1,1,0]);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.001;

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

[d,closest_pt] = NE_dist_pt_line([circle(1:2),0],[P,0; Q,0]);
if d>circle(3)  % no intersection
    return
elseif abs(d-circle(3))<ZERO_THRESH
    status = 1;
    ip1 = closest_pt;
    ip2 = closest_pt;
    return
end

[p,s] = CV_total_LS([P(1),Q(1)],[P(2),Q(2)]);
cx = circle(1);
cy = circle(2);
r = circle(3);
a = p(1);
b = p(2);
c = p(3);

if abs(a)<ZERO_THRESH  % horizontal line
    a = 1;
    b = -2*circle(1);
    c = circle(1)^2 + circle(2)^2 - circle(3)^2;
    ip1 = [(-b-sqrt(b^2-4*a*c))/(2*a),0];
    ip2 = [(-b+sqrt(b^2-4*a*c))/(2*a),0];
    if norm(ip1)>=1
        tmp = ip1;
        ip1 = ip2;
        ip2 = ip1;
    end
    if norm(ip2-ip1)<ZERO_THRESH
        status = 1;
    else
        status = 2;
    end
    return
end

if abs(b)<ZERO_THRESH  % vertical line
    aa = 1;
    bb = -2*circle(2);
    cc = circle(2)^2 + circle(1)^2 - circle(3)^2;
    y1 = (-bb-sqrt(bb^2-4*aa*cc))/(2*aa);
    y2 = (-bb+sqrt(bb^2-4*aa*cc))/(2*aa);
    status = 2;
    if y1<1
        ip1 = [0,y1];
        ip2 = [0,y2];
    else
        ip1 = [0,y2];
        ip2 = [0,y1];
    end
    dirc = ip1 - circle(1:2);
    dirc = dirc/norm(dirc);
    dirc = [-dirc(2),dirc(1)];
    theta = posori(acos(dot([0,1],dirc)));
    return
end

X = closest_pt(1:2);
dir = X - circle(1:2);
dir = dir/norm(dir);
dir = [dir(2),-dir(1)];
b = norm(X-circle(1:2));
a = sqrt(circle(3)^2-b^2);
ip1 = X + a*dir;
ip2 = X - a*dir;
if norm(ip1)>=1
    tmp = ip1;
    ip1 = ip2;
    ip2 = tmp;
end
status = 2;
