gusucode.com > PPML_v1.2 > 2d/rxx_2d_Lshape.m
%%
function [rss,rsp,rps,rpp] = rxx_2d_Lshape(a,L,...
epssup,epssub,epsA,epsB,f1,f2,f3,d,...
halfnpw,k0,kparx,kpary)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 polarization-resolved reflection coefficients.
%
% 2-d L-shaped structure
%
% Simone Zanotto, Pisa, Jan-Feb 2013, Milano, June 2015, and Firenze, nov 2016
% -------------------------------------------------------------------------
% INPUTS
%
% a -> periodicity
% L -> number of layers excluded superstrate and substrate
% epssup, epssub -> super-and sub-strate permittivities
% epsA -> dielectric constant
% for material A in the internal layers
% [vector with L components]
% epsB -> idem, for material B
% f1, f2, f3 -> define the shape of L inclusion (see notes)
% [vector with L components]
% d -> thicknesses of the layers,
% including super- and sub-strate
% [vector with (L+2) components]
% halfnpw -> truncation order (square scheme adopted)
% k0 -> wavevector in vacuum (= 2*pi/lambda0)
% kparx -> x-projection of the incident wavevector
% kpary -> y-projection of the incident wavevector
%
%%% -----------------------------------------------------------------------
%%% OUTPUTS
%%%
%%%
%%% rxx -> polarization resolved reflection amplitudes.
%%% ( zeroth reflection order; phases at the interface
%%% between superstrate and first layer
%%% (or substrate if there are no internal layers) )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% *the lengths (a,d) and inverse lengths (k0,kpar) must be set in the same
% units (e.g. microns and inverse microns, respectively)
%
% *epssup must be real (otherwise the incident waves are ill-defined)
%
% *epssub can be either real or complex.
%
% *halfnpw = 0 sets the number of harmonic waves to 1, e.g. the scattering
% matrix reduces to the ordinary 2x2 formalism for unpatterned multilayers
%
% *thickness of superstrate and substrate do not influence the result.
% The electromagnetic problem is solved de facto assuming semi-infinite
% super- and substrate boundary conditions.
% You find them in view of further development (field profile calculation).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (min(f1-f3) < -1e-3) || (min(f2-f3) < -1e-3)
error('rxx_2d_Lshape:fconstr','Broken constraints on fs')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Building the square lattice
npw = (2*halfnpw+1)^2;
nx = zeros(npw,1); ny = zeros(npw,1);
i = 1;
for ix = -halfnpw:halfnpw
for iy = -halfnpw:halfnpw
nx(i) = ix; ny(i) = iy;
i = i+1;
end
end
normm = nx.^2+ny.^2;
[norm_sort,indices] = sort(normm);
nx_sort = nx(indices);
ny_sort = ny(indices);
nx = nx_sort;
ny = ny_sort;
clear nx_sort ny_sort indices norm_sort normm ix iy i
%%% now nx ny npw are defined
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kx = (2*pi/a)*nx + kparx;
ky = (2*pi/a)*ny + kpary;
% Initializing the work variables
q = zeros(2*npw,2*npw,L+2);
phi = zeros(2*npw,2*npw,L+2);
A = zeros(2*npw,2*npw,L+2);
% Eigensolutions for super- and substrate
q(1:npw,1:npw,1) = diag(sqrt_whittaker(epssup*k0^2 - kx.*kx - ky.*ky));
q(npw+1:2*npw,npw+1:2*npw,1) = q(1:npw,1:npw,1);
phi(:,:,1) = eye(2*npw);
eps = epssup;
kstorto = [(diag(ky)/eps)*diag(ky), -(diag(ky)/eps)*diag(kx);...
-(diag(kx)/eps)*diag(ky), (diag(kx)/eps)*diag(kx)];
A(:,:,1) = ((eye(2*npw)*k0^2 - kstorto)*phi(:,:,1))/q(:,:,1);
q(1:npw,1:npw,L+2) = diag(sqrt_whittaker(epssub*k0^2 - kx.*kx - ky.*ky));
q(npw+1:2*npw,npw+1:2*npw,L+2) = q(1:npw,1:npw,L+2);
phi(:,:,L+2) = eye(2*npw);
eps = epssub;
kstorto = [(diag(ky)/eps)*diag(ky), -(diag(ky)/eps)*diag(kx);...
-(diag(kx)/eps)*diag(ky), (diag(kx)/eps)*diag(kx)];
A(:,:,L+2) = ((eye(2*npw)*k0^2 - kstorto)*phi(:,:,L+2))/q(:,:,L+2);
% Eigensolutions for internal layers
if L > 0
for l = 1:L
F = zeros(npw);
% Quite crazy ordinary Fourier transforms of an L shaped inclusion
for i = 1:npw
for j = 1:npw
F(i,j) = (f1(l)-f3(l))*f3(l)*sinc((f1(l)-f3(l))*(nx(i)-nx(j)))*sinc(f3(l)*(ny(i)-ny(j))) ...
*exp(1i*pi*(f1(l)+f3(l)-1)*(nx(i)-nx(j)))*exp(1i*pi*(f3(l)-1)*(ny(i)-ny(j))) ...
+ (f2(l)-f3(l))*f3(l)*sinc(f3(l)*(nx(i)-nx(j)))*sinc((f2(l)-f3(l))*(ny(i)-ny(j))) ...
*exp(1i*pi*(f3(l)-1)*(nx(i)-nx(j)))*exp(1i*pi*(f2(l)+f3(l)-1)*(ny(i)-ny(j))) ...
+ f3(l)*f3(l)*sinc(f3(l)*(nx(i)-nx(j)))*sinc(f3(l)*(ny(i)-ny(j))) ...
*exp(1i*pi*(f3(l)-1)*(nx(i)-nx(j)))*exp(1i*pi*(f3(l)-1)*(ny(i)-ny(j)));
end
end
eps = (epsB(l) - epsA(l))*F + epsA(l)*eye(npw);
% Crazy crossed mixed Fourier transforms of an L shaped inclusion
F1 = zeros(2*halfnpw+1);
F2 = zeros(2*halfnpw+1);
F3 = zeros(2*halfnpw+1);
G1 = zeros(2*halfnpw+1);
G2 = zeros(2*halfnpw+1);
G13 = zeros(2*halfnpw+1);
G23 = zeros(2*halfnpw+1);
for i = 1:2*halfnpw+1
for j = 1:2*halfnpw+1
F1(i,j) = f1(l)*sinc(f1(l)*(i-j))*exp(1i*pi*(f1(l)-1)*(i-j));
F2(i,j) = f2(l)*sinc(f2(l)*(i-j))*exp(1i*pi*(f2(l)-1)*(i-j));
F3(i,j) = f3(l)*sinc(f3(l)*(i-j))*exp(1i*pi*(f3(l)-1)*(i-j));
G1(i,j) = (1-f1(l))*sinc((1-f1(l))*(i-j))*exp(1i*pi*f1(l)*(i-j));
G2(i,j) = (1-f2(l))*sinc((1-f2(l))*(i-j))*exp(1i*pi*f2(l)*(i-j));
G13(i,j) = (f1(l)-f3(l))*sinc((f1(l)-f3(l))*(i-j))*exp(1i*pi*(f1(l)+f3(l)-1)*(i-j));
G23(i,j) = (f2(l)-f3(l))*sinc((f2(l)-f3(l))*(i-j))*exp(1i*pi*(f2(l)+f3(l)-1)*(i-j));
end
end
c0 = inv(eye(2*halfnpw+1)/epsA(l));
c1 = inv(eye(2*halfnpw+1)/epsA(l) + (1/epsB(l) - 1/epsA(l))*F1);
c2 = inv(eye(2*halfnpw+1)/epsA(l) + (1/epsB(l) - 1/epsA(l))*F2);
c3 = inv(eye(2*halfnpw+1)/epsA(l) + (1/epsB(l) - 1/epsA(l))*F3);
eps_xy = zeros(npw);
eps_yx = zeros(npw);
for i = 1:npw
for j = 1:npw
eps_xy(i,j) = c0(nx(i)+halfnpw+1,nx(j)+halfnpw+1) *G2(ny(i)+halfnpw+1,ny(j)+halfnpw+1) ...
+ c3(nx(i)+halfnpw+1,nx(j)+halfnpw+1)*G23(ny(i)+halfnpw+1,ny(j)+halfnpw+1) ...
+ c1(nx(i)+halfnpw+1,nx(j)+halfnpw+1) *F3(ny(i)+halfnpw+1,ny(j)+halfnpw+1);
eps_yx(i,j) = c0(ny(i)+halfnpw+1,ny(j)+halfnpw+1) *G1(nx(i)+halfnpw+1,nx(j)+halfnpw+1) ...
+ c3(ny(i)+halfnpw+1,ny(j)+halfnpw+1)*G13(nx(i)+halfnpw+1,nx(j)+halfnpw+1) ...
+ c2(ny(i)+halfnpw+1,ny(j)+halfnpw+1) *F3(nx(i)+halfnpw+1,nx(j)+halfnpw+1);
end
end
% Eigensolution calculation
kstorto = [(diag(ky)/eps)*diag(ky), -(diag(ky)/eps)*diag(kx);...
-(diag(kx)/eps)*diag(ky), (diag(kx)/eps)*diag(kx)];
kdritto = [diag(kx)*diag(kx), diag(kx)*diag(ky);...
diag(ky)*diag(kx), diag(ky)*diag(ky)];
epsilone = [eps_yx, zeros(npw); zeros(npw), eps_xy];
[phhi,qq] = eig(epsilone*(eye(2*npw)*k0^2 - kstorto) ...
- kdritto);
q(:,:,l+1) = diag(sqrt_whittaker(diag(qq)));
phi(:,:,l+1) = phhi;
A(:,:,l+1) = (eye(2*npw)*k0^2 - kstorto)*phhi/q(:,:,l+1);
end
end
% Forward propagation
S1 = zeros(2*npw,2*npw,L+2); S2 = zeros(2*npw,2*npw,L+2);
S1(:,:,1) = eye(2*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(2*npw,2*npw,L+2); S4 = zeros(2*npw,2*npw,L+2);
S4(:,:,L+2) = eye(2*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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some useful variables
kparvec = [kparx, kpary, 0];
krefl = [kparx, kpary, -sqrt(epssup*k0^2 - kparx^2 - kpary^2)];
svers = cross(kparvec,[0,0,1]);
svers = svers/norm(svers);
pvers = cross(krefl,svers);
pvers = pvers/norm(pvers);
% s-pol incidence on the 0 order from superstrate %
as = zeros(2*npw,L+2); bs = zeros(2*npw,L+2);
M = [kparx*kpary, epssup*k0^2 - kparx^2; ...
-epssup*k0^2 + kpary^2, -kparx*kpary]/(epssup*q(1,1,1));
abar = linsolve(M,[-kpary; kparx]); %s-pol
Exyz = vertcat(M*abar,([kpary, -kparx]*abar)/epssup);
normm = (Exyz')*Exyz;
abar = abar/sqrt(normm);
as(1,1) = abar(1)*exp(-1i*q(1,1,1)*d(1)); % inizialization of a vectors
as(npw+1,1) = abar(2)*exp(-1i*q(1,1,1)*d(1)); %
% propagation with waves from superstrate
for l = 1:L+2
as(:,l) = (eye(2*npw) - S2(:,:,l)*S3(:,:,l))\(S1(:,:,l)*as(:,1));
bs(:,l) = (eye(2*npw) - S3(:,:,l)*S2(:,:,l))\(S3(:,:,l)*S1(:,:,l)*as(:,1));
end
Ereflx = -(kparx*kpary*bs(1,1) + (epssup*k0^2 - kparx^2)*bs(npw+1,1))/(epssup*q(1,1,1));
Erefly = ((epssup*k0^2 - kpary^2)*bs(1,1) + kparx*kpary*bs(npw+1,1))/(epssup*q(1,1,1));
Ereflz = (kpary*bs(1,1) - kparx*bs(npw+1,1))/epssup;
Erefl = [Ereflx, Erefly, Ereflz];
rss = sum(Erefl.*svers); % s -> s refl coeff
rsp = sum(Erefl.*pvers); % s -> p refl coeff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% p-pol incidence on the 0 order from superstrate %
as = zeros(2*npw,L+2); bs = zeros(2*npw,L+2);
M = [kparx*kpary, epssup*k0^2 - kparx^2; ...
-epssup*k0^2 + kpary^2, -kparx*kpary]/(epssup*q(1,1,1));
abar = linsolve(M,[kparx; kpary]); % p-pol
Exyz = vertcat(M*abar,([kpary, -kparx]*abar)/epssup);
normm = (Exyz')*Exyz;
abar = abar/sqrt(normm);
as(1,1) = abar(1)*exp(-1i*q(1,1,1)*d(1)); % inizialization of a vectors
as(npw+1,1) = abar(2)*exp(-1i*q(1,1,1)*d(1)); %
% propagation with waves from superstrate
for l = 1:L+2
as(:,l) = (eye(2*npw) - S2(:,:,l)*S3(:,:,l))\(S1(:,:,l)*as(:,1));
bs(:,l) = (eye(2*npw) - S3(:,:,l)*S2(:,:,l))\(S3(:,:,l)*S1(:,:,l)*as(:,1));
end
Ereflx = -(kparx*kpary*bs(1,1) + (epssup*k0^2 - kparx^2)*bs(npw+1,1))/(epssup*q(1,1,1));
Erefly = ((epssup*k0^2 - kpary^2)*bs(1,1) + kparx*kpary*bs(npw+1,1))/(epssup*q(1,1,1));
Ereflz = (kpary*bs(1,1) - kparx*bs(npw+1,1))/epssup;
Erefl = [Ereflx, Erefly, Ereflz];
rps = sum(Erefl.*svers); % p -> s refl coeff
rpp = sum(Erefl.*pvers); % p -> p refl coeff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%