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

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));

Answer 1

[Aman Vyas]

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!