function [C,z1,z4] = NE_ptsPD2line(P,Q)
% NE_ptsPD2line - find geodesic through 2 Poincare disk points
% On input:
%     P (1x2 vector): point 1
%     Q (1x2 vector): point2
% On output:
%     C (1x5 vector) center, radius & pt on circle
%     z1 (1x2 vector): point on unit disk boundary closer to P
%     z4 (1x2 vector): point on unit disk boundary closer to Q
%       i.e., order is: z1 P, Q, z4  for cross ratio calculation
% Call:
%     [C,r,z1,z4] = NE_ptsPD2line(P,Q);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.001;
STEP = 0.01;

C = [];
r = 0;
z1 = [];
z4 = [];

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

if d<ZERO_THRESH&norm(P)<ZERO_THRESH
    C = Q;
    r = Inf;
    [z1,z4,st] = NE_int_line_circle([0,0],Q,[0,0,1,1,0]);
    d1P = NE_dist_PD_pt_pt(P,z1);
    d4P = NE_dist_PD_pt_pt(P,z4);
    if d4P<d1P
        tmp = z1;
        z1 = z4;
        z4 = tmp;
    end
    C = [C,r,0,0];
    return
elseif d<ZERO_THRESH&norm(Q)<ZERO_THRESH
    C = P;
    r = Inf;
    [z1,z4,st] = NE_int_line_circle([0,0],P,[0,0,1,1,0]);
    d1P = NE_dist_PD_pt_pt(P,z1);
    d4P = NE_dist_PD_pt_pt(P,z4);
    if d4P<d1P
        tmp = z1;
        z1 = z4;
        z4 = tmp;
    end
    C = [C,r,0,0];
    return
elseif d<ZERO_THRESH
    C = P;
    r = Inf;
    [z1,z4,st] = NE_int_line_circle(P,Q,[0,0,1,1,0]);
    d1P = NE_dist_PD_pt_pt(P,z1);
    d4P = NE_dist_PD_pt_pt(P,z4);
    if d4P<d1P
        tmp = z1;
        z1 = z4;
        z4 = tmp;
    end
    C = [C,r,P];
    return
end

unit_circle = [0,0,1,1,0];
circle = NE_orthog_circle(P,Q);
C = circle(1:2);
r = circle(3);
[z1,z4,st] = NE_int_2circles(unit_circle,[C,r,P]);
z1 = z1';
z4 = z4';
d1P = norm(P-z1);
d4P = norm(P-z4);
if d4P<d1P
    tmp = z1;
    z1 = z4;
    z4 = tmp;
end
C = [C,r,P];
