
Nfc=1024;
L=285;
D=24;
g2=g2norm(bb,{'rect',L,L,'radial'});

bb=(0:.01:15)*24;
db=mean(diff(bb));
dA=diff([0 (bb+db/2).^2]);

Nd=12;
del=2.^(0:Nd-1);

nn=zeros(Nd,length(bb));

Nc=length(sprintf('%d',Nf));
fmt=sprintf('[%%0%dd/%%0%dd]',Nc,Nc);

for idf=1:Nd
  df=del(idf);
  fprintf(1,fmt,0,Nfc);
  for nf=1:Nfc
    x1=repmat(px(1:Npf(nf),nf),1,Npf(nf+df));
    x2=repmat(px(1:Npf(nf+df),nf+df),1,Npf(nf));
    y1=repmat(py(1:Npf(nf),nf),1,Npf(nf+df));
    y2=repmat(py(1:Npf(nf+df),nf+df),1,Npf(nf));
    dx=x2-x1.';
    dy=y2-y1.';
    nn(idf,:)=nn(idf,:)+hist(abs(dx(:)+i*dy(:)),bb);
    if(mod(nf,100)==0);
      backspace(1,fmt,nf,Nfc);
      fprintf(1,fmt,nf,Nfc);
    end
  end
  backspace(1,fmt,Nfc,Nfc);
  fprintf(1,fmt,Nfc,Nfc);
  fprintf(1,'...Done.\n');
  plot(bb,nn./repmat(dA.*g2,Nd,1));
  axis([0 5*D 0 4e7]);
  drawnow;
end

