% Programme permettant:
% 1 - de positionner des points aléatoirement sur un carré selon une loi
% uniforme;
% 2 - d'analyser la distribution précédente en comptant le nombre de points
% dans un disque de rayon h défini par l'utilisateur et de tracer
% un histogramme pour tout les disques de rayon h possibles
% 3 - de comparer la distribution précédente à celle de la loi de Poisson
% de même moyenne

% J.-C. POGGIALE -- Février 2007

% Initialisation
clear;
clc;
n=1000;     % Nombre de points à positionner
L=1;        % Taille du carré sur lequel les points sont positionnés
h=0.1;      % Taille du disque qui permet d'analyser la distribution des points
dh=0.005;   % Pas d'espace défini pour le déplacement du disque de rayon h
x=L*rand(n,1);  % Abscisses aléatoires des points
y=L*rand(n,1);  % Ordonnées aléatoires des points
u=zeros(n,1);   % Initialisation d'un vecteur contenant les distances entre les points et le centre du disque de rayon h
Nbh=fix((L-2*h)/dh); % Nombre de valeurs possibles pour les abscisses et les ordonnées du centre du disque de rayon h
X=zeros(Nbh*Nbh,1); %Initialisation d'un vecteur contenant le nombre d'occurences du nombre de points dans la boule de rayon h

% Comptage des points dans les boules de rayon h
for i=1:Nbh     % Variation de l'abscisse du centre du disque de rayon h
    for j=1:Nbh % Variation de l'ordonnée du centre du disque de rayon h
        xh=(h+(i-1)*dh)*ones(n,1);  % calcul de l'abscisse du centre du disque de rayon h
        yh=(h+(j-1)*dh)*ones(n,1);  % calcul de l'ordonnée du centre du disque de rayon h
        u=sqrt((x-xh).^2+(y-yh).^2); % Distance entre les points et le centre du disque de rayon h
        X(i+(j-1)*Nbh,1)=sum(u<=h); % Nombre de points dans le disque de rayon h
    end
end

% Graphiques
figure(1) 
plot(x,y,'.') % Représentation des points dans le carré de côté L

kmax=max(X); % Nombre maximum de points dans une boule
k=0:kmax;   % Vecteur (0,1,2,...,kmax)
figure(2)   
y=hist(X,k); % y= l'histogramme des occurences 
xlim([0 kmax]);
hist(X) % Tracé de l'histogramme des occurences 


m=mean(X); % Calcul de la moyenne du nombre de points dans un disque de rayon h
for k=0:kmax;
    p(k+1,1)=exp(-m)*m^k/(factorial(k));  % Distribution de Poisson de moyenne m
end
figure(3)
plot(p,p,p,y/(Nbh^2),'+') % Comparaison entre l'histogramme expérimental et la distribution de Poisson
xlabel('Probabilité théorique')
ylabel('Proportion observée')
