gusucode.com > PPML_v1.2 > 1d_tm/SM_1d_tm.m
%%% function [rl,rr,tlr,trl] = SM_1d_tm(a,L,... epssup,epssub,epsxA,epszA,epsxB,epszB,f,d,... halfnpw,k0,kpar) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is free software distributed under the BSD licence (see the % containing folder). % However, shall the results obtained through this code be included % in an academic publication, we kindly ask you to cite the source % website. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculates the S-matrix. % % 1-d structure, TM polarization. % % Simone Zanotto, Orsay Oct. - Dec 2012; Firenze Feb. 2016 % ------------------------------------------------------------------------- % INPUTS % % a -> periodicity % L -> number of layers excluded superstrate and substrate % epssup, epssub -> super-and sub-strate permittivities % epsxA -> in-plane component of dielectric tensor % for material A in the internal layers % [vector with L components] % epsxB -> idem, for material B % epszA, epszB -> idem, out-of-plane components % f -> duty cycle in the internal layers % [vector with L components] % d -> thicknesses of the layers, % including super- and sub-strate % [vector with (L+2) components] % halfnpw -> half number of harmonic waves used in calculation % k0 -> wavevector in vacuum (= 2*pi/lambda0) % kpar -> x-projection of the incident wavevector % % ------------------------------------------------------------------------- % OUTPUTS % % The coefficients connect the amplitudes of the Hy fields, calculated at the interface between % superstrate and first internal layer, or substrate and last internal layer. % % rl -> reflection coefficient from superstrate to superstrate ("left reflectance") % rr -> reflection coefficient from substrate to substrate ("right reflectance") % tlr -> transmission coefficient from superstrate to substrate ("left to right") % trl -> transmission coefficient from substrate to superstrate ("right to left") % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NOTES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % *has been tested extensively only when superstrate and substrate have equal % epsilon and thickness % % *the lengths (a,d) and inverse lengths (k0,kpar) must be set in the same % units (e.g. microns and inverse microns, respectively) % % *halfnpw = 0 sets the number of harmonic waves to 1, e.g. the scattering % matrix reduces to the ordinary 2x2 formalism for unpatterned multilayers % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n = -halfnpw:halfnpw; npw = size(n,2); kx = (2*pi/a)*n + kpar; % Initializing the work variables q = zeros(npw,npw,L+2); phi = zeros(npw,npw,L+2); A = zeros(npw,npw,L+2); epsx = zeros(npw,npw,L+2); etaz = zeros(npw,npw,L+2); % Eigensolutions for super- and substrate q(:,:,1) = diag(sqrt_whittaker(epssup*k0^2 - kx.*kx)); phi(:,:,1) = eye(npw); A(:,:,1) = diag(k0^2 - kx.*kx/epssup)/q(:,:,1); epsx(:,:,1) = eye(npw)/epssup; etaz(:,:,1) = eye(npw)*epssup; q(:,:,L+2) = diag(sqrt_whittaker(epssub*k0^2 - kx.*kx)); phi(:,:,L+2) = eye(npw); A(:,:,L+2) = diag(k0^2 - kx.*kx/epssub)/q(:,:,L+2); epsx(:,:,L+2) = eye(npw)/epssub; etaz(:,:,L+2) = eye(npw)*epssub; % Eigensolutions for internal layers if L > 0 for l = 1:L F = zeros(npw); for i = 1:npw % Fourier transform of square centered at the origin for j = 1:npw if (i == j) F(i,j) = f(l); else F(i,j) = sin(pi*f(l)*(n(i)-n(j)))/(pi*(n(i)-n(j))); end end end epsx(:,:,l+1) = (1/epsxB(l) - 1/epsxA(l))*F + 1/epsxA(l)*eye(npw); etaz(:,:,l+1) = ( epszB(l) - epszA(l))*F + epszA(l)*eye(npw); % Eigensolution calculation [phhi,qq] = eig(epsx(:,:,l+1)\(k0^2*eye(npw) - diag(kx)*(etaz(:,:,l+1)\diag(kx)))); q(:,:,l+1) = diag(sqrt_whittaker(diag(qq))); phi(:,:,l+1) = phhi; A(:,:,l+1) = (k0^2*eye(npw) - diag(kx)*(etaz(:,:,l+1)\diag(kx)))*phhi/q(:,:,l+1); end end % Forward propagation S1 = zeros(npw,npw,L+2); S2 = zeros(npw,npw,L+2); S1(:,:,1) = eye(npw); for l = 1:L+1 [S1(:,:,l+1),S2(:,:,l+1)] = smpropag_fw(S1(:,:,l),S2(:,:,l),... phi(:,:,l),phi(:,:,l+1),... A(:,:,l),A(:,:,l+1),... exp(1i*diag(q(:,:,l)*d(l))),exp(1i*diag(q(:,:,l+1)*d(l+1)))); end % Backward propagation S3 = zeros(npw,npw,L+2); S4 = zeros(npw,npw,L+2); S4(:,:,L+2) = eye(npw); for ll = 1:L+1 l = L+3-ll; [S3(:,:,l-1),S4(:,:,l-1)] = smpropag_bw(S3(:,:,l),S4(:,:,l),... phi(:,:,l),phi(:,:,l-1),... A(:,:,l),A(:,:,l-1),... exp(1i*diag(q(:,:,l)*d(l))),exp(1i*diag(q(:,:,l-1)*d(l-1)))); end rl = S3(halfnpw+1,halfnpw+1,1) *exp(-1i*q(halfnpw+1,halfnpw+1,1)*d(1)); % left -> left reflectance rr = S2(halfnpw+1,halfnpw+1,L+2) *exp(-1i*q(halfnpw+1,halfnpw+1,L+2)*d(L+2)); % right -> right reflectance tlr = S1(halfnpw+1,halfnpw+1,L+2)*exp(-1i*q(halfnpw+1,halfnpw+1,1)*d(1)); % left -> right transmittance trl = S4(halfnpw+1,halfnpw+1,1) *exp(-1i*q(halfnpw+1,halfnpw+1,L+2)*d(L+2)); % right -> left transmittance