Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Register Log in

Solving frequency selective surfaces using method of moments on matlab

Status
Not open for further replies.
A

ahmedatef0

Junior Member level 1
Joined
Oct 16, 2021
Messages
15
Helped
0
Reputation
0
Reaction score
1
Trophy points
1
Activity points
151
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
 

stenzer said:
Hi,

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

BR
Click to expand...
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.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top