openEMS/matlab/examples/microstrip/directional_coupler.m

262 lines
6.8 KiB
Matlab

function directional_coupler
%
% EXAMPLE / microstrip / directional_coupler
%
% Stacked directional coupler in microstrip technology.
%
% This example demonstrates:
% - simple microstrip geometry
% - S-parameter calculation using the ypar-method
% - display of coupler parameters
% - display of S11 (smith chart)
%
%
% Tested with
% - Matlab 2010b
% - Octave 3.2.4
% - openEMS v0.0.17
%
% (C) 2010 Sebastian Held <sebastian.held@gmx.de>
clear
close all
clc
% sim settings
showStructure = 1;
runSimulation = 1;
for n=1:4
if n > 1, showStructure = 0; end
ports{n} = sim( n, showStructure, runSimulation );
end
postprocess( ports );
function ports = sim( simnr, showStructure, runSimulation )
physical_constants
% setup the simulation
drawingunit = 1e-6; % specify everything in um
Sim_Path = ['tmp' int2str(simnr)];
Sim_CSX = 'tmp.xml';
f_max = 100e6;
lambda = c0/f_max;
% specify the coupler
pcb1.w = 147000;
pcb1.h = 54500;
pcb1.t = 1524;
pcb1.epr = 3;
msl1.w = 135000;
msl1.h = 2800;
pcb2.w = 107000;
pcb2.h = 14000;
pcb2.t = 1524;
pcb2.epr = 3;
msl2.w = 95000;
msl2.h = 4000;
CSX = InitCSX();
% create the mesh
mesh.x = [-pcb1.w/2 pcb1.w/2 -pcb2.w/2 pcb2.w/2 -msl1.w/2 msl1.w/2 -msl2.w/2 msl2.w/2];
mesh.x = [mesh.x linspace(-msl2.w/2,-msl2.w/2+msl2.h, 5) linspace(msl2.w/2,msl2.w/2-msl2.h, 5)];
mesh.y = [-pcb1.h/2 pcb1.h/2 -pcb2.h/2 pcb2.h/2 -msl1.h/2 msl1.h/2 -msl2.h/2 msl2.h/2];
mesh.z = [linspace(0,pcb1.t,5) linspace(pcb1.t,pcb1.t+pcb2.t,5)];
mesh.z = [mesh.z mesh.z(end)+10*(mesh.z(end)-mesh.z(1))]; % add space above pcb
res = lambda/sqrt(max([pcb1.epr,pcb2.epr])) / 20 / drawingunit;
mesh.x = SmoothMeshLines2(mesh.x,res);
mesh.y = SmoothMeshLines2(mesh.y,res);
mesh.z = SmoothMeshLines2(mesh.z,res);
mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add space for PML
CSX = DefineRectGrid( CSX, drawingunit, mesh );
%% create the structure
% microstrip
CSX = AddMetal( CSX, 'PEC' );
start = [-msl1.w/2, -msl1.h/2, pcb1.t];
stop = [ msl1.w/2, msl1.h/2, pcb1.t];
priority = 100; % the geometric priority is set to 100
CSX = AddBox( CSX, 'PEC', priority, start, stop );
% ground plane
CSX = AddMetal( CSX, 'PEC_ground' );
start = [-pcb1.w/2, -pcb1.h/2, 0];
stop = [ pcb1.w/2, pcb1.h/2, 0];
CSX = AddBox( CSX, 'PEC_ground', priority, start, stop );
% substrate 1
start = [-pcb1.w/2, -pcb1.h/2, 0];
stop = [ pcb1.w/2, pcb1.h/2, pcb1.t];
priority = 10;
CSX = AddMaterial( CSX, 'substrate1' );
CSX = SetMaterialProperty( CSX, 'substrate1', 'Epsilon', pcb1.epr );
CSX = AddBox( CSX, 'substrate1', priority, start, stop );
% substrate 2
start = [-pcb2.w/2, -pcb2.h/2, pcb1.t];
stop = [ pcb2.w/2, pcb2.h/2, pcb1.t+pcb2.t];
priority = 10;
CSX = AddMaterial( CSX, 'substrate2' );
CSX = SetMaterialProperty( CSX, 'substrate2', 'Epsilon', pcb2.epr );
CSX = AddBox( CSX, 'substrate2', priority, start, stop );
% stripline
start = [-msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t];
stop = [ msl2.w/2, msl2.h/2, pcb1.t+pcb2.t];
priority = 100;
CSX = AddBox( CSX, 'PEC', priority, start, stop );
% connections
start = [-msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t];
stop = [-msl2.w/2+msl2.h, -pcb2.h/2, pcb1.t+pcb2.t];
priority = 100;
CSX = AddBox( CSX, 'PEC', priority, start, stop );
start = [ msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t];
stop = [ msl2.w/2-msl2.h, -pcb2.h/2, pcb1.t+pcb2.t];
priority = 100;
CSX = AddBox( CSX, 'PEC', priority, start, stop );
%% ports
% this project needs 4 simulations
for n=1:4
portexcite{n} = [];
end
portexcite{simnr} = 'excite';
% port 1: input port
start = [-msl1.w/2, 0, pcb1.t];
stop = [-msl1.w/2, 0, 0];
[CSX ports{1}] = AddCurvePort( CSX, 999, 1, 50, start, stop, portexcite{1} );
% port 2: output port
start = [msl1.w/2, 0, pcb1.t];
stop = [msl1.w/2, 0, 0];
[CSX ports{2}] = AddCurvePort( CSX, 999, 2, 50, start, stop, portexcite{2} );
% port 3: coupled port
start = [-msl2.w/2+msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t];
stop = [-msl2.w/2+msl2.h/2, -pcb2.h/2, 0];
[CSX ports{3}] = AddCurvePort( CSX, 999, 3, 50, start, stop, portexcite{3} );
% port 4: isolated port
start = [msl2.w/2-msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t];
stop = [msl2.w/2-msl2.h/2, -pcb2.h/2, 0];
[CSX ports{4}] = AddCurvePort( CSX, 999, 4, 50, start, stop, portexcite{4} );
%% setup FDTD parameters & excitation function
max_timesteps = 50000;
min_decrement = 1e-6;
FDTD = InitFDTD( max_timesteps, min_decrement );
FDTD = SetGaussExcite( FDTD, 0, f_max );
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % faster
FDTD = SetBoundaryCond( FDTD, BC );
%% Write openEMS compatible xml-file
if runSimulation
[dummy,dummy,dummy] = rmdir(Sim_Path,'s');
end
[dummy,dummy,dummy] = mkdir(Sim_Path);
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
if showStructure
CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
end
%% run openEMS
openEMS_opts = '';
openEMS_opts = [openEMS_opts ' --engine=fastest'];
% openEMS_opts = [openEMS_opts ' --debug-material'];
% openEMS_opts = [openEMS_opts ' --debug-boxes'];
if runSimulation
RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts );
end
function postprocess( ports )
f = linspace( 0, 100e6, 201 );
Y = calc_ypar( f, ports{1}, 'tmp' );
R = 50;
S = y2s(Y,R);
% insertion loss
IL_dB = -20 * log10(abs(squeeze(S(2,1,:))));
% coupling factor
CF_dB = -20 * log10(abs(squeeze(S(3,1,:))));
% isolation
I_dB = -20 * log10(abs(squeeze(S(4,1,:))));
% directivity
D_dB = -20 * log10(abs(squeeze(S(4,1,:) ./ S(3,1,:))));
figure
plot( f, [IL_dB CF_dB I_dB D_dB] );
legend( {'insertion loss','coupling factor','isolation','directivity'} );
title( ['performance of the coupler for a termination resistance of R=' num2str(R)] );
grid on
smithchart
S11 = squeeze(S(1,1,:));
plot( real(S11), imag(S11) );
legend( 'S_{11}' );
title( ['performance of the coupler for a termination resistance of R=' num2str(R)] );
axis( [-1 1 -1 1] );
function smithchart
% smith chart
figure
if exist( 'smith', 'file' )
% smith chart
% www.ece.rutgers.edu/~orfanidi/ewa
% or cmt toolbox from git.ate.uni-duisburg.de
smith
else
% poor man smith chart
color = [.6 .6 .6];
h = plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', color );
hg = hggroup;
set( h,'Parent',hg );
hold on
plot( hg, 0.25+0.75*sin(0:0.01:2*pi), 0.75*cos(0:0.01:2*pi), 'Color', color );
plot( hg, 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', color );
plot( hg, 0.75+0.25*sin(0:0.01:2*pi), 0.25*cos(0:0.01:2*pi), 'Color', color );
plot( hg, [-1 1], [0 0], 'Color', color );
axis equal
axis off
end
function s = y2s(y, ZL)
% S = y2s(Y, ZL)
%
% Admittance to Scattering transformation
% for square matrices at multiple frequencies
%
% ZL defaults to 50 Ohm
if nargin < 2
ZL = 50;
end
if size(size(y),2) > 2
nF = size(y,3);
else
nF = 1;
end
I = diag(ones(1, size(y,2)))/ZL;
for i=1:nF
%s(:,:,i) = inv(I+y(:,:,i)) * (I-y(:,:,i));
s(:,:,i) = (I+y(:,:,i)) \ (I-y(:,:,i));
end