How do i get the radiation pattern polar plot for the E_side,E_top,H_side,H_top?(code related)
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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));
0 comentarios
Respuestas (1)
Aman Vyas
el 9 de Jul. de 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.
Hope this helps!
0 comentarios
Ver también
Categorías
Más información sobre Polar Plots en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!