%==========================================================================
%
% Solve for Plane-Wave Propagation Amplitudes in a Layered Medium.
%
% [A B] = solve_layered_medium(n,D,L0,THi,TE) computes the complex
% wave amplitudes (A,B) for plane-wave propagation in a layered medium
% that is excited from the (-z) direction.
%
% n = Vector of complex indices of refraction for each region.
% D = Vector of distances for each planar boundary (m). Note that each
% element in D needs to be increasing.
% L0 = Free-space wavelength of excitation (m).
% THi = Angle of incidence relative to z-axis (rads).
% TE = Boolean flag to indicate either transverse electric (TE == 1) or
% transverse magnetic (TE == 0) polarization.
%
% Note that for TE waves, the incident E-field is Y-polarized, and
% all computed wave amplitudes are the E-FIELD amplitudes. For TM waves,
% (i.e.; TE == 0), the incident H-field is Y-polarized and all the wave
% amplitudes are H-FIELD amplitudes.
%
% [A B kx kz] = solve_layered_medium(...) also returns the X- and Z-
% components to the wave-vectors in each region. Note that due to
% phase-matching, the kx component is identical for ALL REGIONS. That
% means kx is a single, scalar value, while kz is a vector of elements
% representing each region.
%
% [...] = solve_layered_medium(... , lossFactor) specifies a loss factor
% for lossy media. By default, this value is set to 1e-9. Any region
% that attenuates the incident wave by this much is the truncation point
% for further reflections and transmissions beyond. This measure greatly
% helps to stabilize the algorithm in the presence of lossy media. We
% can certaintly make this value smaller, but it increases the risk of
% numerical instability.
%
% SYSTEM GEOMETRY:
%
% z = D(1) z = D(2) ... z = D(M-1)
% | | |
% | | |
% | | |
% A(1) | A(2) | | A(M)
% ------> | ------> | ... | ------> ^ +x
% | | | |
% | | | |
% <------ | <------ | ... | <------ | +z
% B(1) | B(2) | | B(M) ------>
% | | |
% | | |
% | | |
% | | |
% n(1) | n(2) | ... | n(M)
% | | |
%
%
% FUTEHER NOTES:
%
% The computation assumes an incident plane wave propagating in the
% +z-direction. The angle of incidence (THi) is with respect to the
% z-axis. Given the system geometry, it should be apparent that for an
% M-region system, there should be M-1 planar boundaries.
%
% Another important note to keep in mind is the phase relationship
% between the amplitudes in A and B. Each field component is
% phase-referenced against the planar boundary that INJECTS the field.
% So for example, A(M) is phase-referenced against the z = D(M-1) planar
% boundary. This helps to make the computations more numerically stable,
% as well as more meaningful when trying to calculate such parameters
% like the transmittance into Region M. The incident wave in Region 1 is
% referenced against z = 0.
%
% BEWARE: In theory, the values in D can be offest by any arbitrary
% value. However, I have not validated this code under anything other
% than D(1) = 0. The algorithm can certainly be tweaked to make this
% work for arbitrary values, but has not yet been validated. So just be
% safe and assign D(1) = 0 whenever using this code.
%
% Note that the imaginary component to the index can be either
% POSITIVE OR NEGATIVE. By the convention within this code, a POSITIVE
% value for the imaginary index is LOSSY, while a NEGATIVE value
% is AMPLIFYING. BEWARE! THIS CONVENTION IS NOT UNIVERSAL! It is
% therefore important to be sure that the imaginary index fits the
% convention of this program. Another way to think of this is that a
% forward-propagating TE wave has the phasor/time-domain form:
%
% E+(x,z) = A * exp(+j(kx*X + kz*Z))
% E+(x,z,t) = A * cos(kx*X + kz*Z - 2*pi*f*t)
%
% Also note that overly lossy/amplifying media tend to behave in very
% quirky ways. The reason is because the computations become numerically
% unstable and error quickly accumulates. This makes the results
% meaningless if the user attempts to push the materials too far. For
% purely lossy regions, this issue is stabilized by simply truncating the
% waves after a certain point. For regions with amplification, this gets
% very tricky and the code will likely produce unreliable numbers if
% pushed too hard.
%
% This code has been validated against FDTD simulations for both TE and
% TM cases at arbitrary incidence angle in a 4-layered system. If any
% bugs are found, please email the author.
%
% REFERENCE:
% J. A. Kong, "Electromagnetic Wave Theory," pp 385-397 (2000).
%
% Note that I have made my own modifications to the notation and that I
% have also neglected to include arbitrary magnetic media in the
% solution. I have also tweaked the derivation to make the output more
% numerically stable from what Kong derives. See chapter 2 in my PhD
% dissertation for specific details.
%
%==========================================================================
%
% James Nagel
% Department of Electrical and Computer Engineering
% University of Utah, Salt Lake City, Utah
% nageljr@ieee.org
% Copyright August 23, 2011
function [A,B,kx,kz] = solve_layered_medium(n,D,L0,THi,TE,lossFactor)
% Extract parameters from inputs.
M = length(n); % Total number of regions.
k0 = 2*pi/L0; % Free-space wavenumber (rad/m).
k = k0*n; % Complex wavenumber in each region (rad/m).
% Instantaite forward and reverse amplitudes. For TE polarization, the
% units are V/m. For TM, the units are A/m.
A = zeros(M,1); % Forward amplitudes.
B = zeros(M,1); % Reverse amplitudes.
% Set the loss factor if not specified. This helps to stabilize the
% algorithm when dealing with lossy media. Any region where the E-field
% attenuates by this value is basically the cutoff point for all wave
% propagation beyond.
if nargin == 5
lossFactor = 1e-9;
end
% Incident plane wave has unit amplitude.
A(1) = 1;
% x-components of wavenumbers. These are all identical due to the
% phase-matching condition between boundaries.
kx = k0*n(1)*sin(THi);
% Square-root argument for the dispersion relation.
u = k.^2 - kx^2;
% z-components to the wavenumbers must all satisfy the dispersion relation
% within their respective media.
kz = sqrt(u);
% Special case for solving the dispersion relation. This represents an
% incident wave from a lossy medium. The square-root operator experiences
% a branch cut under this condition that needs to be handled.
% if ( real(u) < 0 ) & ( imag(u) < 0 )
% % VERIFY THIS EXPRESSION! - 7/12/11
% kz = -kz;
% end
% Check for length error on inputs.
if length(D) + 1 ~= M
error('D not correct length. Should be length(n) - 1.');
end
% Check to make sure D is properly increasing.
if any(D < 0) || any(diff(D) < 0)
error(['D must be positive and increasing.']);
end
% Solve for the reverse-propagating wave in Region 1. The recursion
% algorithm for this is unconditionally stable for lossy systems. See the
% function description below for more details.
B(1) = solve_layered_reflection(n,kz,D,TE,1);
%==========================================================================
% For the case where M >= 3, we have to be careful about directly applying
% the propagation-matrix method. Instead, begin with the numerically-
% stable amplitudes in Region 1 and iterate through each region directly.
% Do this by enforcing continuity at each boundary and by utilizing the
% reflection/transmission coefficients at each interface.
% When an extremely lossy region is encountered, just truncate the
% series and assume all the wave amplitudes go to zero beyond it. It will
% not be perfectly accurate for a mixed gain/loss system, but such systems
% are rarely encountered (if ever), and a highly-lossy system will become
% much more numerically stable.
%==========================================================================
for m = 1:M-1
% Generate three z-distance references. d0 is the z-value of the
% previous planar boundary. d1 is the current planar boundary, and d2
% is the next planar boundary.
% For Region 1, the distance is referenced at z = 0. Otherwise, use
% the previous distance value.
if m == 1
d0 = 0;
else
d0 = D(m-1);
end
% This is true for all regions.
d1 = D(m);
% For M-1, the "next" planar boundary does not exist. Since B(M) is
% always zero, it makes no difference what value we use for d2.
% Otherwise, just use the next D-value.
if m == M-1
d2 = 0;
else
d2 = D(m+1);
end
% Solve for exponential coefficients. For lossless systems, these
% will just be phase factors. For gain/lossy systems, these can
% easily lead to VERY LARGE or VERY SMALL numbers.
am = exp( +1j*kz(m) *(d1 - d0) );
bm = exp( -1j*kz(m+1)*(d1 - d2) );
% This is the key step that stabilizes the algorithm. Extremely lossy
% regions should effectively kill the signal, but numerical error
% tends to creep in and introduce energy where there should not be
% any. So when am is a very small number, that is our clue to end the
% series and assume that no more signals get through. So simply set
% everything beyond this point to zero and terminate the loop.
if abs(am) < lossFactor
B(m:end) = 0;
A(m+1:end) = 0;
break;
end
% Solve for normalized indices with respect to each planar boundary.
% These are needed in order to properly calculate the reflection
% coefficients.
m12 = n(m+1)/n(m); %
m21 = 1/m12; %
%
% Solve for Fresnel reflection coefficients. Note the convention "12"
% implies "Region 1 looking into Region 2" and vica versa. See the
% nested function provided within this code for details.
r12 = Fresnel_REFL(kz(m),kz(m+1),m12,TE);
r21 = Fresnel_REFL(kz(m+1),kz(m),m21,TE);
% Transmission coefficients.
t12 = 1 + r12;
t21 = 1 + r21;
% For the final iteration, the Bmbm term is zero. Otherwise, follow
% the formula.
if m == M-1
Bmbm = 0;
else
Bmbm = (B(m) - r12*A(m)*am)/t21;
end
% Solve for the wave amplitudes in the next region.
A(m+1) = r21*Bmbm + t12*A(m)*am;
B(m+1) = Bmbm/bm;
end
% End main function.
return;
%==========================================================================
%
% Complex Fresnel reflection coefficient at a planar boundary.
%
% r = Fresnel_REFL(k1z,k2z,m,TE) returns the complex-valued
% reflection coefficient for an EM wave incident on a planar boundary
% between two dielectric media. For TE polarization, the coefficient is
% with respect to the electric fields. For TM, it computes the magnetic
% field reflection coefficient.
%
% INPUTS:
%
% k1z = Z-component to the wavenumber in Region 1.
% k2z = Z-component to the wavenumber in Region 2.
% m = Complex-valued scalar representing the reduced index of
% refraction (n2/n1), where n1 is the index of refraction for
% Region 1 and n2 is the index for Region 2.
% TE = Boolean flag that tells Matlab if the incident polarization is
% transverse electric (TE == 1) or transverse magnetic (TE == 0).
% For TE polarization, the incident E-field is y-polarized. For
% TM, the H-field is y-polarized.
%
% System geometry is indicated below:
%
%
%
% (n1) (n2)
% REGION 1 | REGION 2
% |
% |
% |
% Reflected Wave \ |
% \ |
% \ | ^ x
% \ | |
% \ | |
% \| |
% ---------------| -------> z
% | /|
% THETA_i | / |
% \/ |
% / |
% / |
% / |
% Incident Wave |
% |
%
%
% REFERENCES:
% "Fundamentals of Applied Electromagnetics," 5th Edition, by Fawwaz
% Ulaby. In particular, see Equations (8.60) and (8.68).
%
%==========================================================================
function r = Fresnel_REFL(k1z,k2z,m,TE)
if TE
r = (k1z - k2z) / (k1z + k2z) ;
else
r = (k1z*m^2 - k2z) / (k1z*m^2 + k2z) ;
end
return
%==========================================================================
%
% Complex reflection coefficient at a multi-layered planar boundary.
%
% This function applies the recursive relation outlined in Kong,
% pages 388 - 390. Note that I actually tweaked the expressions to make
% the formula simpler and more numerically robust.
%
% INPUTS:
% n = Complex indices of refraction for each region in the stack.
% kz = z-components of the wavenumbers in each region (rad/m).
% D = Vector of distances for each planar boundary (m).
% TE = Transverse electric flag. Use TE == 1 for transverse electric.
% Otherwise, set TE == 0 for transverse magnetic.
% m = Iteration number in the recusion. Start at 1. When m == M,
% the recursion terminates.
%
%==========================================================================
function R_0 = solve_layered_reflection(n,kz,D,TE,m)
% Solve for fractional coefficients in the recursion relation. These vary
% depending on polarization.
if TE
u1 = kz(m) - kz(m+1) ;
u2 = kz(m) + kz(m+1) ;
else
u1 = kz(m)*n(m+1)^2 - kz(m+1)*n(m)^2 ;
u2 = kz(m)*n(m+1)^2 + kz(m+1)*n(m)^2 ;
end
% Solve for the first exponential term. For m == 1, D(m-1) = 0, and so ai
% is just unity.
if m == 1
ai = 1;
else
ai = exp(+1j*kz(m)*( D(m) - D(m-1) ) );
end
% Recursion terminator.
if m == length(D)
R_0 = ai*u1/u2;
return;
else
% If the recursion has not yet terminated, then we need this
% coefficient as well.
bi = exp(-1j*kz(m+1)*( D(m) - D(m+1) ) );
% Solve for the next recursion term.
r = solve_layered_reflection(n,kz,D,TE,m+1);
% Re-expression of equation 3.4.68 in Kong. Writing it this way is
% simpler and much more numerically stable.
R_0 = ai * ( u1 + u2*bi*r ) / (u2 + u1*bi*r) ;
end
return;