How do i get the radiation pattern polar plot for the E_side,E_t​op,H_side,​H_top?(cod​e related)

2 views (last 30 days)
RUI YAN
RUI YAN on 31 May 2020
Answered: Aman Vyas on 9 Jul 2020
The code is below and description
% FD3D.m
% Sullivan's FDTD method in 3D
% t. kaufmann, 10. 09. 2011
% Initialise
clear all
clc
% Setup
% Antenna Frequency
config.freq=2e9;
% Set Computational domain size (in terms of lambda, smaller is faster)
config.sizexlambda=4;
config.sizeylambda=4;
config.sizezlambda=3;
% Set resolution (should be larger than lambda/10, coarser resolution is
% faster)
config.resolution=1/12;
% Set Simulation time in periods of the simulation frequency
config.Nperiods=15;
% Compute necessary parameters
config.lambda=3e8/config.freq;
config.ddx=config.resolution*config.lambda;
config.sizex=config.sizexlambda*config.lambda;
config.sizey=config.sizeylambda*config.lambda;
config.sizez=config.sizezlambda*config.lambda;
config.Nx=ceil(config.sizex/config.ddx);
config.Ny=ceil(config.sizey/config.ddx);
config.Nz=ceil(config.sizez/config.ddx);
config.sourcepos=ceil([config.Nx/2 config.Ny/2 1]);
config.nsource=ceil(config.lambda/4/config.ddx);
config.dt=0.9*config.ddx/6e8;
config.R0=3/config.freq/config.dt;
config.NSTEPS=config.Nperiods/(config.dt*config.freq);
epsz=8.8e-12;
config.NPML=7;
config.sigma0=1;
% Initialize fields and material parameters
[fields, materials] = FD3D_initialize(config);
% Start time iteration
for tt=1:config.NSTEPS
% update electric and magnetic fields
fields=FD3D_update_fields(fields, materials, config, tt);
% extract two cut planes: top and side view. The function returns the
% MAGNITUDES of the electric and magnetic field vectors.
[E_top E_side H_top H_side]=FD3D_extract_cuts(fields, config);
% plot the electric fields 100 times (not at every time step)
if (mod(tt, config.NSTEPS/100)<1)
disp([int2str(tt/config.NSTEPS*100), '% ', int2str(tt), '/', int2str(config.NSTEPS), ' (',num2str(tt*config.dt*config.freq),' periods)'])
subplot(2,1,1)
pcolor(E_top);
axis equal;
shading flat;
caxis([0 1])
colorbar
title('|E|, top view');
xlabel('x'); ylabel('y');
subplot(2,1,2)
E_side(E_side>1)=NaN;
pcolor(E_side');
axis equal;
shading flat;
caxis([0 1])
colorbar
title('|E|, side view');
xlabel('y'); ylabel('z');
pause(0.05)
end
end
% Implementation of update equations for 3D FDTD code
%
% fields = FD3D_update_fields(fields, materials, config, tt)
% update the electromagnetic fields based on material properties,
% configuration and time step tt
%
% t. kaufmann, 10. 09. 2011
function fields = FD3D_update_fields(fields, materials, config, tt)
curl_h=(fields.Hz(2:end, 2:end, 2:end)-fields.Hz(2:end, 1:end-1, 2:end) - ...
fields.Hy(2:end, 2:end, 2:end) + fields.Hy(2:end, 2:end, 1:end-1));
fields.Iexl(2:end, 2:end, 2:end)=fields.Iexl(2:end, 2:end, 2:end)+curl_h;
fields.Ex(2:end, 2:end, 2:end)= ...
materials.gj3(2:end, 2:end, 2:end).*materials.gk3(2:end, 2:end, 2:end).*fields.Ex(2:end, 2:end, 2:end)+...
materials.gj2(2:end, 2:end, 2:end).*materials.gk2(2:end, 2:end, 2:end).*0.5.*...
(curl_h+materials.gi1(2:end, 2:end, 2:end).*fields.Iexl(2:end, 2:end, 2:end));
curl_h=(fields.Hx(2:end, 2:end, 2:end)-fields.Hx(2:end, 2:end, 1:end-1) - ...
fields.Hz(2:end, 2:end, 2:end) + fields.Hz(1:end-1, 2:end, 2:end));
fields.Ieyl(2:end, 2:end, 2:end)=fields.Ieyl(2:end, 2:end, 2:end)+curl_h;
fields.Ey(2:end, 2:end, 2:end)= ...
materials.gi3(2:end, 2:end, 2:end).*materials.gk3(2:end, 2:end, 2:end).*fields.Ey(2:end, 2:end, 2:end)+...
materials.gi2(2:end, 2:end, 2:end).*materials.gk2(2:end, 2:end, 2:end).*0.5.*...
(curl_h+materials.gj1(2:end, 2:end, 2:end).*fields.Ieyl(2:end, 2:end, 2:end));
curl_h=(fields.Hy(2:end, 2:end, 2:end)-fields.Hy(1:end-1, 2:end, 2:end) - ...
fields.Hx(2:end, 2:end, 2:end) + fields.Hx(2:end, 1:end-1, 2:end));
fields.Iezl(2:end, 2:end, 2:end)=fields.Iezl(2:end, 2:end, 2:end)+curl_h;
fields.Ez(2:end, 2:end, 2:end)=...
materials.gi3(2:end, 2:end, 2:end).*materials.gj3(2:end, 2:end, 2:end).*fields.Ez(2:end, 2:end, 2:end)+...
materials.gi2(2:end, 2:end, 2:end).*materials.gj2(2:end, 2:end, 2:end).*0.5.*...
(curl_h+materials.gk1(2:end, 2:end, 2:end).*fields.Iezl(2:end, 2:end, 2:end));
% xy-plane
fields.Ex(1:end, 1:end, [1 end])=0;
fields.Ey(1:end, 1:end, [1 end])=0;
fields.Ez(1:end, 1:end, end)=0;
fields.Ez(2:end,2:end,1)=fields.Ez(2:end,2:end,1)+0.5*(fields.Hy(2:end, 2:end, 1)-fields.Hy(1:end-1, 2:end, 1) - ...
fields.Hx(2:end, 2:end, 1) + fields.Hx(2:end, 1:end-1, 1));
% xz-plane
fields.Ex(1:end, [1 end], 1:end)=0;
fields.Ey(1:end, [1 end], 1:end)=0;
fields.Ez(1:end, [1 end], 1:end)=0;
% yz-plane
fields.Ex([1 end], 1:end, 1:end)=0;
fields.Ey([1 end], 1:end, 1:end)=0;
fields.Ez([1 end], 1:end, 1:end)=0;
%kz=ones(1,1,config.nsource); kz(:)=((0:config.nsource-1)/config.nsource);
fields.Ez(config.sourcepos(1), config.sourcepos(2), 3:config.nsource)=0;
fields.Ez(config.sourcepos(1), config.sourcepos(2), 1:2)= ...
5*sin(2*pi*config.freq*config.dt*tt)*((tt<=config.R0).*(tt/config.R0)+(tt>config.R0));
curl_e=(fields.Ey(1:end-1, 1:end-1, 2:end)-fields.Ey(1:end-1, 1:end-1, 1:end-1) - ...
fields.Ez(1:end-1, 2:end, 1:end-1) + fields.Ez(1:end-1, 1:end-1, 1:end-1));
fields.Ihxl(1:end-1, 1:end-1, 1:end-1)=fields.Ihxl(1:end-1, 1:end-1, 1:end-1)+curl_e;
fields.Hx(1:end-1, 1:end-1, 1:end-1)= ...
materials.fj3(1:end-1, 1:end-1, 1:end-1).*materials.fk3(1:end-1, 1:end-1, 1:end-1).*fields.Hx(1:end-1, 1:end-1, 1:end-1)+...
materials.fj2(1:end-1, 1:end-1, 1:end-1).*materials.fk2(1:end-1, 1:end-1, 1:end-1).*0.5.*...
(curl_e+materials.fi1(1:end-1, 1:end-1, 1:end-1).*fields.Ihxl(1:end-1, 1:end-1, 1:end-1));
curl_e=(fields.Ez(2:end, 1:end-1, 1:end-1)-fields.Ez(1:end-1, 1:end-1, 1:end-1) - ...
fields.Ex(1:end-1, 1:end-1, 2:end) + fields.Ex(1:end-1, 1:end-1, 1:end-1));
fields.Ihyl(1:end-1, 1:end-1, 1:end-1)=fields.Ihyl(1:end-1, 1:end-1, 1:end-1)+curl_e;
fields.Hy(1:end-1, 1:end-1, 1:end-1)= ...
materials.fi3(1:end-1, 1:end-1, 1:end-1).*materials.fk3(1:end-1, 1:end-1, 1:end-1).*fields.Hy(1:end-1, 1:end-1, 1:end-1)+...
materials.fi2(1:end-1, 1:end-1, 1:end-1).*materials.fk2(1:end-1, 1:end-1, 1:end-1).*0.5.*... %!!!!!!!!!!!!! fk2
(curl_e+materials.fj1(1:end-1, 1:end-1, 1:end-1).*fields.Ihyl(1:end-1, 1:end-1, 1:end-1));
curl_e=(fields.Ex(1:end-1, 2:end, 1:end-1)-fields.Ex(1:end-1, 1:end-1, 1:end-1) - ...
fields.Ey(2:end, 1:end-1, 1:end-1) + fields.Ey(1:end-1, 1:end-1, 1:end-1));
fields.Ihzl(1:end-1, 1:end-1, 1:end-1)=fields.Ihzl(1:end-1, 1:end-1, 1:end-1)+curl_e;
fields.Hz(1:end-1, 1:end-1, 1:end-1)= ...
materials.fi3(1:end-1, 1:end-1, 1:end-1).*materials.fj3(1:end-1, 1:end-1, 1:end-1).*fields.Hz(1:end-1, 1:end-1, 1:end-1)+...
materials.fi2(1:end-1, 1:end-1, 1:end-1).*materials.fj2(1:end-1, 1:end-1, 1:end-1).*0.5.*...
(curl_e+materials.fk1(1:end-1, 1:end-1, 1:end-1).*fields.Ihzl(1:end-1, 1:end-1, 1:end-1));
% BC
% xy-plane
fields.Hz(1:end, 1:end, [1 end])=0;
% xz-plane
fields.Hy(1:end, [1 end], 1:end)=0;
% yz-plane
%fields.Hz([1 end], 1:end, 1:end)=0;
end
% Initialization procedure for 3D FDTD code
%
% [fields materials]=FD3D_initialize(config)
% initialize electromagnetic fields and material parameters based on
% configuration setting.
%
% t. kaufmann, 10. 09. 2011
function [fields materials]=FD3D_initialize(config)
fields.Ex=zeros(config.Nx, config.Ny, config.Nz);
fields.Ey=zeros(config.Nx, config.Ny, config.Nz);
fields.Ez=zeros(config.Nx, config.Ny, config.Nz);
fields.Hx=zeros(config.Nx, config.Ny, config.Nz);
fields.Hy=zeros(config.Nx, config.Ny, config.Nz);
fields.Hz=zeros(config.Nx, config.Ny, config.Nz);
fields.Iexl=zeros(config.Nx, config.Ny, config.Nz);
fields.Ieyl=zeros(config.Nx, config.Ny, config.Nz);
fields.Iezl=zeros(config.Nx, config.Ny, config.Nz);
fields.Ihxl=zeros(config.Nx, config.Ny, config.Nz);
fields.Ihyl=zeros(config.Nx, config.Ny, config.Nz);
fields.Ihzl=zeros(config.Nx, config.Ny, config.Nz);
materials.fi1=zeros(config.Nx, config.Ny, config.Nz); materials.fj1=zeros(config.Nx, config.Ny, config.Nz); materials.fk1=zeros(config.Nx, config.Ny, config.Nz);
materials.fi2=ones(config.Nx, config.Ny, config.Nz); materials.fj2=ones(config.Nx, config.Ny, config.Nz); materials.fk2=ones(config.Nx, config.Ny, config.Nz);
materials.fi3=ones(config.Nx, config.Ny, config.Nz); materials.fj3=ones(config.Nx, config.Ny, config.Nz); materials.fk3=ones(config.Nx, config.Ny, config.Nz);
materials.gi1=zeros(config.Nx, config.Ny, config.Nz); materials.gj1=zeros(config.Nx, config.Ny, config.Nz); materials.gk1=zeros(config.Nx, config.Ny, config.Nz);
materials.gi2=ones(config.Nx, config.Ny, config.Nz); materials.gj2=ones(config.Nx, config.Ny, config.Nz); materials.gk2=ones(config.Nx, config.Ny, config.Nz);
materials.gi3=ones(config.Nx, config.Ny, config.Nz); materials.gj3=ones(config.Nx, config.Ny, config.Nz); materials.gk3=ones(config.Nx, config.Ny, config.Nz);
%pml in x direction
xxn=(config.NPML:-1:1)/config.NPML;
xn=config.sigma0*xxn.^3/3;
materials.fi1(1:config.NPML,:,:)=repmat(xn', [1, config.Ny, config.Nz]);
materials.fi1(config.Nx-config.NPML+1:end,:,:)=repmat(xn(end:-1:1)', [1, config.Ny, config.Nz]);
materials.gi2(1:config.NPML,:,:)=repmat((1./(1+xn))', [1, config.Ny, config.Nz]);
materials.gi2(config.Nx-config.NPML+1:end,:,:)=repmat((1./(1+xn(end:-1:1)))', [1, config.Ny, config.Nz]);
materials.gi3(1:config.NPML,:,:)=repmat(((1-xn)./(1+xn))', [1, config.Ny, config.Nz]);
materials.gi3(config.Nx-config.NPML+1:end,:,:)=repmat(((1-xn(end:-1:1))./(1+xn(end:-1:1)))', [1, config.Ny, config.Nz]);
xxn=((config.NPML:-1:1)-0.5)/config.NPML;
xn=config.sigma0*xxn.^3/3;
materials.gi1(1:config.NPML,:,:)=repmat(xn', [1, config.Ny, config.Nz]);
materials.gi1(config.Nx-config.NPML:end-1,:,:)=repmat(xn(end:-1:1)', [1, config.Ny, config.Nz]);
materials.fi2(1:config.NPML,:,:)=repmat((1./(1+xn))', [1, config.Ny, config.Nz]);
materials.fi2(config.Nx-config.NPML:end-1,:,:)=repmat((1./(1+xn(end:-1:1)))', [1, config.Ny, config.Nz]);
materials.fi3(1:config.NPML,:,:)=repmat(((1-xn)./(1+xn))', [1, config.Ny, config.Nz]);
materials.fi3(config.Nx-config.NPML:end-1,:,:)=repmat(((1-xn(end:-1:1))./(1+xn(end:-1:1)))', [1, config.Ny, config.Nz]);
%pml in y direction
xxn=(config.NPML:-1:1)/config.NPML;
xn=config.sigma0*xxn.^3/3;
materials.fj1(:, 1:config.NPML,:)=repmat(xn, [config.Nx, 1, config.Nz]);
materials.fj1(:, config.Ny-config.NPML+1:end,:)=repmat(xn(end:-1:1), [config.Nx, 1, config.Nz]);
materials.gj2(:, 1:config.NPML,:)=repmat((1./(1+xn)), [config.Nx, 1, config.Nz]);
materials.gj2(:, config.Ny-config.NPML+1:end,:)=repmat((1./(1+xn(end:-1:1))), [config.Nx, 1, config.Nz]);
materials.gj3(:, 1:config.NPML,:)=repmat(((1-xn)./(1+xn)), [config.Nx, 1, config.Nz]);
materials.gj3(:, config.Ny-config.NPML+1:end,:)=repmat(((1-xn(end:-1:1))./(1+xn(end:-1:1))), [config.Nx, 1, config.Nz]);
xxn=((config.NPML:-1:1)-0.5)/config.NPML;
xn=config.sigma0*xxn.^3/3;
materials.gj1(:, 1:config.NPML,:)=repmat(xn, [config.Nx, 1, config.Nz]);
materials.gj1(:, config.Ny-config.NPML:end-1,:)=repmat(xn(end:-1:1), [config.Nx, 1, config.Nz]);
materials.fj2(:, 1:config.NPML,:)=repmat((1./(1+xn)), [config.Nx, 1, config.Nz]);
materials.fj2(:, config.Ny-config.NPML:end-1,:)=repmat((1./(1+xn(end:-1:1))), [config.Nx, 1, config.Nz]);
materials.fj3(:, 1:config.NPML,:)=repmat(((1-xn)./(1+xn)), [config.Nx, 1, config.Nz]);
materials.fj3(:, config.Ny-config.NPML:end-1,:)=repmat(((1-xn(end:-1:1))./(1+xn(end:-1:1))), [config.Nx, 1, config.Nz]);
%pml in z direction
xxn=(config.NPML:-1:1)/config.NPML;
xn=zeros(1,1,config.NPML);
xn(1,1,:)=config.sigma0*xxn.^3/3;
%fk1(:, :, 1:config.NPML)=repmat(xn, [config.Nx, config.Ny, 1]);
materials.fk1(:, :, end:-1:end-config.NPML+1)=repmat(xn, [config.Nx, config.Ny, 1]);
%gk2(:, :, 1:config.NPML)=repmat((1./(1+xn)), [config.Nx, config.Ny, 1]);
materials.gk2(:, :, end:-1:end-config.NPML+1)=repmat((1./(1+xn)), [config.Nx, config.Ny, 1]);
%gk3(:, :, 1:config.NPML)=repmat(((1-xn)./(1+xn)), [config.Nx, config.Ny, 1]);
materials.gk3(:, :, end:-1:end-config.NPML+1)=repmat(((1-xn)./(1+xn)), [config.Nx, config.Ny, 1]);
xxn=((config.NPML:-1:1)-0.5)/config.NPML;
xn=zeros(1,1,config.NPML);
xn(1,1,:)=config.sigma0*xxn.^3/3;
%gk1(:, :, 1:config.NPML)=repmat(xn, [config.Nx, config.Ny, 1]);
materials.gk1(:, :, end-1:-1:end-config.NPML)=repmat(xn, [config.Nx, config.Ny, 1]);
%fk2(:, :, 1:config.NPML)=repmat((1./(1+xn)), [config.Nx, config.Ny, 1]);
materials.fk2(:, :, end-1:-1:end-config.NPML)=repmat((1./(1+xn)), [config.Nx, config.Ny, 1]);
%fk3(:, :, 1:config.NPML)=repmat(((1-xn)./(1+xn)), [config.Nx, config.Ny, 1]);
materials.fk3(:, :, end-1:-1:end-config.NPML)=repmat(((1-xn)./(1+xn)), [config.Nx, config.Ny, 1]);
% Extraction of magnetic and electric field values on cut planes
%
% [E_top E_side H_top H_side]=FD3D_extract_cuts(fields, config)
% take the 3D matrices fields.Ex, fields.Ey, fields.Ez and
% fields.Hx, fields.Hy, fields.Hz, and extract the values at the cut planes
% top view (z=1) and side view (x=1). Return the length of the vectors.
%
% t. kaufmann, 10. 09. 2011
function [E_top E_side H_top H_side]=FD3D_extract_cuts(fields, config)
E_top=sqrt(squeeze(fields.Ex(:,:,1).^2+fields.Ey(:,:,1).^2+fields.Ez(:,:,1).^2));
E_side=sqrt(squeeze(fields.Ex(config.sourcepos(1),:,:).^2+fields.Ey(config.sourcepos(1),:,:).^2+fields.Ez(config.sourcepos(1),:,:).^2));
H_top=sqrt(squeeze(fields.Hx(:,:,1).^2+fields.Hy(:,:,1).^2+fields.Hz(:,:,1).^2));
H_side=sqrt(squeeze(fields.Hx(config.sourcepos(1),:,:).^2+fields.Hy(config.sourcepos(1),:,:).^2+fields.Hz(config.sourcepos(1),:,:).^2));

Sign in to answer this question.

Answers (1)

Aman Vyas
Aman Vyas on 9 Jul 2020
Hi,
You can try patternCustom , fieldCustom functions to visualise the data in whichever format you like.
You might need to define wavelength to set the values for the 3 dipoles.
For more info,you can refer the documentation pages on radiation pattern polar plot.
Hope this helps!

Categories

Find more on Polar Plots in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!