function theta = NE_angle_2circles(C1,C2)
% NE_angle_2circles - angle between tangents of two PD circles
% On input:
%     C1 (1x5 vector): center, radius, pt on circle 1
%     C2 (1x5 vector): center, radius, pt on circle 2
% On output:
%     theta (float): angle between tangents of circle intersection
% Call:
%     theta = NE_angle_2circles(C1,C2);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

theta = [];

if isinf(C1(3))&isinf(C2(3))   % both straight
    dir1 = [C1(4)-C1(1),C1(5)-C1(2)];
    dir1 = dir1/norm(dir1);
    dir2 = [C2(4)-C2(1),C2(5)-C2(2)];
    dir2 = dir2/norm(dir2);
    theta = posori(acos(dot(dir1,dir2)));
elseif isinf(C1(3))  % C1 straight; C2 curved
    [p1,p2,s] = NE_int_line_circle(C1(1:2),C1(4:5),C2);
    if s==0  % no intersection
        return
    end
    if s==1    % straight line is tangent to circle
        theta = 0;
        return
    end
    if norm(p1)<1  % get intersection point inside PD
        p = p1;
    else
        p = p2;
    end
    dir1 = [C1(4)-C1(1),C1(5)-C1(2)];
    dir1 = dir1/norm(dir1);
    cdir = p - C2(1:2);
    cdir = cdir/norm(cdir);
%    cdir = [cdir(2),-cdir(1)];
    theta = posori(acos(dot(dir1,cdir)));
elseif isinf(C2(3))  % C2 straight; C1 curved
    [p1,p2,s] = NE_int_line_circle(C2(1:2),C2(4:5),C1);
    if s==0   % no intersection
        return
    end
    if s==1    % straight line is tangent to circle
        theta = pi/2;
        return
    end
    if norm(p1)<1  % get intersection point inside PD
        p = p1;
    else
        p = p2;
    end
    dir1 = [C2(4)-C2(1),C2(5)-C2(2)];
    dir1 = dir1/norm(dir1);
    cdir = p - C1(1:2);
    cdir = cdir/norm(cdir);
    theta = posori(acos(dot(dir1,cdir)));
else  % both curved
    [ip1,ip2,s] = NE_int_2circles(C1,C2);
    if s==0
        return
    end
    if norm(ip1)<1
        p = ip1;
    else
        p = ip2;
    end
    dir1 = p - C1(1:2);
    dir1 = dir1/norm(dir1);
    dir2 = p - C2(1:2);
    dir2 = dir2/norm(dir2);
    theta = posori(acos(dot(dir1,dir2)));
end
if theta>pi/2
    theta = theta - pi/2;
end
tch = 0;
