function circle = NE_orthog_circle(P,Q)
% NE_orthog_circle - find orthogonal circle to unit circle through P & Q
% On input:
%     P (1x2 vector): point inside unit disk
%     Q (1x2 vector): point inside unit disck distinct from P
% On output:
%     circle (1x5 vector): center, radius & circle pt of orthogonal circle
% Call:
%     circle = NE_orthog_circle([0.2,0.8],[0.4,0.6]);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.0001;

line = [P,0;Q,0];
d = cv_dist_pt_line([0,0,0],line);

if d<ZERO_THRESH    % straight lines
    if norm(P)<ZERO_THRESH&norm(Q)<ZERO_THRESH
        circle = [0,0,0,0,0];
        return
    elseif norm(P)<ZERO_THRESH
        circle = [Q,Inf,0,0];
        return
    else
        circle = [P,Inf,0,0];
        return
    end
end

A = [P;Q];
b = [(1+norm(P)^2)/2;(1+norm(Q)^2)/2];
C = (A\b)';
r2 = norm(P-C);
circle = [C,r2,P];
