top of page
exercise 1.5
ZIP with MATLAB scripts and note:
example 1.1 notes:
eps0=8.854e-12 % [F/m] electric permittivity
mu0=4*pi*1e-7 % [H/m] magnetic permeability
c0=1/(eps0*mu0)^.5 % light velocity
epsr1=1 % relative permittivity air
epsr2=2 % relative permittivity dielectric slice
c1=c0/(epsr1)^.5 % wave velocity in air = c0
c2=c0/(epsr2)^.5 % wave velocity in dielectric < c0
D=1 % x size of overall scenario
tmax=3*4.5e-9 % to visualise pulses
% tmax=1.2*14.5e-9 % to measure reflection coeff
r=1
NX=1e3 % space resolution
dx=D/(NX-1)
dt=r*dx/c0 % time resolution
NT=ceil(tmax/dt) % amount time samples
x=[0:dx:D];
t0=60*dt
s=10*dt % pulse width
u1=zeros(NX,1);
u2=zeros(NX,1);
u3=zeros(NX,1);
a1=(c1*dt/dx)^2*ones(NX,1); % update eq parameter
% 2 interfaces : slice of dielectric
d1=D/3;
d2=2*D/3;
idx2=find(x>d1);
idx3=find(x<d2);
idx=intersect(idx2,idx3);
a1(idx)=(c2*dt/dx)^2;
% coefficients needed to approximate fields
a2=(sqrt(a1(NX))-1)/(sqrt(a1(NX))+1); % absorbing boundary coefficient, right hand side
a3=(sqrt(a1(1))-1)/(sqrt(a1(1))+1); % absorbing boundary coefficient, left hand side
hf1=figure(1)
hf1.OuterPosition=[672 550 776 513]
ax1=gca
axis([0 1 -2 2])
rectangle('Position',[d1 -2 d2-d1 4],...
'FaceColor',[.5 .8 .8],...
'EdgeColor','none')
grid on
hold on
ax1.DataAspectRatio=[1 4 1]
text(.1,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
text(.4,1.5,['\epsilon_r_2= ' num2str(epsr2)],'FontSize',14)
text(.8,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
title('1D FDTD : dielectric slice')
xlabel('x [m]')
ylabel('u','Rotation',0)
% generate initial value
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
t=dt;
u3(1)=exp(-(t-t0)^2/(2*s^2)); % generator left hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
L1=line(x,u3)
L2=zeros(1,numel(x)); % logging traces
L2=[L2;u3'];
L1.LineWidth=1.5
G1=text(1.2,1.5,[num2str(t/1e-9,2) ' ns'])
drawnow
hold off
for n=1:1:NT
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
L2=[L2;u3'];
t=n*dt;
u3(1)=exp(-(t-t0)^2/(2*s^2)); % source left hand side
% absorbing boundary right hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
% absorbing boundary left hand side
% turn it on after whole transmitted pulse has left
if n>120
u3(1)=u2(2) + a3*(u3(2)-u2(1));
end
L1.YData=u3
G1.String=[num2str(t/1e-9,2) ' ns'];
drawnow
% variables time shift
u1=u2;
u2=u3;
end
OK, now we have a way to simulate a slice of dielectric. But I chose an arbitrary thickness, and what the question is asking is about the effective reflection coefficient, from generator, on left hand side.
how to calculate reflection coefficient
in example 1.1 I calculated R and T out of simulation choosing t span so that
there were exactly and only 2 pulses present right at the end of the simulation, just
next on the dielectric interface. Therefore one could then use max just to find
wave peak values, and from there get reflection and transmission
coefficients. And the results matched the expected analytical calculation.
The problem is that to choose tmax one has to know the incident signal.
When measuring power, many times one doesn't really know the shape of the signal, and actually, when handling power values, there's no need to really know the exact shape
of incident and reflected signals.
So let's calculate the Reflection and Transmission coefficients just as ratio pf how much power went right and then measuring and integrating how much power has eventually bounced back.
This can be achieved by logging data on each loop and then processing logged data.
L2(1,:)=[]; % remove initial trace
L21=L2(:,[1:min(idx2)]); % traces points in zone 1 : air left side of dielectric
L22=L2(:,[min(idx2)+1:max(idx3)]); % traces points in zone 2 : in dielectric
L23=L2(:,[max(idx3)+1:end]); % traces points in zone 3 : air right side of dielectric
this is the generated pulse
t1=dt*[1:100];
v1=exp(-(t1-t0).^2/(2*s^2));
figure(2);plot(t1,v1)
grid on
title('transmitted pulse')
xlabel('t')
radiated EH signals do not have DC.
Instead of generating a pulse with all positive values, one could have generated and + - pulse, which are the actual pulses produced by any antenna.
But for this particular exercise, we can save the AC to DC step.
and this is the overall signal sent and received, on the left hand side:
s2=L21(:,100); % detector position
figure(3);plot(s2)
ax3=gca
grid on
title('detected pulses')
xlabel('t/dt : time steps')
s3=L21(:,10); % detector position
hold on
plot(s3)
legend({'detector position x=100','detector position x=10'})
the closer to the transmitter than we put the detector the more accurate would be distance measurements to interfaces.
Let's work with power
p1=s3.^2;
figure(4)
ax4=gca
plot(p1)
title('signal power')
grid on
when changing to dBW we even catch the last echo that in linear goes unnoticed, to bare eye.
p2=10*log10(p1);
figure(5)
ax5=gca
plot(p2)
title('signal power dBW')
grid on
[pks,locs]=findpeaks(p2,'MinPeakProminence',40,'SortStr','descend')
hold on
plot(locs,pks,'rd')
text(locs+100,pks,num2str((1:numel(pks))'))
pks =
0
-15.280112110706945
-15.582918729180538
-46.246646755489962
-76.940383695043707
locs =
70
719
1659
2599
3539
prefl_dB =
-15.280112110706945
-15.582918729180538
-46.246646755489962
-76.940383695043707
prefl =
0.029647548550376
0.027650827131676
0.000023732053814
0.000000020228405
prx =
0.057322127964270
ptx_dB =
0
ptx =
1
R =
0.239420400058704
reflection_coefficient^2 [dB] = transmitter power[dB] - bounced back power[dB]
reflection_coefficient = 10^(sum(reflected)/10)
prefl_dB=pks([2:end])
prefl=10.^(prefl_dB/10)
% received
prx=sum(prefl)
% transmitted
ptx_dB=pks(1)
ptx=10^(ptx_dB/10)
R=(prx/ptx)^.5
the analytical solution for a slice of dielectric uses single tone calculations:
1. the slice of dielectric is illuminated with single frequency pulse,
2. the slab is lambda/4 thick, lambda=lambda0/er^.5
3. the analytical reflection coefficient for a lambda/4 thick slice is
(etha^2-etha0^2)/(etha^2+etha0^2)
leaving the theoretical result to (1/er-1)/(1/er+1) that for er=2
means R=1/3
I have obtained R=.233 , but I chose a randomly thick dielectric
and I have used an equivalent all + gaussian pulse.
here we have 2 options,
1.- variating slice thickness, one can find the optimal dielectric
thickness, that would reflect highest power, equivalent to lambda/4
transformer.
2.- illuminate with single frequency pulse, large enough amount of cycles
setting dielectric thickness to lambda/4, and compare results with
theoretical.
Let's check both options
1.-
f0=1e9
lambda0=c0/f0
lambda=lambda0/epsr2^.5
u3(1)=exp(-(t-t0)^2/(2*s^2)); % generator left hand side
dD=lambda0/(4*epsr2^.5) % dielectric slice thickness
u1=zeros(NX,1);
u2=zeros(NX,1);
u3=zeros(NX,1);
a1=(c1*dt/dx)^2*ones(NX,1);
% slice of dielectric
d1=D/3; % location of left side of dielectric interface
d2=d1+dD;
idx2=find(x>d1);
idx3=find(x<d2);
idx=intersect(idx2,idx3);
a1(idx)=(c2*dt/dx)^2;
% coefficients needed to approximate fields
a2=(sqrt(a1(NX))-1)/(sqrt(a1(NX))+1); % absorbing boundary coefficient, right hand side
a3=(sqrt(a1(1))-1)/(sqrt(a1(1))+1); % absorbing boundary coefficient, left hand side
%%
hf1=figure(6)
hf1.OuterPosition=[672 550 776 513]
ax1=gca
axis([0 1 -2 2])
rectangle('Position',[d1 -2 dD 4],...
'FaceColor',[.5 .8 .8],...
'EdgeColor','none')
grid on
hold on
ax1.DataAspectRatio=[1 4 1]
text(.1,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
text(.4,1.5,['\epsilon_r_2= ' num2str(epsr2)],'FontSize',14)
text(.8,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
title('1D FDTD : dielectric slice')
xlabel('x [m]')
ylabel('u','Rotation',0)
% generate initial value
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
t=dt;
u3(1)=cos(t) % generator
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
L1=line(x,u3)
L2=zeros(1,numel(x)); % logging traces
L2=[L2;u3'];
L1.LineWidth=1.5
G1=text(1.2,1.5,[num2str(t/1e-9,2) ' ns'])
drawnow
hold off
for n=1:1:NT
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
L2=[L2;u3'];
t=n*dt;
u3(1)=exp(-(t-t0)^2/(2*s^2)); % source left hand side
% absorbing boundary right hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
% absorbing boundary left hand side
% turn it on after whole transmitted pulse has left
if n>120
u3(1)=u2(2) + a3*(u3(2)-u2(1));
end
L1.YData=u3
% L2=[L2;u3'];
G1.String=[num2str(t/1e-9,2) ' ns'];
drawnow
% variables time shift
u1=u2;
u2=u3;
end
L2(1,:)=[]; % remove initial trace
L21=L2(:,[1:min(idx2)]); % traces points in zone 1 : air left side of dielectric
L22=L2(:,[min(idx2)+1:max(idx3)]); % traces points in zone 2 : in dielectric
L23=L2(:,[max(idx3)+1:end]);
s3=L21(:,10);
p1=s3.^2;
p2=10*log10(p1);
figure
ax4=gca
plot(p2)
title('signal power')
grid on
[pks,locs]=findpeaks(p2,'MinPeakProminence',10,'SortStr','descend')
hold on
plot(locs,pks,'rd')
text(locs,pks+10,num2str((1:numel(pks))'))
% reflection_coefficient^2 [dB] = transmitter power[dB] - bounced back power[dB]
% reflection coefficient = 10^(sum(reflected)/10)
prefl_dB=pks([2:end])
prefl=10.^(prefl_dB/10)
% received
prx=sum(prefl)
% transmitted
ptx_dB=pks(1)
ptx=10^(ptx_dB/10)
R=(prx/ptx)^.5
pks =
1.0e+02 *
0
-0.152801121107069
-0.155456524101453
-0.461149817767802
-0.766702293357774
-1.071936256983315
-1.364720290822214
locs =
70
719
866
1014
1161
1308
1455
prefl_dB =
1.0e+02 *
-0.152801121107069
-0.155456524101453
-0.461149817767802
-0.766702293357774
-1.071936256983315
-1.364720290822214
prefl =
0.029647548550376
0.027889116662274
0.000024462555406
0.000000021526681
0.000000000019083
0.000000000000023
prx =
0.057561149313842
ptx_dB =
0
ptx =
1
R =
0.239919047417752
As I am showing a bit further down, again this reflection coefficient is not the expected one, but it's ok because the solution supplied in the solutions manual is for a single frequency incident field, yet gaussian pulses are broadband compared to single frequency bursts.
Before going for option 2, checking reflection R for different slice thickness
ND=50 % resolution sweep dD
dD=linspace(lambda/10,D/2,ND)
% left side of interface remains in same place
% d1=D/3; % location of left side of dielectric interface
R2=[] % log obtained reflection coefficients for each slice thickness
a1=(c1*dt/dx)^2*ones(NX,1);
for k2=1:1:ND
k2 % this loop is time consuming, this way one knows on what lap is the loop on
% dD=lambda0/(4*epsr2^.5) % dielectric slice thickness
u1=zeros(NX,1);
u2=zeros(NX,1);
u3=zeros(NX,1);
% slice of dielectric
d2=d1+dD(k2);
idx2=find(x>d1);
idx3=find(x<d2);
idx=intersect(idx2,idx3);
a1(idx)=(c2*dt/dx)^2;
% coefficients needed to approximate fields
a2=(sqrt(a1(NX))-1)/(sqrt(a1(NX))+1); % absorbing boundary coefficient, right hand side
a3=(sqrt(a1(1))-1)/(sqrt(a1(1))+1); % absorbing boundary coefficient, left hand side
% generate initial value
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
t=dt;
u3(1)=exp(-(t-t0)^2/(2*s^2)); % generator left hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
L1=line(x,u3)
L2=zeros(1,numel(x)); % logging traces
L2=[L2;u3'];
for n=1:1:NT
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
L2=[L2;u3'];
t=n*dt;
u3(1)=exp(-(t-t0)^2/(2*s^2)); % source left hand side
% absorbing boundary right hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
% absorbing boundary left hand side
% turn it on after whole transmitted pulse has left
if n>120
u3(1)=u2(2) + a3*(u3(2)-u2(1));
end
L1.YData=u3;
u1=u2;
u2=u3;
end
L2(1,:)=[]; % remove initial trace
L21=L2(:,[1:min(idx2)]); % traces points in zone 1 : air left side of dielectric
L22=L2(:,[min(idx2)+1:max(idx3)]); % traces points in zone 2 : in dielectric
L23=L2(:,[max(idx3)+1:end]);
s3=L21(:,10);
p1=s3.^2;
p2=10*log10(p1);
[pks,locs]=findpeaks(p2,'MinPeakProminence',10,'SortStr','descend');
prefl_dB=pks([2:end]);
prefl=10.^(prefl_dB/10);
% received
prx=sum(prefl);
% transmitted
ptx_dB=pks(1);
ptx=10^(ptx_dB/10);
R=(prx/ptx)^.5;
R2=[R2 R];
end
figure
plot(dD,R2)
grid on
title('reflection coeff over dielectric thickness')
xlabel('dD [m]')
2.- pulse is now burst single frequency
% now this is the pulse
t1=dt*[1:ceil(1/f0/dt)];
v3=zeros(1,numel(t1));
v3([1:ceil(1/f0/dt)])=sin(2*pi*f0*(t1+dt));
figure;plot(t1,v3)
grid on
title('transmitted pulse : not adjusted')
xlabel('t')
hold on
v4=exp(-(t1-t0).^2/(2*s^2));
plot(t1,v4)
legend({'single freq burst','gauss pulse'})
to avoid sharp cuts on the signal line, one sometimes has to adjust as follows:.
%% adjusting frequency to avoid sharp changes
N0=ceil(1/f0/dt)
f0=1/(dt*N0)
lambda0=c0/f0
lambda=lambda0/epsr2^.5
dD=lambda/4
a1=(c1*dt/dx)^2*ones(NX,1);
u1=zeros(NX,1);
u2=zeros(NX,1);
u3=zeros(NX,1);
% slice of dielectric
d2=d1+dD;
idx2=find(x>d1);
idx3=find(x<d2);
idx=intersect(idx2,idx3);
a1(idx)=(c2*dt/dx)^2;
% coefficients needed to approximate fields
a2=(sqrt(a1(NX))-1)/(sqrt(a1(NX))+1); % absorbing boundary coefficient, right hand side
a3=(sqrt(a1(1))-1)/(sqrt(a1(1))+1); % absorbing boundary coefficient, left hand side
t1=dt*[1:N0+1];
v3=zeros(1,numel(t1));
v3([1:numel(t1)])=sin(2*pi*f0*(t1))-sin(2*pi*f0*dt);
figure;plot(t1,v3)
grid on
title('transmitted pulse')
xlabel('t')
hold on
v4=exp(-(t1-t0).^2/(2*s^2));
plot(t1,v4)
legend({'single freq burst','gauss pulse'})
v5=zeros(1,NT);
v5([1:numel(v3)])=v3
figure;plot(v5)
N0 =
300
f0 =
9.983194735184102e+08
lambda0 =
0.300300300300300
lambda =
0.212344378734699
dD =
0.053086094683675
NT=NT/2 % on 2nd half of time span there are no echoes
u1=zeros(NX,1);
u2=zeros(NX,1);
u3=zeros(NX,1);
hf1=figure
hf1.OuterPosition=[672 550 776 513];
ax1=gca
axis([0 1 -2 2])
rectangle('Position',[d1 -2 dD 4],...
'FaceColor',[.5 .8 .8],...
'EdgeColor','none')
grid on
hold on
ax1.DataAspectRatio=[1 4 1];
text(.1,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
text(.4,1.5,['\epsilon_r_2= ' num2str(epsr2)],'FontSize',14)
text(.8,1.5,['\epsilon_r_1= ' num2str(epsr1)],'FontSize',14)
title('1D FDTD : dielectric slice')
xlabel('x [m]')
ylabel('u','Rotation',0)
% generate initial value
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
t=dt;
u3(1)=v5(1)
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
L1=line(x,u3);
L2=zeros(1,numel(x)); % logging traces
L2=[L2;u3'];
L1.LineWidth=1.5;
G1=text(1.2,1.5,[num2str(t/1e-9,2) ' ns'])
drawnow
hold off
for n=1:1:NT
u3([2:NX-1])=a1([2:NX-1]).*(u2([3:NX]) - 2*u2([2:NX-1]) + ...
u2([1:NX-2])) + 2*u2([2:NX-1]) - u1([2:NX-1]);
L2=[L2;u3'];
t=n*dt;
u3(1)=v5(n);
% absorbing boundary right hand side
u3(NX)=u2(NX-1) + a2*(u3(NX-1)-u2(NX));
% absorbing boundary left hand side
% turn it on after whole transmitted pulse has left
if n>N0+1
u3(1)=u2(2) + a3*(u3(2)-u2(1));
end
L1.YData=u3
G1.String=[num2str(t/1e-9,2) ' ns'];
drawnow
% variables time shift
u1=u2;
u2=u3;
end
prx =
0.113063428646975
ptx_dB =
0.180024980188956
ptx =
1.042323424717285
R =
0.329351637190975
pretty close to R=1/3 obtained when using the theoretical formula
(etha^2-etha0^2)/(etha^2+etha0^2)
suggested in the solutions manual
(etha^2-etha0^2)/(etha^2+etha0^2) = (1/er-1)/(1/er+1)
that for er2=2, er1=1 makes R=1/3
Therefore (etha^2-etha0^2)/(etha^2+etha0^2) is a single frequency
calculation, it does not apply to broadband signals.
Closing comments:
1.- the thickness of the ideal (lossless) dielectric does not really affect reflection coefficient.
2.- effect of large er : when er >> 5 and dD>.5 the pulse inside the dielectric , when
dielectric thick enough shows distortion : accumulated error, derivative
approximation is then not good enough.
3.- If for instance c0/10 hopping to simplify things, or improve accuracy : it doesn't
work. It doesn't remove wave distortion when slab thick enough,
and one has to apply same inverted factor to tmax, or the simulation ends
up way too early.
4.- The thickness of a metal would really affect reflection coefficient
but again beyond a few skin depths, all field that gets through 1st
interface cannot make it back.
5.- since theoretical gives R(1/er) the next logic thing to do would be to sweep er, for a fixed dielectric slice
6.- mlx live script gets too busy and does not produce any results for way too much time.
bottom of page