Continue to Site

# Solving frequency selective surfaces using method of moments on matlab

Status
Not open for further replies.

#### ahmedatef0

##### Junior Member level 1
I have been trying to simulate a narrow strip fss (frequency selective surface) using matlab, but I'm getting wrong results can someone help me understand why.
Code:
%% Initialization
clear;
clc;
close all;

%% Single element information
% Strip Size 15mm*0.5mm
% Cell Size 7.5mm*20
Dimensions.strip = [0.5 15]*1e-3;
Dimensions.dx = 1e-3;
Dimensions.dy = 20*1e-3;
%% Constants and Sweep Variables
Constants.eps = 8.854187812813*1e-12;
Constants.mu = 1.25663706*1e-6;
Constants.eta = sqrt(Constants.mu/Constants.eps);
Constants.c = 1/sqrt(Constants.mu*Constants.eps);
Constants.M = 30;
Constants.N = 30;
w = 10:0.5:20;
w = 2*pi*w*1e9;
%% Basis Data
% Basis.size = [0.5 2.5/2]*1e-3;
% Basis.numOfBasis = [1 0.03/Basis.size(2)-1];
Basis.numOfBasis = [1 20];
Basis.size = [0.5*1e-3 2*Dimensions.strip(2)/(1+Basis.numOfBasis(2))];
%% Incidence Data
Incidence.phi = 0;
Incidence.theta = pi/4;
Incidence.k = w/Constants.c;
Incidence.kx = -Incidence.k*sin(Incidence.theta)*cos(Incidence.phi);
Incidence.ky = -Incidence.k*sin(Incidence.theta)*sin(Incidence.phi);
%% Incidence Fields
% TM
Incidence.E_tm = 1i*w;%*cos(Incidence.theta);
Incidence.E_tm = [cos(Incidence.theta)*cos(Incidence.phi);cos(Incidence.theta)*sin(Incidence.phi)]...
*Incidence.E_tm;
% TE
Incidence.E_te = -1i*w*Constants.eta;%*cos(Incidence.theta);
Incidence.E_te = [sin(Incidence.phi);-cos(Incidence.phi)]...
*Incidence.E_te;
%% Basis Function
% Fourier of pulse of amplitude = 1 from -width/2 to width/2
% Fourer of Triangle of amplitude = 1 from -length/2 to length/2
% Basis(kx,ky) = width*sinc(kx*width/2)*(length/2)*sinc^2(ky*length/4)
Incidence.Ey = Incidence.E_tm(2,:)+Incidence.E_te(2,:);

I = zeros(Basis.numOfBasis(2),1,length(w));

%% Create Voltage vector
V = Basis.size(1)*sinc(Incidence.kx*Basis.size(1)/2);
V = V*(Basis.size(2)/2).*(sinc(Incidence.ky*Basis.size(2)/4)).^2;
V = V.*Incidence.Ey;
V = reshape(V,[1 1 length(w)]);
V = repelem(V,Basis.numOfBasis(2),1,1);
% Vshift & vshift2 are used for the exponenetial caused y the phase shift
Vshift = Incidence.ky;
Vshift = reshape(Vshift,[1 1 length(w)]);
Vshift = repelem(Vshift,Basis.numOfBasis(2),1,1);
Vshift2 = 1:Basis.numOfBasis(2);
Vshift2 = reshape(Vshift2,[length(Vshift2) 1 1]);
Vshift2 = repelem(Vshift2,1,1,length(w));
V = V.*exp(1i*Vshift.*Vshift2*Basis.size(2)/2);

%% Create Z matrix
% kyn = ky_inc+2npi/dy
%                           p q   freq.   m n
ky = reshape(Incidence.ky,[1 1 length(w) 1 1]);
ky = repelem(ky,Basis.numOfBasis(2),Basis.numOfBasis(2),1,Constants.M+1,Constants.N+1);
periodOfN = -Constants.N/2:Constants.N/2;
ky2 = 2*pi*periodOfN/Dimensions.dy;
ky2 = reshape(ky2,[1 1 1 1 length(periodOfN)]);
ky2 = repelem(ky2,Basis.numOfBasis(2),Basis.numOfBasis(2),length(w),Constants.M+1,1);
kyn = ky + ky2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% kxm = kx_inc+2mpi/dx
%                           p q   freq.   m n
kx = reshape(Incidence.kx,[1 1 length(w) 1 1]);
kx = repelem(kx,Basis.numOfBasis(2),Basis.numOfBasis(2),1,Constants.M+1,Constants.N+1);
periodOfM = -Constants.M/2:Constants.M/2;
kx2 = 2*pi*periodOfM/Dimensions.dx;
kx2 = reshape(kx2,[1 1 1 length(periodOfM) 1]);
kx2 = repelem(kx2,Basis.numOfBasis(2),Basis.numOfBasis(2),length(w),1,Constants.N+1);
kxm = kx + kx2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% basis
%                           p q   freq.   m n
basis_x = Basis.size(1)*sinc(kxm*Basis.size(1)/2);
basis_x = basis_x.^2;
basis_y = (Basis.size(2)/2)*(sinc(kyn*Basis.size(2)/4)).^2;
basis_y = basis_y.^2;
basis = basis_x.*basis_y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Green's Function
%                           p q   freq.   m n
k = reshape(Incidence.k,[1 1 length(w) 1 1]);
k = repelem(k,Basis.numOfBasis(2),Basis.numOfBasis(2),1,Constants.M+1,Constants.N+1);
g = 0.5./sqrt(kxm.^2+kyn.^2-k.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase Shift
%                           p q   freq.   m n
P = 1:Basis.numOfBasis(2);
Q = P;
P = reshape(P,[length(P) 1 1 1 1]);
P = repelem(P,1,Basis.numOfBasis(2),length(w),Constants.M+1,Constants.N+1);
Q = reshape(Q,[1 length(Q) 1 1 1]);
Q = repelem(Q,Basis.numOfBasis(2),1,length(w),Constants.M+1,Constants.N+1);
shift = exp(-1i*kyn.*(Q-P).*Basis.size(2)/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constants
Z = 1./(Dimensions.dx*Dimensions.dy*1i*w*Constants.eps);
%              p q   freq.   m n
Z = reshape(Z,[1 1 length(w) 1 1]);
Z = repelem(Z,Basis.numOfBasis(2),Basis.numOfBasis(2),1,Constants.M+1,Constants.N+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adding it all
Z = Z.*(k.^2-kyn.^2).*g.*basis.*shift;
Z = sum( sum( Z , 5 ) , 4 );
for i = 1:length(w)
I(:,:,i) = Z(:,:,i)\V(:,:,i);
end

%% Reflected Floquet Harmonics
% M*N*omega*P
%Constant
Esmn = 1./(Dimensions.dx*Dimensions.dy*1i*w*Constants.eps);
Esmn = reshape(Esmn,[1 1 length(w) 1]);
Esmn = repelem(Esmn,Constants.M+1,Constants.N+1,1,Basis.numOfBasis(2));
% k
k = reshape(Incidence.k,[1 1 length(w) 1]);
k = repelem(k,Constants.M+1,Constants.N+1,1,Basis.numOfBasis(2));
% kxm
kx = reshape(Incidence.kx,[1 1 length(w) 1]);
kx = repelem(kx,Constants.M+1,Constants.N+1,1,Basis.numOfBasis(2));
periodOfM = -Constants.M/2:Constants.M/2;
kx2 = 2*pi*periodOfM/Dimensions.dx;
kx2 = reshape(kx2,[length(periodOfM) 1 1 1]);
kx2 = repelem(kx2,1,Constants.N+1,length(w),Basis.numOfBasis(2));
kxm = kx + kx2;
% kyn
ky = reshape(Incidence.ky,[1 1 length(w) 1]);
ky = repelem(ky,Constants.M+1,Constants.N+1,1,Basis.numOfBasis(2));
periodOfN = -Constants.N/2:Constants.N/2;
ky2 = 2*pi*periodOfN/Dimensions.dy;
ky2 = reshape(ky2,[1 length(periodOfN) 1 1]);
ky2 = repelem(ky2,Constants.M+1,1,length(w),Basis.numOfBasis(2));
kyn = ky + ky2;
g = 0.5./sqrt(kxm.^2+kyn.^2-k.^2);
f = Basis.size(1)*(Basis.size(2)/2).*sinc(kxm*Basis.size(1)/2).*(sinc(kyn*Basis.size(2)/4)).^2;
Q = 1:Basis.numOfBasis(2);
Q = reshape(Q,[1 1 1 length(Q)]);
Q = repelem(Q,Constants.M+1,Constants.N+1,length(w),1);
shift = exp(-1i*kyn.*Q.*Basis.size(2)/2);
I2 = reshape(I,1,1,[],length(w));
I2 = permute(I2,[1 2 4 3]);
Esmn = Esmn.*(k.^2-kyn.^2).*g.*f.*shift.*I2;
Esmn = sum(Esmn,4);
%% Caculate Reflection & Transmission
%% RTM
RTM = -1i*w/((Constants.c)^2);
RTM = reshape(RTM,[1 1 length(w)]);
RTM = repelem(RTM,Constants.M+1,Constants.N+1,1);
% k
k = reshape(Incidence.k,[1 1 length(w)]);
k = repelem(k,Constants.M+1,Constants.N+1,1);
% kxm
kx = reshape(Incidence.kx,[1 1 length(w)]);
kx = repelem(kx,Constants.M+1,Constants.N+1,1);
periodOfM = -Constants.M/2:Constants.M/2;
kx2 = 2*pi*periodOfM/Dimensions.dx;
kx2 = reshape(kx2,[length(periodOfM) 1 1 ]);
kx2 = repelem(kx2,1,Constants.N+1,length(w));
kxm = kx + kx2;
% kyn
ky = reshape(Incidence.ky,[1 1 length(w)]);
ky = repelem(ky,Constants.M+1,Constants.N+1,1);
periodOfN = -Constants.N/2:Constants.N/2;
ky2 = 2*pi*periodOfN/Dimensions.dy;
ky2 = reshape(ky2,[1 length(periodOfN) 1]);
ky2 = repelem(ky2,Constants.M+1,1,length(w));
kyn = ky + ky2;
kzmn = sqrt(k.^2-kxm.^2-kyn.^2);
RTM = RTM.*kyn./(kzmn.*(kxm.^2+kyn.^2));
RTM = RTM.*Esmn;
%% RTE
RTE = -1i*kxm.*Esmn*Constants.eps./(kxm.^2+kyn.^2);
%% TTM
%% TTE
%% Plots
% RTM_MAG = sum(abs(RTM).^2,[1 2]);
% RTM_MAG = reshape(sqrt(RTM_MAG),length(w),1);
% RTE_MAG = sum(abs(RTE).^2,[1 2]);
% RTE_MAG = reshape(sqrt(RTE_MAG),length(w),1);
RTM_MAG = abs(RTM(Constants.M/2+1,Constants.N/2+1,:));
RTM_MAG = reshape(RTM_MAG,length(w),1);
RTE_MAG = abs(RTE(Constants.M/2+1,Constants.N/2+1,:));
RTE_MAG = reshape(RTE_MAG,length(w),1);
subplot(2,1,1)
plot(w*1e-9/(2*pi),20*log10(RTM_MAG))
ylabel("R_T_M")
xlabel("f/GHz")
grid on;
subplot(2,1,2)
plot(w*1e-9/(2*pi),20*log10(RTE_MAG))
ylabel("R_T_E")
xlabel("f/GHz")
grid on;

Last edited:

Hi,

what do you mean with wrong results? What do you expect?

BR

Hi,

what do you mean with wrong results? What do you expect?

BR
I compared the results with simulation on FSS and I didn't get the same result and when I change some parameters as the basis function length or the number of floquet harmonics the Z matrix becomes singular.

Status
Not open for further replies.