clear all

%% 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=300;                      % 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 mrad
m=1;                             % resolution du calcul (nbre de pm par pixels sur l'écran)

%% definition du reseau de graphene
Res=zeros(n);  % reseau graphene
RES=zeros(n);  % reseau graphene bicouche
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/6;      % 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/2:n
    for jt1=1:D/2: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


%% bicouche AB
for it3=1:d:n
    for jt3=1:D:n
    RES(it3+d/4:(it3+R)-1+d/4,jt3+D/4:(jt3+R)-1+D/4)=ROND;
    end
end
for it4=d/2:d:n
    for jt4=D/2:D:n
    RES(it4+d/4:(it4+R)-1+d/4,jt4+D/4:(jt4+R)-1+D/4)=ROND;
    end
end
Res=Res(1:n,1:n);
RES=RES(1:n,1:n);
Res=RES+Res;

%% fonction d'onde du faisceau d hydrogene

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*5,'XData', xdet*10^3-r, 'YData', ydet*10^3-r);
xlabel('Angle (mrad)', 'interpreter', 'latex');
ylabel('Angle (mrad)', 'interpreter', 'latex');
axis on

%% 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])