function etalement_gaussien(Nph) % Nph: nbre de particule
%
% Cette fonction permet de vérifier l'étalement Gaussien des particules après une marche aléatoire
% Le début reprend le script densite_des_particules.m, légèrement modifié pour supprimer le film et l'affichage
% des position à chaque cycle, et le nombre de pas est diminué, tout ça pour gagner du temps.
close all
N=100; %nbre de pas, doit être >> 1 pour être en régime diffusif.
%% Marche aléatoire
x=zeros(Nph,N);
y=zeros(Nph,N); 
for jj=1:Nph % boucle sur les particules
for ii=2:N % boucle sur le temps
teta=2*pi*rand;
x(jj,ii)= x(jj,ii-1)+cos(teta);
y(jj,ii)= y(jj,ii-1)+sin(teta);
end
end
% Ici on enregistre toutes les positions intermédiaires, mais on ne va regarder que la distribution à la fin
X = x(:,N);
Y = y(:,N);
clear x y
figure(1), plot(X,Y, '.'), xlabel('x/l'), ylabel('y/l'), title(['Position des ' num2str(Nph) ' particules apres ' num2str(N) ' pas']);
%% Construction d'une image de la densité locale
% Pour vérifier que c'est gaussien, la première étape est de construire un histogramme 2D, càd une image 
représentant la % densité locale de particule. Il faut définir des pixels et compter le nombre de particule dans chaque pixel.
R = [X,Y];
N_hist = 80; % Nombre de points dans l'image N_hist x N_hist
[densite, C] = hist3(R, [N_hist, N_hist]); % la fonction hist3 fait le travail de comptage; 
xx = C{1};
yy = C{2};
IMG = flipud(rot90(densite)); % Construction d'une image avec la bonne orientation
figure(2), imagesc(xx, yy, IMG), colorbar, xlabel('x/l'), ylabel('y/l'), title('Densite de particules (N/pixel)');
%% Maintenant il faut faire un ajustement gaussien.
pix_x = xx(2)-xx(1); % pixel size
pix_y = yy(2)-yy(1);
x=0:(N_hist-1);
y=0:(N_hist-1);
[XX,YY] = meshgrid(x,y);
%% Il faut des valeurs de départ ("guess")
A_g = max(max(IMG)); % amplitude
[y0_g, x0_g] = find(IMG==A_g); % centre
sigmaX_g = sqrt( mean(mean((XX-x0_g(1)).^2.*IMG))); % standard deviation
sigmaY_g = sqrt( mean(mean((YY-y0_g(1)).^2.*IMG))); 
GUESS = [x0_g(1) y0_g(1) A_g sigmaX_g/4 sigmaY_g/4]; % Guess
%% Fit
SOLUS = fminsearch(@gauss2D, GUESS, [], XX, YY, IMG);
[khi2, FIT] = gauss2D(SOLUS, XX, YY, IMG);
x0 = round(SOLUS(1)); % pour faire une coupe
y0 = round(SOLUS(2));
%% Affichage du résultat
figure(3), 
subplot(2,3,1), imagesc(xx, yy, IMG), colorbar, title('Data'), xlabel('x/l'), ylabel('y/l');
subplot(2,3,2), imagesc(xx, yy, FIT), colorbar, title('Fit'), xlabel('x/l'), ylabel('y/l');
subplot(2,3,3), imagesc(xx, yy, IMG-FIT), colorbar, title('Residu'), xlabel('x/l'), ylabel('y/l');
subplot(2,3,4), plot(xx, IMG(y0,:), '.r', xx, FIT(y0,:), 'k'), ...
xlabel('x (mm)'), ylabel('Intensity (u.a.)'), title('Cut along y_0');
subplot(2,3,5), plot(yy, IMG(:,x0), '.r', yy, FIT(:,x0), 'k'), ...
xlabel('y (mm)'), ylabel('Intensity (u.a.)'), title('Cut along x_0');
subplot(2,3,6), axis off, text(0.2, .8, ['\sigma_x/l = ' num2str(SOLUS(4)*pix_x)]),...
text(0.2, 0.6, ['\sigma_y/l = ' num2str(SOLUS(5)*pix_y)]), ...
text(0.2, 0.4, ['x_0/l = ' num2str(SOLUS(1)*pix_x+xx(1))]),...
text(0.2, 0.24, ['y_0 /l = ' num2str(SOLUS(2)*pix_y+yy(1))]);


%%%%%%%%%%%%%%%% Fonction de fit
function [khi2, FIT] = gauss2D(COEFF, X, Y, DATA)
x0 = COEFF(1);
y0 = COEFF(2);
A = COEFF(3);
sigmaX = COEFF(4);
sigmaY = COEFF(5);
FIT = A * exp(-(X-x0).^2/(2*sigmaX^2)-(Y-y0).^2/(2*sigmaY^2) ) ; % Statistical gaussian function
% Function to be minimised :
khi2 = mean(mean ( (FIT-DATA).^2 ));

                  
                  

                    
                      
% déclaration + initialisation 
clear all
close all
N=1000; %nbre de pas
Nph=10000; %nbre de photon
t=[1:N]; % le temps
%d=2; %nbre de dimension
% teta=0; % angle d'émission
% r=zeros(N,1); % vecteur position
% moypos=zeros(N,1); % position moyenne
R2=zeros(Nph,N);
for jj=1:Nph
x=0; y=0;
for ii=1:N
teta=2*pi*rand;
x= x+cos(teta);
y= y+sin(teta);
R2(jj,ii) = x*x + y*y;
end
end
moypos= mean(R2);
%moypos(ii)= sum(r)/ii;
plot(t,moypos);
xlabel('t');
ylabel('moypos');