function [n rad N z] = landscapegeneration(n, rad)

N  = poissrnd(n);
r = randn(N,2);
z = r./repmat(sqrt(sum(r.^2,2)),1,2);
found = 0;
while ~found
    u = sqrt(rand(N,1));
    found = ~(sum(u < rad));
end
z = z.*repmat(u,1,2);

fig = figure;
hold on
plot([0,0],[0,0],'r*')
rectangle('Position',[-1,-1,2,2],'Curvature',[1,1],'EdgeColor','r','LineWidth',2)
for i=1:N
    rectangle('Position',[z(i,1) - rad ,z(i,2) - rad,  2*rad,   2*rad], ...
        'Curvature',[1,1],'FaceColor','k');
    plot(z(i,1),z(i,2),'.')
end
axis equal
hold off

end
