new feature: near-field to far-field transform
Sebastian Held
parent
661410cd66
commit
1aae16dc95
|
|
@ -0,0 +1,208 @@
|
|||
function [E_theta,E_phi,Prad,Dmax] = AnalyzeNFFF2( Sim_Path, filenames_E, filenames_H, f, theta, phi, r )
|
||||
% [E_theta,E_phi,Prad,Dmax] = AnalyzeNFFF2( Sim_Path, filenames_E, filenames_H, f, theta, phi, r )
|
||||
%
|
||||
% calculates the farfield via a near field to far field transformation
|
||||
%
|
||||
% input:
|
||||
% Sim_Path: simulation directory
|
||||
% filenames_E: cell array of filenames for the time domain fields on the NFFF contour (6 E-planes; hdf5)
|
||||
% filenames_H: cell array of filenames for the time domain fields on the NFFF contour (6 H-planes; hdf5)
|
||||
% f: frequency (Hz) for far field calculation
|
||||
% theta: (degrees) vector of discrete theta values to calculate the far field for
|
||||
% phi: (degrees) vector of discrete phi values to calculate the far field for
|
||||
% r: (optional) Radius (m) at which the E-fields are calculated (default: 1 m)
|
||||
%
|
||||
% output:
|
||||
% E_theta: E_theta(theta,phi); theta component of the electric field strength at radius r
|
||||
% E_phi: E_phi(theta,phi); phi component of the electric field strength at radius r
|
||||
% Prad: time averaged radiated power
|
||||
% Dmax: maximum directivity
|
||||
%
|
||||
% example:
|
||||
% see examples/NF2FF/infDipol.m
|
||||
%
|
||||
% (C) 2010 Sebastian Held <sebastian.held@gmx.de>
|
||||
|
||||
% check arguments
|
||||
error( nargchk(7,7,nargin) );
|
||||
if ~isscalar(f)
|
||||
error 'Currently only one frequency is supported. Call this function multiple times.'
|
||||
end
|
||||
|
||||
|
||||
|
||||
% read time domain field data and transform into frequency domain
|
||||
for n=1:numel(filenames_E)
|
||||
[Ef{n}, E_mesh{n}] = ReadHDF5Dump( [Sim_Path '/' filenames_E{n}], 'Frequency', f );
|
||||
[Hf{n}, H_mesh{n}] = ReadHDF5Dump( [Sim_Path '/' filenames_H{n}], 'Frequency', f );
|
||||
|
||||
% reshape mesh into row vector
|
||||
mesh{n}.x = reshape( E_mesh{n}.lines{1}, 1, [] );
|
||||
mesh{n}.y = reshape( E_mesh{n}.lines{2}, 1, [] );
|
||||
mesh{n}.z = reshape( E_mesh{n}.lines{3}, 1, [] );
|
||||
end
|
||||
|
||||
% create a normal vector for every plane
|
||||
% FIXME!!! this is dependent upon the order of filenames_*
|
||||
n = {};
|
||||
for a=1:6
|
||||
temp = [(a<=2), ((a>=3)&&(a<=4)), (a>=5)];
|
||||
n{a} = temp - 2*mod(a,2)*temp;
|
||||
end
|
||||
|
||||
|
||||
physical_constants
|
||||
|
||||
k = 2*pi*f/c0;
|
||||
center = [0 0 0];
|
||||
Umax = 0;
|
||||
|
||||
phi_idx = 0;
|
||||
for phi_deg_aufpunkt = phi
|
||||
phi_rad_aufpunkt = phi_deg_aufpunkt/180*pi; % radiant
|
||||
phi_idx = phi_idx + 1;
|
||||
|
||||
theta_idx = 0;
|
||||
for theta_deg_aufpunkt = theta
|
||||
theta_rad_aufpunkt = theta_deg_aufpunkt/180*pi; % radiant
|
||||
theta_idx = theta_idx + 1;
|
||||
|
||||
N_theta = 0;
|
||||
N_phi = 0;
|
||||
L_theta = 0;
|
||||
L_phi = 0;
|
||||
for a=1:6
|
||||
[N_theta_,N_phi_,L_theta_,L_phi_] = process_plane( k, n{a}, center, mesh{a}, Ef{a}.FD.values{1}, Hf{a}.FD.values{1}, theta_rad_aufpunkt, phi_rad_aufpunkt );
|
||||
N_theta = N_theta + N_theta_; N_phi = N_phi + N_phi_;
|
||||
L_theta = L_theta + L_theta_; L_phi = L_phi + L_phi_;
|
||||
end
|
||||
|
||||
% E-fields
|
||||
erg_E_theta = -1i*k*exp(-1i*k*r) / (4*pi*r)*(L_phi+Z0*N_theta);
|
||||
erg_E_phi = 1i*k*exp(-1i*k*r) / (4*pi*r)*(L_theta-Z0*N_phi);
|
||||
|
||||
% output
|
||||
E_theta(theta_idx,phi_idx) = erg_E_theta;
|
||||
E_phi(theta_idx,phi_idx) = erg_E_phi;
|
||||
|
||||
% directivity
|
||||
U = r^2/(2*Z0) * sum(abs([erg_E_theta erg_E_phi]).^2);
|
||||
Umax = max( [Umax U] );
|
||||
end
|
||||
end
|
||||
|
||||
% power
|
||||
Prad = 0;
|
||||
for a=1:6
|
||||
[~,~,~,~,P] = process_plane( k, n{a}, center, mesh{a}, Ef{a}.FD.values{1}, Hf{a}.FD.values{1}, theta_rad_aufpunkt, phi_rad_aufpunkt );
|
||||
Prad = Prad + P;
|
||||
end
|
||||
|
||||
% directivity
|
||||
Dmax = 4*pi*Umax / Prad;
|
||||
|
||||
|
||||
% integrate over one plane
|
||||
function [N_theta,N_phi,L_theta,L_phi,Prad] = process_plane( k, n, center, mesh, E_field, H_field, theta_rad_aufpunkt, phi_rad_aufpunkt )
|
||||
% [N_theta,N_phi,L_theta,L_phi,Prad] = process_plane( k, n, center, mesh, E_field, H_field, theta_rad_aufpunkt, phi_rad_aufpunkt )
|
||||
%
|
||||
% k: wave number
|
||||
% n: normal vector of the plane
|
||||
% center: correction coordinates for the center of the antenna
|
||||
% mesh: mesh info
|
||||
% E_field: E field array ?x?x?x3
|
||||
% H_field: H field array ?x?x?x3
|
||||
|
||||
% speed up
|
||||
sin__theta_rad_aufpunkt = sin(theta_rad_aufpunkt);
|
||||
cos__theta_rad_aufpunkt = cos(theta_rad_aufpunkt);
|
||||
sin__phi_rad_aufpunkt = sin(phi_rad_aufpunkt);
|
||||
cos__phi_rad_aufpunkt = cos(phi_rad_aufpunkt);
|
||||
|
||||
if abs(n(1)) == 1
|
||||
% x-plane
|
||||
x = mesh.x(1);
|
||||
[y z] = ndgrid( mesh.y, mesh.z );
|
||||
coord1 = mesh.y.';
|
||||
coord2 = mesh.z.';
|
||||
Ex = squeeze( E_field(1,:,:,1) );
|
||||
Ey = squeeze( E_field(1,:,:,2) );
|
||||
Ez = squeeze( E_field(1,:,:,3) );
|
||||
Hx = squeeze( H_field(1,:,:,1) );
|
||||
Hy = squeeze( H_field(1,:,:,2) );
|
||||
Hz = squeeze( H_field(1,:,:,3) );
|
||||
elseif abs(n(2)) == 1
|
||||
% y-plane
|
||||
y = mesh.y(1);
|
||||
[x z] = ndgrid( mesh.x, mesh.z );
|
||||
coord1 = mesh.x.';
|
||||
coord2 = mesh.z.';
|
||||
Ex = squeeze( E_field(:,1,:,1) );
|
||||
Ey = squeeze( E_field(:,1,:,2) );
|
||||
Ez = squeeze( E_field(:,1,:,3) );
|
||||
Hx = squeeze( H_field(:,1,:,1) );
|
||||
Hy = squeeze( H_field(:,1,:,2) );
|
||||
Hz = squeeze( H_field(:,1,:,3) );
|
||||
elseif abs(n(3)) == 1
|
||||
% z-plane
|
||||
z = mesh.z(1);
|
||||
[x y] = ndgrid( mesh.x, mesh.y );
|
||||
coord1 = mesh.x.';
|
||||
coord2 = mesh.y.';
|
||||
Ex = squeeze( E_field(:,:,1,1) );
|
||||
Ey = squeeze( E_field(:,:,1,2) );
|
||||
Ez = squeeze( E_field(:,:,1,3) );
|
||||
Hx = squeeze( H_field(:,:,1,1) );
|
||||
Hy = squeeze( H_field(:,:,1,2) );
|
||||
Hz = squeeze( H_field(:,:,1,3) );
|
||||
end
|
||||
|
||||
Jx = n(2) .* Hz - n(3) .* Hy;
|
||||
Jy = n(3) .* Hx - n(1) .* Hz;
|
||||
Jz = n(1) .* Hy - n(2) .* Hx;
|
||||
Mx = -n(2) .* Ez + n(3) .* Ey;
|
||||
My = -n(3) .* Ex + n(1) .* Ez;
|
||||
Mz = -n(1) .* Ey + n(2) .* Ex;
|
||||
r_cos_psi = x*sin__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + y*sin__theta_rad_aufpunkt*sin__phi_rad_aufpunkt + z*cos__theta_rad_aufpunkt;
|
||||
e_fkt = exp( +1i*k*r_cos_psi );
|
||||
N_theta = dbltrapz( ( Jx*cos__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + Jy*cos__theta_rad_aufpunkt*sin__phi_rad_aufpunkt - Jz*sin__theta_rad_aufpunkt) .* e_fkt, coord1, coord2 );
|
||||
N_phi = dbltrapz( (-Jx*sin__phi_rad_aufpunkt + Jy*cos__phi_rad_aufpunkt) .* e_fkt, coord1, coord2 );
|
||||
L_theta = dbltrapz( ( Mx*cos__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + My*cos__theta_rad_aufpunkt*sin__phi_rad_aufpunkt - Mz*sin__theta_rad_aufpunkt) .* e_fkt, coord1, coord2 );
|
||||
L_phi = dbltrapz( (-Mx*sin__phi_rad_aufpunkt + My*cos__phi_rad_aufpunkt) .* e_fkt, coord1, coord2 );
|
||||
|
||||
if nargout > 4
|
||||
% Prad requested
|
||||
|
||||
% this is crap! recode it!
|
||||
EH = zeros(size(Ex));
|
||||
for i1 = 1:numel(coord1)
|
||||
for i2 = 1:numel(coord2)
|
||||
E = [Ex(i1,i2) Ey(i1,i2) Ez(i1,i2)];
|
||||
H = [Hx(i1,i2) Hy(i1,i2) Hz(i1,i2)];
|
||||
EH(i1,i2) = real( dot(cross(E,conj(H)),n) );
|
||||
end
|
||||
end
|
||||
Prad = 0.5 * dbltrapz( EH, coord1, coord2 );
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
function Q = dbltrapz(matrix,a,b)
|
||||
%DBLTRAPZ Trapezoidal numerical integration in two dimensions.
|
||||
% Z = DBLTRAPZ(MATRIX,A,B) computes an approximation of the double integral
|
||||
% of MATRIX via the trapezoidal method (with respect to A and B). A and B must be
|
||||
% column vectors of the same length.
|
||||
% index like this: MATRIX(A,B)
|
||||
|
||||
if nargin < 3, error('MATLAB:dblquad:NotEnoughInputs',...
|
||||
'Requires at least three inputs.'); end
|
||||
if size(a,2) ~= 1, error('column vectors required'); end
|
||||
if size(b,2) ~= 1, error('column vectors required'); end
|
||||
|
||||
temp = zeros(size(b));
|
||||
for i = 1:length(b)
|
||||
temp(i) = trapz( a, matrix(:,i) );
|
||||
end
|
||||
|
||||
Q = trapz( b, temp );
|
||||
|
|
@ -0,0 +1,212 @@
|
|||
function infDipol
|
||||
%
|
||||
% infinitesimal dipole example
|
||||
%
|
||||
|
||||
close all
|
||||
clear
|
||||
clc
|
||||
drawnow
|
||||
|
||||
|
||||
postprocessing_only = 0;
|
||||
|
||||
|
||||
|
||||
|
||||
global g
|
||||
setup
|
||||
|
||||
dipole_orientation = 3; % 1,2,3: x,y,z
|
||||
CSX = createStructure(dipole_orientation);
|
||||
if ~postprocessing_only
|
||||
writeCSX( CSX );
|
||||
CSXGeomPlot( [g.Sim_Path '/' g.Sim_CSX] )
|
||||
run;
|
||||
end
|
||||
postprocess;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
function setup
|
||||
global g
|
||||
physical_constants
|
||||
|
||||
% setup the simulation
|
||||
g.drawingunit = 1e-6; % specify everything in um
|
||||
g.Sim_Path = 'tmp';
|
||||
g.Sim_CSX = 'tmp.xml';
|
||||
|
||||
g.f_max = 1e9;
|
||||
g.lambda = c0/g.f_max;
|
||||
|
||||
% setup geometry values
|
||||
g.dipole_length = g.lambda/50 /g.drawingunit;
|
||||
|
||||
|
||||
|
||||
function CSX = createStructure(dipole_orientation)
|
||||
global g
|
||||
physical_constants
|
||||
|
||||
CSX = InitCSX();
|
||||
|
||||
% create an equidistant mesh
|
||||
mesh.x = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10;
|
||||
mesh.y = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10;
|
||||
mesh.z = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10;
|
||||
mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add space for PML
|
||||
CSX = DefineRectGrid( CSX, g.drawingunit, mesh );
|
||||
|
||||
% excitation
|
||||
ex_vector = [0 0 0];
|
||||
ex_vector(dipole_orientation) = 1;
|
||||
start = ex_vector * -g.dipole_length/2;
|
||||
stop = ex_vector * g.dipole_length/2;
|
||||
CSX = AddExcitation( CSX, 'infDipole', 1, ex_vector );
|
||||
CSX = AddBox( CSX, 'infDipole', 1, start, stop );
|
||||
|
||||
% NFFF contour
|
||||
s1 = [-4.5, -4.5, -4.5] * g.dipole_length/2;
|
||||
s2 = [ 4.5, 4.5, 4.5] * g.dipole_length/2;
|
||||
CSX = AddBox( AddDump(CSX,'Et_xn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_xn', 0, s1, [s1(1) s2(2) s2(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Et_xp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_xp', 0, [s2(1) s1(2) s1(3)], s2 );
|
||||
CSX = AddBox( AddDump(CSX,'Et_yn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_yn', 0, s1, [s2(1) s1(2) s2(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Et_yp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_yp', 0, [s1(1) s2(2) s1(3)], s2 );
|
||||
CSX = AddBox( AddDump(CSX,'Et_zn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_zn', 0, s1, [s2(1) s2(2) s1(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Et_zp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_zp', 0, [s1(1) s1(2) s2(3)], s2 );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_xn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_xn', 0, s1, [s1(1) s2(2) s2(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_xp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_xp', 0, [s2(1) s1(2) s1(3)], s2 );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_yn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_yn', 0, s1, [s2(1) s1(2) s2(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_yp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_yp', 0, [s1(1) s2(2) s1(3)], s2 );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_zn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_zn', 0, s1, [s2(1) s2(2) s1(3)] );
|
||||
CSX = AddBox( AddDump(CSX,'Ht_zp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_zp', 0, [s1(1) s1(2) s2(3)], s2 );
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
function writeCSX(CSX)
|
||||
global g
|
||||
% setup FDTD parameters & excitation function
|
||||
max_timesteps = 2000;
|
||||
min_decrement = 1e-6;
|
||||
FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling',10 );
|
||||
FDTD = SetGaussExcite( FDTD, g.f_max/2, g.f_max/2 );
|
||||
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
|
||||
FDTD = SetBoundaryCond( FDTD, BC );
|
||||
|
||||
% Write openEMS compatible xml-file
|
||||
[~,~,~] = rmdir(g.Sim_Path,'s');
|
||||
[~,~,~] = mkdir(g.Sim_Path);
|
||||
WriteOpenEMS([g.Sim_Path '/' g.Sim_CSX],FDTD,CSX);
|
||||
|
||||
|
||||
|
||||
|
||||
function run
|
||||
global g
|
||||
% define openEMS options and start simulation
|
||||
openEMS_opts = '';
|
||||
openEMS_opts = [openEMS_opts ' --engine=fastest'];
|
||||
RunOpenEMS( g.Sim_Path, g.Sim_CSX, openEMS_opts );
|
||||
|
||||
|
||||
|
||||
|
||||
function postprocess
|
||||
global g
|
||||
|
||||
disp( ' ' );
|
||||
disp( ' ********************************************************** ' );
|
||||
disp( ' ' );
|
||||
|
||||
% NFFF contour
|
||||
filenames_E = {'Et_xn.h5','Et_xp.h5','Et_yn.h5','Et_yp.h5','Et_zn.h5','Et_zp.h5'};
|
||||
filenames_H = {'Ht_xn.h5','Ht_xp.h5','Ht_yn.h5','Ht_yp.h5','Ht_zn.h5','Ht_zp.h5'};
|
||||
|
||||
% calculate the far field at phi=0 degrees and at phi=90 degrees
|
||||
thetaRange = 0:2:359;
|
||||
r = 1; % evaluate fields at radius r
|
||||
disp( 'calculating far field at phi=[0 90] deg...' );
|
||||
[E_far_theta,E_far_phi,Prad,Dmax] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.f_max, thetaRange, [0 90], r );
|
||||
|
||||
|
||||
% display power and directivity
|
||||
disp( ['radiated power: Prad = ' num2str(Prad)] );
|
||||
disp( ['directivity: Dmax = ' num2str(Dmax)] );
|
||||
|
||||
|
||||
% calculate the e-field magnitude for phi = 0 deg
|
||||
E_phi0_far = zeros(1,numel(thetaRange));
|
||||
for n=1:numel(thetaRange)
|
||||
E_phi0_far(n) = norm( [E_far_theta(n,1) E_far_phi(n,1)] );
|
||||
end
|
||||
|
||||
% display polar plot
|
||||
figure
|
||||
polar( thetaRange/180*pi, E_phi0_far );
|
||||
ylabel( 'theta / deg' );
|
||||
title( ['electrical far field (V/m) @r=' num2str(r) ' m phi=0 deg'] );
|
||||
legend( 'e-field magnitude', 'Location', 'BestOutside' );
|
||||
|
||||
|
||||
% calculate the e-field magnitude for phi = 90 deg
|
||||
E_phi90_far = zeros(1,numel(thetaRange));
|
||||
for n=1:numel(thetaRange)
|
||||
E_phi90_far(n) = norm([E_far_theta(n,2) E_far_phi(n,2)]);
|
||||
end
|
||||
|
||||
% display polar plot
|
||||
figure
|
||||
polar( thetaRange/180*pi, E_phi90_far );
|
||||
ylabel( 'theta / deg' );
|
||||
title( ['electrical far field (V/m) @r=' num2str(r) ' m phi=90 deg'] );
|
||||
legend( 'e-field magnitude', 'Location', 'BestOutside' );
|
||||
|
||||
|
||||
|
||||
|
||||
% calculate the far field at theta=90 degrees
|
||||
phiRange = 0:2:359;
|
||||
r = 1; % evaluate fields at radius r
|
||||
disp( 'calculating far field at theta=90 deg...' );
|
||||
[E_far_theta,E_far_phi] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.f_max, 90, phiRange, r );
|
||||
|
||||
E_theta90_far = zeros(1,numel(phiRange));
|
||||
for n=1:numel(phiRange)
|
||||
E_theta90_far(n) = norm([E_far_theta(1,n) E_far_phi(1,n)]);
|
||||
end
|
||||
|
||||
% display polar plot
|
||||
figure
|
||||
polar( phiRange/180*pi, E_theta90_far );
|
||||
ylabel( 'phi / deg' );
|
||||
title( ['electrical far field (V/m) @r=' num2str(r) ' m theta=90 deg'] );
|
||||
legend( 'e-field magnitude', 'Location', 'BestOutside' );
|
||||
|
||||
|
||||
|
||||
|
||||
% calculate 3D pattern
|
||||
phiRange = 0:15:360;
|
||||
thetaRange = 0:10:180;
|
||||
r = 1; % evaluate fields at radius r
|
||||
disp( 'calculating 3D far field...' );
|
||||
[E_far_theta,E_far_phi] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.f_max, thetaRange, phiRange, r );
|
||||
E_far = sqrt( abs(E_far_theta).^2 + abs(E_far_phi).^2 );
|
||||
E_far_normalized = E_far / max(E_far(:));
|
||||
[theta,phi] = ndgrid(thetaRange/180*pi,phiRange/180*pi);
|
||||
x = E_far_normalized .* sin(theta) .* cos(phi);
|
||||
y = E_far_normalized .* sin(theta) .* sin(phi);
|
||||
z = E_far_normalized .* cos(theta);
|
||||
figure
|
||||
surf( x,y,z, E_far_normalized );
|
||||
axis equal
|
||||
xlabel( 'x' );
|
||||
ylabel( 'y' );
|
||||
zlabel( 'z' );
|
||||
Loading…
Reference in New Issue