clear all

%% definitions
data = load('alldata.txy');                     %donnees experimentales
 x= data(1:size(data,1),1)'-2;
 y= data(1:size(data,1),2)'-5; 
 z=0.93*10^3;                                   % distance detecteur graphène en mm
 
% figure(1)                                       %representation brut des donnees experimentales
%  plot(x,y,'.')              
 
%% passage coord polaires
 [theta,R] = cart2pol(x,y);
 
%% passage en coord cartesiennes
[X,Y]=pol2cart(theta,R);

%% enleve le bord du detecteur
   theta=theta(R<=330);
    R = R(R<=330);
    
 %% histogramme 2D du signal
A_t=theta(R>=50);
A_R=R(R>=50);
[X_A,Y_A]=pol2cart(A_t,A_R);        %ordre 1 de diffraction

B_t=theta(R<50);
B_R=R(R<50);
[X_B,Y_B]=pol2cart(B_t,B_R);        %ordre 0 de diffraction (tache centrale)

edges_x=[-500:2:500];   
edges_y=[-500:2:500];

nb_A=histcounts2(X_A,Y_A,edges_x,edges_y);
nb_B=histcounts2(X_B,Y_B,edges_x,edges_y);
nb_tot=nb_A+nb_B/500;

X_axe = (edges_x(1:end-1)+edges_x(2:end))/2;
Y_axe = (edges_y(1:end-1)+edges_y(2:end))/2;

[bin_X,bin_Y]=meshgrid(X_axe,Y_axe);
bin_X=bin_X*40/900;                             % transformation de px à mm
bin_Y=bin_Y*40/900;
bin_X=atan(bin_X/z)*10^3;                       % transformation de mm à mrad
bin_Y=atan(bin_Y/z)*10^3;

figure(2)
hold on
graph = pcolor(bin_X,bin_Y,nb_tot);
colormap(bone(10));
graph.LineStyle = 'none';
pbaspect manual;
pbaspect([1,1,1]);
axis([-20 20 -20 20])
xlabel('Diffraction angle (mrad)', 'interpreter', 'latex');
ylabel('Diffraction angle (mrad)', 'interpreter', 'latex');
hold off

%% enleve grosse tache centrale
%    theta=theta(R>=50);
%     R = R(R>=50);


%% histogramme des coups en fonction du rayon
    R=R(2:end);
    theta=theta(2:end);
    
edges=[0:2:1000];
R_axes = (edges(1:end-1)+edges(2:end))/2;

 nb_R = histcounts(R,edges);
 err_R = sqrt(nb_R);

 
 R_axes = R_axes(nb_R>0);
 err_R  = err_R(nb_R>0);
 nb_R   = nb_R(nb_R>0);
 
 R_axes=R_axes*40/900;
 R_axes=atan(R_axes/z)*10^3;
 
 figure(3)
 errorbar(R_axes, nb_R./R_axes, err_R./R_axes, '.');
 axis([0 16 200 500])
 xlabel('$R$ (mrad)', 'interpreter', 'latex');
 ylabel('Counts', 'interpreter', 'latex');
 
%% calcul de la diffusion
Dt1=theta(R<=320);                  %tranche de cercle apres la diffraction
DR1=R(R<=320);
Dt1=Dt1(DR1>260);                  %tranche de cercle apres la diffraction
DR1=DR1(DR1>260);

Dt2=theta(R<=200);                  %tranche de cercle avant la diffraction
DR2=R(R<=200);
Dt2=Dt2(DR2>140);                  %tranche de cercle avant la diffraction
DR2=DR2(DR2>140);
    
Dt3=theta(R<=260);                  %tranche de cercle de la diffraction
DR3=R(R<=260);
Dt3=Dt3(DR3>200);                  %tranche de cercle de la diffraction
DR3=DR3(DR3>200);

%% histogramme des coups en fonction de l'angle sans diffusion
edges=[-pi:0.1:pi];
R_axes = (edges(1:end-1)+edges(2:end))/2;

 nb_t1 = histcounts(Dt1,edges);
 nb_t2 = histcounts(Dt2,edges);
err_t1 = sqrt(nb_t1);
err_t2 = sqrt(nb_t2);

  nb_tdiffract = histcounts(Dt3,edges);
err_tdiffract = sqrt(nb_tdiffract);
 figure(4)
 hold on
 %diffraction sans diffusion
plot(R_axes*180/pi, nb_tdiffract-(nb_t1+nb_t2)/2, 'Color','black');
  xlabel('Angle $\theta(^\circ$) ', 'interpreter', 'latex');
 ylabel('Nombre de coups', 'interpreter', 'latex');
 axis([-180 180 -100 1000])
 %diffusion moyenne
% plot(R_axes, (nb_t1+nb_t2)/2, 'b');

for j=1:6                                   %taches d'hexagone les plus brillantes
A(j,:)=zeros(1,100)-2.39+pi*(j-1)*60/180;
end
for j=1:6
B(j,:)=(1:100)*9;
end
line(A'*180/pi,B','Color','red')

for j=1:6                                   %taches d'hexagone les moins brillantes
C(j,:)=zeros(1,100)-2.69+pi*(j-1)*60/180;
end
for j=1:6
D(j,:)=(1:100)*9;
end
line(C'*180/pi,D','Color','blue')

 hold off