
function plot_2D(domain_structure,pts1,pts2,pts3)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% 1. Plotting domain
% 2. Plotting pointset "pts1" (matrix N1 x 2), with small dots
% 3. Plotting pointset "pts2" (matrix N2 x 2), with larger magenta circles
% 4. Plotting pointset "pts3" (matrix N3 x 2), with larger blues asterisks
%--------------------------------------------------------------------------
% Note: 2,3,4 are not mandatory.
%--------------------------------------------------------------------------


if nargin < 2, pts1=[]; end
if nargin < 3, pts2=[]; end
if nargin < 4, pts3=[]; end

domain_str=domain_structure.domain;
bright_gray=[236, 236, 236]/256;

plot_dbox=0;

switch domain_str

    case {'polygon','square','rectangle','triangle',...
            'unit-square[0,1]x[0,1]','unit-simplex'}
       
        % ...................... decode variables .........................
        
        P=domain_structure.polyshape;
        plot(P,'FaceColor',bright_gray,'LineWidth',2);


    case {'NURBS','circular-annular-sector','sector',...
            'asymmetric-circular-sector','lune'}

        plot_NURBS(domain_structure);


%     case 'lune'
% 
%         centers=domain_structure.centers;
%         rV=domain_structure.radii;
% 
%         center1=centers(1,:); center2=centers(2,:);
%         r1=rV(1); r2=rV(2);
% 
%         [rM,thM]=meshgrid(0:0.01:r1, 0:0.01:2*pi);
%         rV=rM(:); thV=thM(:);
%         X0=center1(1)+rV.*cos(thV);
%         Y0=center1(2)+rV.*sin(thV);
% 
%         dist_2=sqrt((X0-center2(1)).^2+(Y0-center2(2)).^2);
% 
%         in=find(dist_2 >= r2);
% 
%         plot(X0(in),Y0(in),'.','Color',bright_gray,'LineWidth',2);
%         hold on; axis equal; 
%         plot_NURBS(domain_structure);
%         hold off;

    case 'disk'
        
        % ...................... decode variables .........................
        center=domain_structure.center;
        r=domain_structure.radius;

        th=linspace(0,2*pi,1000);
        xv=r*cos(th); xv(end+1)=xv(1);
        yv=r*sin(th); yv(end+1)=yv(1);

        plot(xv,yv,'k-','LineWidth',2);
        hold on;
        fill(xv,yv,bright_gray);
        hold off

    case 'asymmetric-annulus'
        
        % ...................... decode variables .........................

        centerV=domain_structure.centers;
        rV=domain_structure.radii;

        a=centerV(1,1); b=centerV(1,2); c=centerV(2,1); d=centerV(2,2);
        r1=rV(1); r2=rV(2);

        % ........................ plotting domain ......................

        A=[r1 0; r2 0]; B=[0 r1;0 r2]; C=[a b;c d];
        alpha=-pi; beta=pi;

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

    case 'vertical-circular-zone'

        center=domain_structure.center;
        r=domain_structure.radius;
        anglesV=domain_structure.angles;

        a=center(1); b=center(2);  
        alpha=anglesV(1); beta=anglesV(2);

        A=[r 0;r 0]; B=[0 r;0 -r]; C=[a b;a b];

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

   case 'horizontal-circular-zone'

        center=domain_structure.center;
        r=domain_structure.radius;
        anglesV=domain_structure.angles;

        a=center(1); b=center(2);  
        alpha=anglesV(1); beta=anglesV(2);

        A=[r 0;-r 0]; B=[0 r;0 r]; C=[a b;a b];

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

   case 'asymmetric-circular-sector'

        centerV=domain_structure.centers;
        r=domain_structure.radius;
        anglesV=domain_structure.angles;

        a=centerV(1,1); b=centerV(1,2); c=centerV(2,1); d=centerV(2,2);  
        alpha=anglesV(1); beta=anglesV(2);

        A=[r 0;0 0]; B=[0 r;0 0]; C=[a b;c d];

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


    case 'circular-segment'

        center=domain_structure.center;
        r=domain_structure.radius;
        anglesV=domain_structure.angles;

        a=center(1); b=center(2);
        alpha=anglesV(1);
        beta=anglesV(2);
        if beta < alpha, beta=beta+2*pi; end

        A=[r 0;r*cos(alpha+beta) r*sin(alpha+beta)];
        B=[0 r;r*sin(alpha+beta) -r*cos(alpha+beta)];
        C=[a b;a b];

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


    case 'symmetric-lens'

        % ....................... decode variables .......................
        a=domain_structure.center; a=abs(a);
        r=domain_structure.radius;

        if a > r
            error('MATLAB: "define_cub_rule"',getString(message(...
                'a must be strictly smaller than b')));
        end

        beta=acos(a/r); alpha=-beta;
        A=[r 0;-r 0]; B=[0 r;0 r]; C=[-a 0;a 0];
        plot_ellblend(A,B,C,alpha,beta);



    case 'candy'

        
        % ...................... decode variables .........................

        a=domain_structure.center;
        r=domain_structure.radius;
        alpha=domain_structure.angle;
        beta=-alpha;

        acos_a_r=acos(a/r);
        if beta <= acos_a_r
            error('MATLAB: "define_cub_rule"',getString(message(...
                'it must be "-alpha > acos(a/r)"')));
        end
        A=[r 0;-r 0]; B=[0 r;0 r]; C=[-a 0;a 0];
        plot_ellblend(A,B,C,alpha,beta);


   case 'butterfly'

        % "BUTTERFLY" domain: linear blending of two opposite circular arcs
        % of a disk with center in (a,b) and radius r
        % A=[r 0;-r 0]; B=[0 r;0 -r]; C=[a b;a b];
        %
        % In this case "parms" is a vector [a b r alpha beta]

        
        center=domain_structure.center;
        r=domain_structure.radius;
        anglesV=domain_structure.angles;

        a=center(1); b=center(2);
        alpha=anglesV(1); beta=anglesV(2);

        A=[r 0;-r 0]; B=[0 r;0 -r]; C=[a b;a b];

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


    case 'union-disks'
        
        % ...................... decode variables .........................

        centerV=domain_structure.centers;
        rV=domain_structure.radii;

        [centers,rs,active_arcs,disk_sequence,arcs_sizes]=...
            determine_uniondisks_boundary(centerV,rV);
        plot_figure_arcs(centers,rs,active_arcs,...
            disk_sequence,arcs_sizes);

    case {'polygcirc'}
       
        plot_polygcirc(domain_structure);

    otherwise
        warning('routine not implemented for this domain');

end

hold on;
axis equal;

dbox=domain_structure.dbox;

if isempty(dbox) == 0 & plot_dbox
    xLimit=dbox(1,:); yLimit=dbox(2,:);
    xm=xLimit(1); xM=xLimit(2); ym=yLimit(1); yM=yLimit(2);

    % platinum_RGB=[229, 228, 226]/256;
    bright_gray=[236, 236, 236]/256;
    plot([xm xM xM xm xm],[ym ym yM yM ym],':','color', bright_gray,...
        'LineWidth',1)
end


if not(isempty(pts1))
    plot(pts1(:,1),pts1(:,2),'o',...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor','k',...
        'MarkerSize',2);
end

if not(isempty(pts2))
    plot(pts2(:,1),pts2(:,2),'o',...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor','m',...
        'MarkerSize',6);
end

if not(isempty(pts3))
    plot(pts3(:,1),pts3(:,2),'*',...
        'MarkerEdgeColor','b',...
        'MarkerFaceColor','b',...
        'MarkerSize',4);
end

axis equal
%axis off
%set(gca,'xtick',[])
%set(gca,'ytick',[])
%title(domain_str);
hold off;
























