clc;
clear;

%Blazed grating metasurface simulation

%Provided by Markus Schake as supplemental material for the JOSA A article
%"Examining and explaining the “Generalized Laws of Reflection and Refraction” at metasurface gratings"

%This simulation covers the interference image, generated by the
%superposition of Q*P elemental waves which are generated equidistantly along
%the x_1 axis, each with an adjustable phase shift.

%The intensity distribution is recorded on a virtual detector on the x_2 axis a
%specified distance z away and parallel to the x_1 axis.

%%%%%%%%%%%%%%%%%%%%%
%%
%Parameters

Theta_m=0;  %The angle of the incidence wave in degree,
%           it will manipulate the phase offset of the respective elemental
%           waves along x_1 direction.

n_si=3.4;   %Estimated refrective index of Silicon at a wavelength of 8 µm
n_1=n_si;   %Refractive index of the 1. medium
n_2=1;      %Refrective index of the second medium (air)

lambda=8;   %wavelength [µm]
P=16;       %number of periodes in x_1 direction

Gamma_array=[11,13,15,17,19,21]; %grating period of the phase shift pattern in [µm]
%The phase shift pattern repeats periodically after this distance
Q=8;            %Number of elemental waves per periode, the phase shift is equidistantly distributed on these waves
p_B=-1;          %Design diffraction order



% z0=lambda*logspace(0,6,100)'; %Positions of the virtual detector on the z-axis
z0=lambda*10^6; %Positions of the virtual detector on the z axis (just the last position on the screen is employed)

x_2=linspace(-11^(6)*lambda,11^(6)*lambda,2000); % x_2 positions to evaluate in the distance z0

%Incorporate Z points, which are not occupied by the antennas and therefore
%act abiding basic Snell's law
ratio=2;
Z=ratio*Q; %number of unaffected elemental waves


%% Rigorous


I_x_2_T_R_array=zeros(length(z0),length(x_2),length(Gamma_array)); %Intensity distribution on the screen for rigorous simulation of the transmission grating
I_x_2_R_R_array=zeros(length(z0),length(x_2),length(Gamma_array)); %Intensity distribution on the screen for rigorous simulation of the reflection grating

for m=1:length(Gamma_array)
    
    Gamma=Gamma_array(m);
    
    %The vector x0 contains the elemental waves origin positions
    x0=linspace(0,P*Gamma,P*(Q+Z));
    
    %The vector Phi0 contains the elemental waves phase offset (introduced by the antennas)
    %assuming a linear phase distribution along the whole grid
    Dphi_Gamma=linspace(0,p_B*2*pi+(2*pi/Q),Q); %Phase shift on one period Gamma
    Phi0=[];
    for indm=1:P
        Phi0=[Phi0,Dphi_Gamma];
    end
    
    
    %Include the elementary waves unaffected by the resonators
    Phi0_expanded=zeros(1,P*(Q+Z));
    subind=1;
    for ind=1:ratio+1:P*(Q+Z)
        Phi0_expanded(ind)=Phi0(subind);
        subind=subind+1;
    end
    
    
    %Phi_Thetai contains the incidence angle specific phase offset for the
    %respective elemental waves, considering zero offset for the first
    %elemental wave.
    Phi_thetai=linspace(0,sind(Theta_m)*P*Gamma*(2*pi*n_1/lambda),P*(Q+Z));
    %sind(Theta_m)*P*Gamma is the maximum OPLD caused by the angle of
    %incidence as can be deducted from the triangle between the grating,
    %screen and distance z.
    
    
    %Calculate the intensity of the interferogram assuming coherent illumination by a laser source.
    %Only the dependency of the electric field is considered for the
    %simulation.
    
    E0=1;   %Wave Amplitude;
    
    %% Transmission
    
    %The field is calculated for each (z0,x_2) position considering each source x0
    r_transmission=zeros(length(z0),length(x_2),length(x0));
    
    %Matrix containing the z0 value for each x position
    zM=z0*ones(1,length(x_2));
    xM=ones(length(z0),1)*x_2;
    
    for k=1:length(x0)
        r_transmission(:,:,k)=sqrt(zM.^2+(xM-x0(k)).^2);
    end
    
    E_transmission=E0*exp(j*(2*pi*n_2/lambda)*r_transmission);
    
    %Add The phase offset to the respective Elemental waves
    for k=1:length(Phi0_expanded)
        E_transmission(:,:,k)=E_transmission(:,:,k)*exp(j*(Phi0_expanded(k)+Phi_thetai(k)));
    end
    
    Esum_transmission=sum(E_transmission,3);
    
    I_transmission=Esum_transmission.*conj(Esum_transmission);
    
    %% Reflection
    
    r_reflection=zeros(length(z0),length(x_2),length(x0));
    z_reflection=-z0*ones(1,length(x_2));
    for k=1:length(x0)
        r_reflection(:,:,k)=sqrt(z_reflection.^2+(xM-x0(k)).^2);
    end
    E_reflection=E0*exp(j*(2*pi*n_1/lambda)*r_reflection);
    %Add The phase offset to the respective Elemental waves
    for k=1:length(Phi0_expanded)
        E_reflection(:,:,k)=E_reflection(:,:,k)*exp(j*(Phi0_expanded(k)+Phi_thetai(k)));
    end
    Esum_reflection=sum(E_reflection,3);
    I_reflection=Esum_reflection.*conj(Esum_reflection);
    %%%%%%%%%%%%%%%%%%%%
    
    %Save arrays
    I_x_2_T_R_array(:,:,m)=I_transmission;
    I_x_2_R_R_array(:,:,m)=I_reflection;
    
end


%Scale x_2-axis as angle representation
xangle=atan2(x_2,max(z0))*180/pi;


%% Analytic

%Simulation employing analytic grid equation and Fraunhofer diffraction.
%The analytical solution employs the frequency shifting
%theorem of the fourier transform and thus requires a linear phase
%distribution of constant slope.

z=max(z0);             %Observation distance
Gamma=Gamma_array;     %Grating period
lambda_B=1.0*lambda;   %Design wavelength
w=(max(x0)-min(x0));   %Total width of the grid

theta_m_rad=Theta_m*pi/180;      %Angle of incidence

x_2_analytic=linspace(x_2(1),x_2(end),10000); %Lateral screen coordinate


for m=1:length(Gamma)
    
    %% Transmission Analytic
    
    korder=-10:10;
    for k=1:length(korder)
        xi_analytic_transmission=x_2_analytic/((lambda/n_2)*z);
        x_2_transmission=xi_analytic_transmission*(lambda/n_2)*z;
        argsinc_motif_transmission=pi*Gamma(m)*((x_2_transmission/((lambda/n_2)*z))-(p_B*(lambda_B/n_2)/((lambda/n_2)*Gamma(m))+theta_m_rad*n_1/lambda));
        sinc_motif_transmission=sin(argsinc_motif_transmission)./argsinc_motif_transmission;
        argsinc_interferencefunction_transmission=pi*w*((x_2_transmission/((lambda/n_2)*z))-(korder(k)/Gamma(m)+theta_m_rad*n_1/lambda));
        sinc_interferencefunction_transmission(:,k)=sin(argsinc_interferencefunction_transmission)./argsinc_interferencefunction_transmission;
    end
    I_x2_analytic_transmission(m,:)=(sum(sinc_interferencefunction_transmission.*sinc_motif_transmission',2)).^2;
    
    
    %% Reflection Analytic
    
    korder=-10:10;
    for k=1:length(korder)
        xi_analytic_reflection=x_2_analytic/((lambda/n_1)*z);
        x_2_reflection=xi_analytic_reflection*(lambda/n_1)*z;
        argsinc_motif_reflection=pi*Gamma(m)*((x_2_reflection/((lambda/n_1)*z))-(p_B*(lambda_B/n_1)/((lambda/n_1)*Gamma(m))+theta_m_rad*n_1/lambda));
        sinc_motif_reflection=sin(argsinc_motif_reflection)./argsinc_motif_reflection;
        argsinc_interferencefunction_reflection=pi*w*((x_2_reflection/((lambda/n_1)*z))-(korder(k)/Gamma(m)+theta_m_rad*n_1/lambda));
        sinc_interferencefunction_reflection(:,k)=sin(argsinc_interferencefunction_reflection)./argsinc_interferencefunction_reflection;
    end
    I_x2_analytic_reflection(m,:)=(sum(sinc_interferencefunction_reflection.*sinc_motif_reflection',2)).^2; %Coherent superposition
    
    %%
    
    %For comparison with the transmission/reflection law from [1]
    regular_snell_angle_transmission(m)=real(asind(sind(Theta_m)*n_1/n_2));
    
    blazed_grating_angle_transmission(m)=real(asind(sind(Theta_m)*n_1/n_2+p_B*lambda/(Gamma_array(m)*n_2)));
    
    blazed_grating_angle_reflection(m)=real(asind(sind(Theta_m)+p_B*lambda/(Gamma_array(m)*n_1)));
    
end


%% Transmission plot

fig=figure(7),clf
for m=1:length(Gamma_array)
    [Mval,Mind]=max(I_x_2_T_R_array(end,:,m));
    [Mvalminus,Mindminus]=max(I_x_2_R_R_array(end,:,m));
    subplot(length(Gamma_array),1,m)
    plot(xangle,I_x_2_T_R_array(end,:,m)*(1/Mval),'-','LineWidth',1.0)
    hold on
    plot((x_2_analytic/z)*180/pi,I_x2_analytic_transmission(m,:)*1,'--k','LineWidth',1.0)
    plot([blazed_grating_angle_transmission(m)],[1],'vg','LineWidth',1.0)
    xlim([-60,60]);
    ylim([-0.1,1.1]);
    txt = ['$\Gamma=\unit[',num2str(Gamma(m)),']{\text{\textmu}m}$'];
    text(20,0.5,txt,'interpreter','none')
    
end
han=axes(fig,'visible','off');
han.Title.Visible='on';
han.XLabel.Visible='on';
han.YLabel.Visible='on';
xlabel(han,'$\theta_{m+1}$/(deg)','interpreter','none');
ylabel(han,'$\frac{I(x_2)_T}{\text{max}\left(I(x_2)_T \right)}$','interpreter','none');

%% Reflection plot

fig=figure(8),clf
for m=1:length(Gamma_array)
    [Mval,Mind]=max(I_x_2_T_R_array(end,:,m));
    [Mvalminus,Mindminus]=max(I_x_2_R_R_array(end,:,m));
    subplot(length(Gamma_array),1,m)
    plot(xangle,I_x_2_R_R_array(end,:,m)*(1/Mvalminus),'-r','LineWidth',1)
    hold on
    plot((x_2_analytic/z)*180/pi,I_x2_analytic_reflection(m,:)*1,'--k')
    plot([blazed_grating_angle_reflection(m)],[1],'vb')
    xlim([-60,60]);
    ylim([-0.1,1.1]);
    txt = ['$\Gamma=\unit[',num2str(Gamma(m)),']{\text{\textmu}m}$'];
    text(20,0.5,txt,'interpreter','none')
end
han=axes(fig,'visible','off');
han.Title.Visible='on';
han.XLabel.Visible='on';
han.YLabel.Visible='on';
xlabel(han,'$\grave{\theta_{m}}$/(deg)','interpreter','none');
ylabel(han,'$\frac{I(x_2)_R}{\text{max}\left(I(x_2)_R \right)}$','interpreter','none');

%% Save your results in a *.m file
%xangle_analytic=(x_2_analytic/z)*180/pi;
%save('Fig_S1_a_Gratingperiod_variation.mat','ratio','Theta_m','n_1','n_2','lambda','P','Q','Gamma_array','p_B','xangle','I_x_2_T_R_array','xangle_analytic','I_x2_analytic_transmission','blazed_grating_angle_transmission')

%Imortant notes:
%The analytical model employs the paraxial approximation and is therefore
%unreliable for high angles of diffraction.

%The sinc(x) function is defined as sin(pi*x)/pi*x


