function fraunhoffertest()
%% definitions

z=0.93*10^12;                    % distance à l'ecran en pm
n=2^13;                          % taille totale du reseau reseau  en pm
I_0=100;                         % intensité initiale
Energy=150;                      % energie du faisceau en eV
lambda=3.2*sqrt(80/Energy);      % longueur d'onde de de broglie en pm
resol=0.0005;                    % resolution angulaire en rad
m=1;                             % resolution du calcul (nbre de pm par pixels sur l'écran)

%% definition du reseau de graphene
Res=zeros(n);  % reseau graphene
D=426/m;       % distance entre 2 centres de mailles cotes a cote (haut et bas) 426 pm
d=246/m;       % distance entre 2 centres de mailles cotes a cote (droite et gauche) 246 pm

%% definition d un trou de passage des atomes
R=round(d/2);  % taille du cote des carres des trous (passage des atomes)
x=1:R;
y=1:R;
[X,Y]=meshgrid(x, y);
Diam=d/3;      % diametre des trous en pm

ROND=(sqrt((2*X/Diam-R/Diam).^2+(2*Y/Diam-R/Diam).^2))*0.4; %trou de passage de l hydrogene
ROND=(ones(R)-ROND);
ROND(ROND<0.6)=0;       %enleve la diminution des proba au bord du cercle
ROND(ROND>0.6)=1;       %enleve la diminution des proba au bord du cercle

%% construction du reseau de trous de graphene
for it1=1:d:n
    for jt1=1:D:n
    Res(it1:(it1+R)-1,jt1:(jt1+R)-1)=ROND;
    end
end
for it2=d/2:d:n
    for jt2=D/2:D:n
    Res(it2:(it2+R)-1,jt2:(jt2+R)-1)=ROND;
    end
end

%% fonction d'onde du faisceau d hydrogene
Res=Res(1:n,1:n);
a = (1:n)-n/2;
b = (1:n)'-n/2;
p=3;                               % proportionalité pour réduire la taille des matrices
sigma = lambda/sin(resol)/p; 
Res = Res.*exp(-(a.^2+b.^2)/(2*sigma^2))/(2*pi*sigma^2);

%% inclinaison du reseau
theta=0*pi/180;                    % angle d'inclinaison 
Res = Res.*exp(-i*2*pi*a*sin(theta)/lambda);

%% transformee de fourier
E=fftshift(fft2(fftshift(Res)));   % transformee de fourrier du reseau = figure de diffraction
I=conj(E).*E;

%% Mise a l echelle pour la transfo de fourier
 xdet = linspace(-(n/2), n/2, n)*lambda/n/m;        %en radiants
 ydet = linspace(-(n/2), n/2, n)*lambda/n/m;
%  
%  xdet = atan(xdet)*z;                             %en metres
%  ydet = atan(ydet)*z;

%% diffraction 2D
figure(2);
r=0.15;   %rectification des axes

imshow(I*I_0^2*10,'XData',xdet*10^3-r, 'YData',ydet*10^3-r);
xlabel('Distance (pm)', 'interpreter', 'latex');
ylabel('Distance (pm)', 'interpreter', 'latex');

% graph=pcolor(xdet*10^3-r, ydet*10^3-r,I*I_0^2*10);
% graph.LineStyle='none';
% pbaspect manual;
% pbaspect([1,1,1]);

xlabel('Angle (mrad)', 'interpreter', 'latex');
ylabel('Angle (mrad)', 'interpreter', 'latex');
axis on
colormap('gray')
axis([-100 100 -100 100]);
caxis([0 0.1])
%% reseau
% figure(3);
% Res=Res*2*pi*sigma^2;
% axesreseau=(-n/2:10:n/2)*m;
% imshow(I_0*Res,'XData', axesreseau, 'YData', axesreseau);
% xlabel('Distance (pm)', 'interpreter', 'latex');
% ylabel('Distance (pm)', 'interpreter', 'latex');
% axis on
 %caxis([0 0.01])
end