clc ; clear all ; close all;

figure;
a = 0;
b = 50;
pts=[50 100 200];
gap=[2 5];
lrmax=[1,1,1;0,0,0;0,0,0];
nmax=zeros(3,3);
expmax=zeros(3,3);
xcmax=zeros(3,3);
ycmax=zeros(3,3);
r0max=zeros(3,3);
r1max=zeros(3,3);
for lv1=1:3
    r = (b-a).*rand(pts(lv1),2) + a;
    for r0=5:((b/2)-2)
        for r1=(r0+2):(b/2)
            for xc=r1:(b-r1)
                for yc=r1:(b-r1)
                    % circles
                    th = linspace(0,2*pi) ;
                    x0 = r0*sin(th)+xc ; y0 = r0*cos(th)+yc ;
                    x1 = r1*sin(th)+xc ; y1 = r1*cos(th)+yc ;
                    
                    x = [x1 NaN fliplr(x0)] ;
                    y = [y1 NaN fliplr(y0)] ;
                    exp=pi*(r1^2-r0^2)*pts(lv1)/(b^2);
                    n = 0 ;   % number of points lying inside
                    for i=1:pts(lv1)
                        idx = inpolygon(r(i,1),r(i,2),x,y) ;
                        if idx
                            %plot(r(i,1),r(i,2),'.g')
                            n = n+1;
                        else
                            %plot(r(i,1),r(i,2),'.b')
                        end
                        %drawnow
                    end
                    lr=0;
                    if(n>exp)
                        lr=log(((n/exp)^n)*((pts(lv1)-n)/(pts(lv1)-exp))^(pts(lv1)-n));
                    end
                    if lr>lrmax(3,lv1)
                        lrmax(3,lv1)=lr;
                        nmax(lv1,3)=n;
                        expmax(lv1,3)=exp;
                        xcmax(lv1,3)=xc;
                        ycmax(lv1,3)=yc;
                        r0max(lv1,3)=r0;
                        r1max(lv1,3)=r1;
                    end
                    if lr<lrmax(1,lv1) && lr~=0
                        lrmax(1,lv1)=lr;
                        nmax(lv1,1)=n;
                        expmax(lv1,1)=exp;
                        xcmax(lv1,1)=xc;
                        ycmax(lv1,1)=yc;
                        r0max(lv1,1)=r0;
                        r1max(lv1,1)=r1;
                    end
                    if lr<1.005 && lr>0.995
                        lrmax(2,lv1)=lr;
                        nmax(lv1,2)=n;
                        expmax(lv1,2)=exp;
                        xcmax(lv1,2)=xc;
                        ycmax(lv1,2)=yc;
                        r0max(lv1,2)=r0;
                        r1max(lv1,2)=r1;
                    end
                end
            end
        end
    end
    
    subplot(4,3,lv1);
    hold on;
    for i=1:pts(lv1)
        plot(r(i,1),r(i,2),'.k','MarkerSize',10);
        drawnow
    end
    title(sprintf('No. of Points = %d',pts(lv1)));
    hold off;
    
    for lv2=1:3
        x0 = r0max(lv1,lv2)*sin(th)+xcmax(lv1,lv2) ;
        y0 = r0max(lv1,lv2)*cos(th)+ycmax(lv1,lv2) ;
        x1 = r1max(lv1,lv2)*sin(th)+xcmax(lv1,lv2) ;
        y1 = r1max(lv1,lv2)*cos(th)+ycmax(lv1,lv2) ;
        
        x = [x1 NaN fliplr(x0)] ;
        y = [y1 NaN fliplr(y0)] ;
        
        subplot(4,3,(lv2*3+lv1));
        hold on;
        
        plot(x0,y0,'-r','LineWidth',1.5);
        plot(x1,y1,'-r','LineWidth',1.5) ;
        
        for i=1:pts(lv1)
            idx = inpolygon(r(i,1),r(i,2),x,y) ;
            if idx
                plot(r(i,1),r(i,2),'.k','MarkerSize',10);
            else
                plot(r(i,1),r(i,2),'.k','MarkerSize',10);
            end
            drawnow
        end
        if(lv2==1)
            xlabel(sprintf('B=%2.0f, c=%d, LLR=%1.1e',expmax(lv1,lv2),nmax(lv1,lv2),lrmax(lv2,lv1)));
        else
            xlabel(sprintf('B=%2.3f, c=%d, LLR=%2.3f',expmax(lv1,lv2),nmax(lv1,lv2),lrmax(lv2,lv1)));
        end
        hold off;
        % savefig(sprintf('Ring_Pts=%d_LLR=%2.2f.fig',pts,lrmax));
    end
end
savefig(sprintf('Ring_Pts=%d_LLR=%2.2f.fig',pts(1),lrmax(3,1)));
%set(h,'facealpha',.5)