Squared magnitude of the scattered field versus observation angle
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.sstatic.net/3atos.png [3]: https://i.sstatic.net/Fk004.png
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');
- 10e5 == 10^5 is not true in MATLAB, why?
- 'fplot' plots hidden parameters figure in Matlab
- Using spfun with two identically ordered sparse matrices
- Scipy filter coefficients and filtering in MATLAB
- RANSAC using SIFT in Computer Vision
- Train Cascade Object Detector bounding box out of bound?
- vision.OpticalFlow is not showing the video nor the value of the flow
- Matlab Image Histogram Analysis: how do I test for an underlying bimodal distribution?
- MATLAB: Corr Function not working for Spearman
- Is there a programmatic way to generate s-functions directly out of Simulink blocks
- Matlab mldivide returns NaN when there are multiple solutions
- How do I access a MATLAB UI function in another .m file?
- Matlab: Is there a way to get the path of the current script?
- training neural network for user recognition in MATLAB
- how to make this function periodic in MATLAB?
- Automatically simplify redundant arithmetic relations
- Matlab fittype of a logged power law, with a y-intercept
- Matlab cross validation and K-NN
- KNN Classifier using cross validation
- Choosing k for KNN in Matlab
- LPR with MATLAB: how to find only one rectangle?
- Obtain origin of warning issued by m-functions with code generation
- numpy.dtype size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject
- Reverse 2-D table lookup
- Implementation of Absolute Differences Stereo Matching Algorithm
- Image Based Visual Servoing algorithm in MATLAB
- Reinforcement Learning
- Frangi filter working as an edge detector
- Matlab: test a class that takes user keyboard input
- Matlab code to generate the Phase and angle of the sinusoid?