clc;
clear;

%Blazed grating metasurface simulation with not equidistant spacing
%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 along
%the x_1 axis, each with an adjustable phase shift and various spacing.

%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_quartz=1.46;   %Refrective index of quartz substrate at 800 nm
n_1=n_quartz;    %Refractive index of the first medium
n_2=1;           %Refrective index of the second medium (air)

lambda=0.8;     %wavelength [µm]
P=16;           %number of periodes in x-direction

%diameter_array=[400,350,300,250,200,150,100,50];
%Containing the diameters of the resonators employed for the spacing on the
%x_1 axis.
diameter_array=[]
for ind_M=1:P
    diameter_array=[diameter_array,[400,350,300,250,200,150,100,50]];
end

radial_array=diameter_array/2;
distance_array=100:100:length(diameter_array)*100; %additional distance by spacing d=100 nm between the resonators

Gamma_array=[2.6]; %Total latice constant L_x, periodicity of the phase shift pattern [µm]
                %The phase shift pattern repeats after this distance
Q=8;            %Number of elemental waves per periode, the phase shift is equidistantly distributed on these waves
p_B=-1;


% 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

ratio=0; %Not employed in this simulation

%% Rigorous

I_x_2_T_R_array=zeros(length(z0),length(x_2),length(Gamma_array));
I_x_2_R_R_array=zeros(length(z0),length(x_2),length(Gamma_array));

for m=1:length(Gamma_array)
    
    Gamma=Gamma_array(m);

    %The vector x0 contains the elemental waves origin positions, here it
    %is generated based on the spatial distribution of the resonators,
    %which may be non equidistant
    x0=zeros(1,length(radial_array));
    for ind=1:length(radial_array)
    x0(ind)=(radial_array(ind)+distance_array(ind)+sum(diameter_array(1:ind-1)))*10^(-3);
    end

    %The vector Phi0 contains the elemental waves phase offset (introduced by the antennas)
    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 (in this simulation not applicable)
    Phi0_expanded=zeros(1,P*(Q));
    subind=1;
    for ind=1:1:P*(Q)
        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));
    %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) 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/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));
    zM_reflection=-z0*ones(1,length(x_2));
    for k=1:length(x0)
        r_reflection(:,:,k)=sqrt(zM_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;
    
    %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


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

%% Transmission plot

fig=figure(7),clf
for m=1:length(Gamma_array)
    [Mval_transmission,Mind_transmission]=max(I_x_2_T_R_array(end,:,m));
    [Mval_reflection,Mind_reflection]=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_transmission),'-','LineWidth',1.0)
    hold on
    %plot([blazed_grating_angle_transmission(m)],[1],'vg','LineWidth',1.0)
    xlim([-60,60]);
    ylim([-0.1,1.1]);
    txt = ['$\Gamma=\unit[',num2str(Gamma_array(m)),']{\text{\textmu}m}$'];
    text(15,0.5,txt,'interpreter','none')
end
txt = ['$p_d=-2$'];
text(-50,0.2,txt,'interpreter','none')
txt = ['$p_d=-1$'];
text(-15,0.85,txt,'interpreter','none')
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_transmission,Mind_transmission]=max(I_x_2_T_R_array(end,:,m));
    [Mval_reflection,Mind_reflection]=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/Mval_reflection),'-r','LineWidth',1)
    hold on
    %plot([blazed_grating_angle_reflection(m)],[1],'vb')
    xlim([-60,60]);
    ylim([-0.1,1.1]);
    txt = ['$\Gamma=\unit[',num2str(Gamma_array(m)),']{\text{\textmu}m}$'];
    text(15,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
%save('Fig_S2_b_Nonlinear_phase_gradient_side_lobes.mat','ratio','Theta_m','n_1','n_2','lambda','P','Q','Gamma_array','p_B','xangle','I_x_2_T_R_array')







