ZIP with MATLAB scripts and note:
exercise 6.14
exercise 6.14 notes:
Support function J_bessel_zeros supplied that calculates all zeros of 1st kind Bessel functions and zeros of its 1st derivative within input interval.
The resolution of the input interval is proportional to the accuracy of the roots
The solutions manual gets a from resonant frequency expression for cylindrical cavity applying n=1;m=1;L=1 and obtaining a=0.02107
Then uses the peak of the relevant mode TE_111 in the normalised unloaded Q graph above to obtain d .
Double check:
load bessel_roots.mat
p1_11=p1{2}(1)
c0/(2*pi*er^.5)*((p1_11/.02107)^2+(pi/(2*.02107/1.7))^2)^.5
= 5.977795786452106e+09
Run M=5;x=[0:.001:15]; [p p1]=J_Bessel_zeros(M,x) to calculate cells p and p1 stored in bessel_roots.mat . bessel_roots.mat and J_Bessel_zeros are included in ZIP download from Dropbox link right corner of this page.
To find radius a one can also use solve
aovd=1.7 % aovd = a/d . [COL] graph is Qc(2*a/d)
syms a_sym
eq1=fc_TE_111==c0/(2*pi*(mur*er)^.5)*((p1_11/a_sym)^2+(aovd*L*pi/(2*a_sym))^2)^.5
a1=double(solve(eq1,a_sym))
a_radius=a1(2)
a1 =
-0.020992026203424
0.020992026203424
a1 = 0.020992026203424
I tried 3 times to obtain the normalised unloaded Q graphs from [POZAR] Qc expression, but no way.
Either I did not apply the correct expression, what at this point seems to be the cause, [COL] and [POZAR] expressions for unloaded Qc(2*d/a) are different.
Since Q(a/d) is just a slice of Q(a,d,n,m,L) I generated surfaces Q(a,d) but they all show continuous downhill slopes, no peaks.
[POZAR] also briefly mentions just underneath figure 6.10 that assumes k*a and beta both frequency independent.
But both k*a and beta=L*pi/d are frequency dependent.
Also, using MATLAB instead of simple visual graph peak eye spotting helps avoid 2*aovd assumed to be 1.7 , that is not this value , as you will see after reading the exact value out of calculated [COL] graph. In any case exact peaks detection should be the choice to avoid errors propagation.
clear all;close all;clc
c0 = 298792586 % m/s
mu0 = 1.256637061435917e-06
eps0 = 8.854187817e-12 % [F/m]=[C^2/N*m^2] permittivity in empty space
etha0 = (mu0/eps0)^.5 % ~ 120*pi [ohm] characteristic impedance of empty space
er=1.5 % relative permittivity of filler
mur=1
etha=etha0/er^.5
f0=6e9
lambda0=c0/f0;
k0=2*pi/lambda0
k1=2*pi*f0*er^.5/c0
% same as
k1=k0*er^.5
n=1 % mode TE_111 only
m=1
L=1
tand=.0005 % dielectric attenuation per metre
% metal attenuation
sigma_Au=4.1e7 % S/m
Rs=(2*pi*f0*mu0/(2*sigma_Au))^.5
lambda=lambda0/er^.5
a3=1 % [POZAR] writes Collin assumes k*a frequency independent yet there's an a^3 hanging term left undefined
aovd=[0:.001:3]; % this is a/d
beta=L*pi*aovd % beta=L*pi/d=L*pi*a/d*1/a ~ L*pi*aovd/a3=L*pi*aovd
% plot Qc(a/d) for different TE TM modes
% calculating p and p1, J and J1 roots respectively.
M=5;x=[0:.001:15];
% [p p1]=J_Bessel_zeros(M,x)
% save bessel_roots.mat p p1
load bessel_roots.mat
Since MATLAB doesn't use null matrix index, yet Bessel zeros p and p1 are defined using null index, to retrieve particular zeros p_nm or p1_nm
out of cells p and p1, where p_nm means Bessel 1st kind order n-th, m-th zero, one has to add one (+1) to n input to cells p. Retrieve zeros using the following format: p{n+1}(m) or p1{n+1}(m).
For instance to get p1_11 needed in this exercise, one has to do the following:
p1_11=p1{2}(1)
% p1{1}(1) is p1_11 with L=2 TE_012
% p1{1}(1) % TE_011
% p1{2}(1) % TE_111
% p{2}(1) % TM_111
% p{1}(1) % TM_010
z1=[p1{1}(1) p1{1}(1) p1{2}(1) p{2}(1) p{1}(1)]
L1=[2 1 1 1 0]
n1=[0 0 1 1 0]
m1=[1 1 1 1 1]
EM=[1 1 0 1 0] % 1 for Electric, 0 for Magnetic
s1=1/pi*k1^3./(4*z1.^2)
s2=1-(n1./z1).^2
s3=.5*pi^2*n1.^2./z1.^4.*L1.^2
s4=L1.^2.*pi^2.*1./z1.^2.*(1-n1.^2./z1.^2)
Let's check what comes up of [POZAR] expression, that doesn't really look like [COL] expression for same parameter:
Q_norm_graph=zeros(1,numel(aovd));
cell_legend={}
for r=1:1:numel(L1)
Qc_norm1=s2(r).*a3./(1./aovd.*...
(.5+s3(r).*L1(r).^2.*aovd.^2)+...
L1(r).^2.*s4(r).*aovd.^2);
Q_norm_graph(r,:)=Qc_norm1;
switch EM(r)
case 1
cell_legend{r}=...
case 0
cell_legend{r}=['TM' num2str(n1(r))
num2str(m1(r)) num2str(L1(r))]
otherwise
end
end
h10=figure;
hold all
grid on
for r=1:1:numel(L1)
plot(aovd,10*log10(Q_norm_graph(r,:)))
end
legend(cell_legend)
title('Qc(a/b) in dB')
h11=figure;
hold all
grid on
for r=1:1:numel(L1)
plot(aovd,Q_norm_graph(r,:))
end
legend(cell_legend)
title('Qc(a/b) linear')
for r=1:1:numel(L1)
Qc_norm1=s2(r).*a3./...
((.5+s3(r).*L1(r).^2.*aovd.^2)+...
L1(r).^2.*s4(r).*aovd.^2);
Q_norm_graph(r,:)=Qc_norm1;
end
h12=figure;
hold all
grid on
for r=1:1:numel(L1)
plot(aovd,10*log10(Q_norm_graph(r,:)))
end
legend(cell_legend)
title('Qc(a/b) in dB')
h14=figure;
hold all
grid on
for r=1:1:numel(L1)
plot(aovd,Q_norm_graph(r,:))
end
legend(cell_legend)
title('Qc(a/b) linear')
These are not the expected functions that [POZAR] shows in page 292.
Attempting to reproduce [POZAR] normalised Qc as Qc(a,d) with surface
a_max=.075 % max radius 2.5cm
da=a_max/100
a=[0:da:a_max]; % a range 0 to 2.5cm
d_max=.75 % d max 1/4 metre
dd=d_max/100
d=[0:dd:d_max];
[A,D]=meshgrid(a,d);
Qc_norm_surf=zeros([size(A,1) size(A,2) size(L1,2)]);
for r=1:1:numel(L1)
beta=L1(r)*pi./D
Qc_norm2=...
1/pi*(k1*A).^3*A.*D./(4*z1(r)^2).*(1-(n1(r)/z1(r))^2)./...
(A.*D/2.*(1+(beta.*A*n1(r)/(z1(r))^2).^2)+...
(beta*A.^2/z1(r)).^2*(1-n1(r)^2/z1(r)^2))
Qc_norm_surf(:,:,r)=Qc_norm2;
end
figure;hs1=surf(Qc_norm_surf(:,:,2))
hs1.EdgeColor='none'
% fixing beta
for r=1:1:numel(L1)
beta=L1(r)*pi./D
Qc_norm2=...
1/pi*(k1*A).^3*A.*D./(4*z1(r)^2).*(1-(n1(r)/z1(r))^2)./...
(A.*D/2.*(1+(beta.*A*n1(r)/(z1(r))^2).^2)+...
(beta*A.^2/z1(r)).^2*(1-n1(r)^2/z1(r)^2))
Qc_norm_surf(:,:,r)=Qc_norm2;
end
figure;hs1=surf(Qc_norm_surf(:,:,2))
hs1.EdgeColor='none'
All obtained surfaces look similar: continuous downhill slope, no peaks.
Now using [COL] expressions for TE
aovd=[0:.001:3];
z1=[p1{1}(1) p1{1}(1) p1{2}(1)]
L1=[2 1 1]
n1=[0 0 1]
m1=[1 1 1]
EM=[1 1 1] % 1 for Electric, 0 for Magnetic
% TE
Qc_col_TE_nmL=zeros(numel(L1),numel(aovd));
Qc_col_TM_nmL=zeros(numel(L1),numel(aovd));
for r=1:1:numel(L1)
Qc_col_TE=...
(1-(n1(r)./z1(r)).^2).*(z1(r).^2+(L1(r).*pi*aovd).^2).^1.5./...
(2*pi*(z1(r).^2+2*aovd.*(L1(r).*pi*aovd).^2 +...
(1-2*aovd).*(n1(r)*L1(r)*pi*aovd./z1(r)).^2));
Qc_col_TE_nmL(r,:)=Qc_col_TE;
switch EM(r)
case 1
cell_legend{r}=['TE' num2str(n1(r)) num2str(m1(r)) num2str(L1(r))]
case 0
cell_legend{r}=['TM' num2str(n1(r)) num2str(m1(r)) num2str(L1(r))]
otherwise
end
end
h24=figure;
hold all
grid on
for r=1:1:numel(L1)
plot(aovd,Qc_col_TE_nmL(r,:))
end
legend(cell_legend)
title('Qc(a/b) [RCOL]')
If any reader finds a way to connect [POZAR] pg291 expression 6.57 that supposedly produces figure 6.10 in next pg292, and kind enough to let know it would be appreciated.
To exactly locate TE111 ||Qc|| peak in such graph, TE111 is the 3rd row in Qc_col_TE_nmL
Qc_E111=Qc_col_TE_nmL(3,:);
n_peak=find(Qc_E111==max(Qc_E111))
% and the sought a over d ratio is:
aovd(n_peak)
% for the other 2 calculated TE modes the peaks are located in the
% following a/d ratios:
Qc_E012=Qc_col_TE_nmL(1,:);
n_peak_TE012=find(Qc_E012==max(Qc_E012))
aovd(n_peak_TE012)
Qc_E011=Qc_col_TE_nmL(2,:);
n_peak_TE011=find(Qc_E011==max(Qc_E011))
aovd(n_peak_TE011)
n_peak_TE111 = 658
= 0.657000000000000
n_peak_TE012 = 501
= 0.500000000000000
n_peak_TE011 = 501
= 0.500000000000000
So the reason why the solutions manual chooses 2*a/d=1.7 is that such ratio was probably just visual assessment.
Remark: Since I have put ratio a/d in variable aovd, while leaving any factor 2 out of aovd, such factor 2 has to be taken into account when comparing literature Qc graphs and these produced with MATLAB.
% skin depth at f0 for given material
ds=(2/(2*pi*f0*mu0*sigma_Au))^.5
% obtain Qc peak by inverting Qc_norm=Qc/lambda0*ds
Qc_peak=Qc_E111(n_peak_TE111)*lambda0/ds
ds = 1.014734854843133e-06
Qc_peak = 1.356381041390617e+04
Let's compare Qc_peak obtained with [COL] Qc_norm=Qc*ds/lambda0 with an apparently equal normalisation briefly mentioned by [POZAR] that return same value:
Qc_E111(n_peak_TE111)*Rs/(pi*etha)
= 1.111186527621784e+04
Qc*Rs/(pi*etha) should be same as Qc*atand/lambda0
1e2*(Qc_peak- Qc_E111(n_peak_TE111)*Rs/(pi*etha))/Qc_peak
99.99..
100% error
Q caused by resonance in dielectric
Qd=1/tand
And the total Q
Q_f0=1/(1/Qc_peak+1/Qd)
Qd = 2000
Q_f0 = 1.742993528344060e+03