Squared magnitude of the scattered field versus observation angle

Question

I am trying to generate the figures in the paper [Intelligent Reflecting Surfaces at Terahertz Bands: Channel Modeling and Analysis][1]. Figure 2 is as below. [![enter image description here][2]][2]

What I get from my code is as below. [![enter image description here][3]][3]

The following is the code used.

% Intelligent reflecting surfaces at terahertz bands: channel modelling and
% analysis
% figure 2

close all
clear all
clc

% setting latex inerpreter
set(groot,'defaulttextinterpreter','latex'); 
set(groot,'DefaultTextFontname', 'CMU Serif');
set(groot,'DefaultAxesFontName', 'CMU Serif');
set(groot,'DefaultTextFontSize',14)
set(groot,'DefaultAxesFontSize',14)

% phi: azimuth angle, theta: polar angle
% t: transmitter, r: receiver
% d: distance from origin
theta_r = 0:.1:90;   % receiver polar angel
theta_t = 30;        % transmitter polar angle
phi_r = 60;          % receiver azimuth angle
% amplitude of electric (E-field) of the incident plane
amplitudeE_i_squared = 1;            
f = 300e9;          % carrier frequency
D_r = 4;            % distances measured from the (0,0)th IRS element to Rx
c = 3e8;            % velocity of light
lambda = c/f;       % wavelength

%% 1. lx = ly = lambda/2, eq. 9
L_x = lambda/2;     % size of reflecting element in x
L_y = lambda/2;     % size of reflecting element in y

X = ((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = ((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + ...
    sind(phi_r)^2);

E_s_NormSquared_eq9 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F .* sinc(X).^2 .* sinc(Y).^2);

plot(theta_r,E_s_NormSquared_eq9,'b');
hold on

%% 2. lx = ly = lambda/2, eq. 10
E_s_NormSquared_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F);

plot(theta_r,E_s_NormSquared_eq10,'k-.');

%% 3. lx = ly = 5*lambda, eq. 9
L_x = 5*lambda;
L_y = 5*lambda;

% lambda = 2*lambda;

X = ((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = ((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + ...
    sind(phi_r)^2);

E_s_NormSquared_eq9 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F .* ...
    sinc(X).^2 .* sinc(Y).^2);

% E_s_NormSquared_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
% (amplitudeE_i_squared/D_r^2) * F);

plot(theta_r,E_s_NormSquared_eq9,'r');
% plot(theta_r,E_s_NormSquared_eq10,'r');

xlabel('$\theta_r [degrees]$');
ylabel('$||E_s||^2 [dB]$');
hl = legend('$L_x = L_y = \lambda/2$ (Eq. 9)',...
    '$L_x = L_y = \lambda/2$ (Eq. 10)',...
    '$L_x = L_y = 5\lambda$');
set(hl, 'Interpreter','latex')
xlim([min(theta_r) max(theta_r)]);

% saveas(gca, ['f2_', datestr(now,'ddmmyyyy'), '.png']);

My plot and the z-axis range are not similar to the one in the paper.

Could you please point out the mistake I am doing. Thank you.

[1]: https://arxiv.org/abs/2103.15239#:~:text=Intelligent%20Reflecting%20Surfaces%20at%20Terahertz%20Bands%3A%20Channel%20Modeling%20and%20Analysis,-Konstantinos%20Dovelos%2C%20Stylianos&text=An%20intelligent%20reflecting%20surface%20(IRS,for%20the%20severe%20propagation%20losses. [2]: https://i.stack.imgur.com/3atos.png [3]: https://i.stack.imgur.com/Fk004.png

Answer 1

1.- Cardinal sin or sinc has 2 definitions :

While the obvious definition is sinc(x)=sin(x)/x in signal processing sinc is often defined as sinc(x)=sin(pi*x)/(pi*x)

Please understand that my main objective here was to get to the code producing the expected curves with certain coherence, but further refinement may be needed.

The following changes in the code supplied in the question bring all 3 graphs to the expected curves.

close all;clear all;clc

theta_r = [0:.1:90]; % [degree] receiver polar angle
theta_t =30;            % [degree] transmitter polar angle range
phi_r= 60;              % [degree] receiver azimuth angle range

% amplitude of electric (E-field) of the incident plane
Ei2 = 1;            
f = 300e9;              % [Hz]  carrier frequency
D_r = 4;                 % [m]  distances measured from the (0,0)th IRS element to Rx
c = 3e8;                  % [m/s] velocity of light
lambda = c/f;          % [m]  wavelength

%% 1. lx = ly = lambda/2, eq. 9
L_x = lambda/2;     % [m]  size of reflecting element in x 
L_y = lambda/2;     % [m] size of reflecting element in y

X = pi*((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = pi*((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + sind(phi_r)^2);

Es2_eq9 = 10*log10( ((L_x*L_y)/lambda^2)^2 * (Ei2/D_r^2) * F .*sinc(X).^2 .* sinc(Y).^2*1/pi^4);

plot(theta_r,Es2_eq9,'b');
hold on
axis([0 90 -140 -40])

%% 2. Lx = Ly = lambda/2, eq. 10
Es2_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * (Ei2/D_r^2) * abs(F));

plot(theta_r,Es2_eq10,'k-.');

%% 3. Lx = Ly = 5*lambda, eq. 9
N=50
L_x = lambda/N;
L_y = lambda/N;

X = pi*((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = pi*((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + sind(phi_r)^2);

Es2_eq9_2 = 10*log10(N*((L_x*L_y)/lambda^2)^2 * (Ei2/D_r^2) * F .* sinc(X).^2 .* sinc(Y).^21/pi^4);

plot(theta_r,Es2_eq9_2,'r');
grid on; grid minor;

xlabel('theta_r (degree)');ylabel('|E_s|^2 (dB)');

legend({'L_x = L_y = \lambda/2 (Eq. 9)',...
                  'L_x = L_y = \lambda/2 (Eq. 10)',...
                  ['L_x = L_y = lambda /' num2str(N) ]}, ...
                  'Location','northeast');

2.- X and Y modifyed

Heading pi to comply with sinc definition sinc(x)=sin(pix)/(pix)

A factor 1/pi^4 is added when calculating |Es2| in dB, again tot satisfy the sinc definition including pi.

If you would like to know why in signal processing sinc(x)=sin(pi*x)/(pi*x) please post another question.

3.- Electric field have to be [V/m]

Otherwise it's electric field no more.

Yet as defined in the paper, right after Figure.2

The units of Es are not [V/m] but [V/m^3]

May be the opening factor (Lx*Ly/lambda) should be (Lx*Ly/lambda^2) .

However it's more probable that the opening Lx Ly factor should be ((Lx/lambda)*(Ly/lambda)/lambda).

This squares electric field units to [V/m] .

5.- 3rd curve in paper may be wrong

When Lx Ly both rise to 5*lambda the resulting abs(Es(theta_r)) cannot remain as smooth as as compared to when Lx Ly a fraction of lambda.

The red graph in the paper seems closer to Lx Ly containing a fraction of lambda rather than several lambda.

So perhaps for the 3rd curve, either the 3rd graph in the paper Figure.2 is not correct, or instead of Lx=5*lambda Ly=5*lambda it should be Lx=lambda/50 Ly=lambda/50 then meeting Lx>>lambda and Ly>>lambda.

Consider that the 3rd curve on Figure.2 may not be correct.

close all;clear all;clc

theta_r = [0:.1:90]; % [degree] receiver polar angle
theta_t =30;            % [degree] transmitter polar angle range
phi_r= 60;              % [degree] receiver azimuth angle range

% amplitude of electric (E-field) of the incident plane
Ei2 = 1;            
f = 300e9;              % [Hz]  carrier frequency
D_r = 4;                 % [m]  distances measured from the (0,0)th IRS element to Rx
c = 3e8;                  % [m/s] velocity of light
lambda = c/f;          % [m]  wavelength

% 1. lx = ly = lambda/2, eq.9
L_x = lambda/2;     % [m]  size of reflecting element in x 
L_y = lambda/2;     % [m] size of reflecting element in y

X = pi*((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = pi*((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + sind(phi_r)^2);

Es2_eq9 = 10*log10( ((L_x*L_y)/lambda^2)^2 * (Ei2/D_r^2) * F .*sinc(X).^2 .* sinc(Y).^2*1/pi^4);

plot(theta_r,Es2_eq9,'b');
hold on; grid on
axis([0 90 -140 -40])

% 2. Lx = Ly = lambda/2 , eq.10
Es2_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * (Ei2/D_r^2) * abs(F));

plot(theta_r,Es2_eq10,'k-.');

% 3. Lx = Ly = 5*lambda , eq.10
N=1/5
L_x = lambda/N;
L_y = lambda/N;

X = pi*((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = pi*((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + sind(phi_r)^2);

Es2_eq9_2 = 10*log10(N*((abs(L_x*L_y)/lambda^2)^2 * (Ei2/D_r^2) * F .* sinc(X).^2 .* sinc(Y).^2*1/pi^4));

plot(theta_r,Es2_eq9_2,'r');

xlabel('theta_r (degree)');ylabel('|E_s|^2 (dB)');

legend({'L_x = L_y = \lambda/2 (Eq. 9)',...
            'L_x = L_y = \lambda/2 (Eq. 10)',...
            ['L_x = L_y =' num2str(1/N) ' lambda' ]}, ...
            'Location','northeast');