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

4 次查看（过去 30 天）
RUI YAN
RUI YAN 2020-5-31
回答： Aman Vyas 2020-7-9
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));

请先登录，再回答此问题。

回答（1 个）

Aman Vyas
Aman Vyas 2020-7-9
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!

类别

在 Help Center 和 File Exchange 中查找有关 Analysis, Benchmarking, and Verification 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by