
function plot_ellblend(A,B,C,alpha,beta)

% Function that fills specific circular domains define

A1=A(1,:); B1=B(1,:); C1=C(1,:);
A2=A(2,:); B2=B(2,:); C2=C(2,:);

% P(theta)=A1*cos(theta)+B1*sin(theta)+C1
% Q(theta)=A2*cos(theta)+B2*sin(theta)+C2
% R = {(x,y)=t*P(theta)+(1-t)*Q(theta), t in [0,1], theta in [alpha,beta],
% 0<beta-alpha<=2*pi}

t=linspace(0,1,500); t=t';
theta=linspace(alpha,beta,1000);
Px=(A1(1))*cos(theta)+(B1(1))*sin(theta)+C1(1);
Py=(A1(2))*cos(theta)+(B1(2))*sin(theta)+C1(2);
Qx=(A2(1))*cos(theta)+(B2(1))*sin(theta)+C2(1);
Qy=(A2(2))*cos(theta)+(B2(2))*sin(theta)+C2(2);

X=[]; Y=[];
XB=[]; YB=[];

% points in the domain

for k=1:length(Px)

    PxL=Px(k); QxL=Qx(k);
    xL=t*PxL+(1-t)*QxL;
    X=[X; xL];

    PyL=Py(k); QyL=Qy(k);
    yL=t*PyL+(1-t)*QyL;
    Y=[Y; yL];

    XB=[XB; xL(1); xL(end)]; 
    YB=[YB; yL(1); yL(end)]; 

end


bright_gray=[236, 236, 236]/256;
plot(X,Y,'.','Color',bright_gray);
hold on; 
plot(XB,YB,'k.','LineWidth',2);

if abs(beta-alpha) < 2*pi
    plot([XB(1) XB(2)],[YB(1) YB(2)],'k-','LineWidth',2);
    plot([XB(end-1) XB(end)],[YB(end-1) YB(end)],'k-','LineWidth',2);
end

axis equal; 
hold off;