function pts = NE_PD_seg_pts(pt1,pt2)
% NE_PD_seg_pts - return set of points along PD segment from pt1 to pt2
% On input:
%     pt1 (1x2 vector): first point
%     pt2 (1x2 vector): second point
% On output:
%     pts (nx2 array): points along segment
% Call:
%     pt = NE_PD_seg_pts([-0.5,.5],[0.5,0.5]);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.001;
DX = 0.01;

pts = [];

if norm(pt2-pt1)<ZERO_THRESH
    pts = pt1;
    return
end

[C,z1,z4] = NE_ptsPD2line(pt1,pt2);
c = C(1:2);
r = C(3);
if r==Inf  % diameter line
    dir = pt2 - pt1;
    dir = dir/norm(dir);
    dist = norm(pt2-pt1);
    num_steps = floor(dist/DX);
    pts = zeros(num_steps,2);
    pts(1,:) = pt1;
    for s = 2:num_steps
        pts(s,:) = pt1 + (s-1)*DX*dir;
    end
else
    theta1 = posori(atan2(pt1(2)-c(2),pt1(1)-c(1)));
    theta2 = posori(atan2(pt2(2)-c(2),pt2(1)-c(1)));
    if theta1<theta2
        alpha = theta1;
        beta = theta2;
    else
        alpha = theta2;
        beta = theta1;
    end
    t = [alpha:0.01:beta];
    x = c(1) + r*cos(t);
    y = c(2) + r*sin(t);
    if abs(theta2-theta1)>pi
        t1 = [beta:0.01:2*pi];
        if t1(end)==2*pi
            t1(end) = [];
        end
        t2 = [0:0.01:alpha];
        x1 = c(1) + r*cos(t1);
        y1 = c(2) + r*sin(t1);
        x2 = c(1) + r*cos(t2);
        y2 = c(2) + r*sin(t2);
        x = [x1,x2];
        y = [y1,y2];
    end
    pts = [x',y'];
end

tch = 0;
