Загрузить файлы ''

fmmtd.m

@@ -0,0 +1,122 @@
+%{
+Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
+
+This file is part of GratingFMM.
+
+GratingFMM is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+GratingFMM is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
+%}
+%% description:
+% calculation of a grating S-matrix by the Fourier Modal Method
+%  in the case of the diffraction by 2D gratings being periodic in x and y
+%  dimensions of the Cartesian coordinates
+%% input:
+% xno, yno: numbers of Fourier harmonics
+% kx0, ky0: incident plane wave x and y projections (Bloch wavevector projections)
+% kgx, kgy: wavelength-to-period ratios (grating vector projections)
+% kh: grating depth multiplied by the vacuum wavenumber
+% eps1: permittivity of the substrate
+% eps2: permittivity of the superstrate
+% FE: Fourier matrix of the grating profile
+%% output:
+% SM: scattering matrix of size (2*no,2*no,2,2) where no = xno*yno
+%  is the total number of Fourier harmonics
+% block SM(:,:,1,1) corresponds to refelection from substrate to substrate
+% block SM(:,:,2,2) corresponds to refelection from superstrate to superstrate
+% block SM(:,:,2,1) corresponds to transmission from substrate to superstrate
+% block SM(:,:,1,2) corresponds to transmission from superstrate to substrate
+% first (no) components in each of the two first dimensions of the S-matrix
+%  correspond to the TE polarization, and indeces from (no+1) to (2*no) 
+%  correspond to the TM polarization
+% central harmonic index is ind_0 = (ceil(xno/2)-1)*yno+ceil(yno/2)
+% for example, an ampitude of the TE-polarized transmitted wave to i-th diffraction order
+%  from the substrate to the superstrate under the TM plane wave illumination
+%  with unit amplitude is SM(ind_0+i, no+ind_0, 2, 1)
+%% implementation
+function [SM] = fmmtd(xno, yno, kx0, ky0, kgx, kgy, kh, eps1, eps2, FE)
+	no = xno*yno;
+	ib1 = 1:no;
+	ib2 = no+1:2*no;
+	ib3 = 2*no+1:3*no;
+	ib4 = 3*no+1:4*no;
+
+		% wavevector projections
+	[kz1, kz2, kx, ky, kxy] = fmmtd_kxyz(xno, yno, kx0, ky0, kgx, kgy, eps1, eps2);
+	
+	ME = toeplitz2(FE{1,1},xno,yno);
+	MU = toeplitz2(FE{1,2},xno,yno);
+
+		% solve the eigenvalue problem:
+	EV = zeros(2*no,2*no);
+	HV = zeros(2*no,2*no);
+		% matrix for the electric field
+	EV(ib1,ib1) = ((kx').*MU).*ky;
+	EV(ib1,ib2) = -((kx').*MU).*kx;
+	EV(2*no*no+1:2*no+1:end) = EV(2*no*no+1:2*no+1:end) + 1;
+	EV(ib2,ib1) = ((ky').*MU).*ky;
+	EV(no+1:2*no+1:2*no*no) = EV(no+1:2*no+1:2*no*no) - 1;
+	EV(ib2,ib2) = -((ky').*MU).*kx;
+		% matrix for the magnetic field
+	HV(2*no*no+no+1:2*no+1:end) = HV(2*no*no+no+1:2*no+1:end) + kx.*ky;
+	HV(1:2*no+1:2*no*no) = -HV(2*no*no+no+1:2*no+1:end);
+	HV(ib1,ib2) = -ME;
+	HV(2*no*no+1:2*no+1:end) = HV(2*no*no+1:2*no+1:end) + kx.^2;
+	HV(ib2,ib1) = ME;
+	HV(no+1:2*no+1:2*no*no) = HV(no+1:2*no+1:2*no*no) - ky.^2;
+
+		% solve the eigenvalue problem for the electric field
+	[EV,MB] = eig(EV*HV);
+	beta = transpose(sqrt(diag(MB))); % row of eigenvalues
+	ind = angle(beta) < -1e-7; % check the branch of the square root
+	beta(ind) = -beta(ind);
+
+		% calculate the magnetic field amplitude eigen vectors
+	HV = (HV*EV).*(1./beta);
+
+		% calculate T-matrices
+	TS = fmmtd_calc_T(no,EV,HV,kx,ky,kxy,kz1,eps1); % substrate-grating T-matrix
+	TC = fmmtd_calc_T(no,EV,HV,kx,ky,kxy,kz2,eps2); % grating-cover T-matrix
+
+		% combine T-matrices
+	bexp = exp((1i*kh)*beta);
+	M1 = zeros(4*no,4*no);
+	M2 = zeros(4*no,4*no);
+
+	M1([ib1,ib2],[ib1,ib2]) = TS([ib3,ib4],[ib1,ib2]);
+	M1([ib3,ib4],[ib3,ib4]) = TC([ib1,ib2],[ib3,ib4]);
+	M1(ib1,ib3) = TS(ib3,ib3).*bexp(ib1);
+	M1(ib2,ib3) = TS(ib4,ib3).*bexp(ib1);
+	M1(ib1,ib4) = TS(ib3,ib4).*bexp(ib2);
+	M1(ib2,ib4) = TS(ib4,ib4).*bexp(ib2);
+	M1(ib3,ib1) = TC(ib1,ib1).*bexp(ib1);
+	M1(ib4,ib1) = TC(ib2,ib1).*bexp(ib1);
+	M1(ib3,ib2) = TC(ib1,ib2).*bexp(ib2);
+	M1(ib4,ib2) = TC(ib2,ib2).*bexp(ib2);
+
+	M2([ib1,ib2],[ib1,ib2]) = TS([ib1,ib2],[ib1,ib2]);
+	M2([ib3,ib4],[ib3,ib4]) = TC([ib3,ib4],[ib3,ib4]);
+	M2(ib1,ib3) = TS(ib1,ib3).*bexp(ib1);
+	M2(ib2,ib3) = TS(ib2,ib3).*bexp(ib1);
+	M2(ib1,ib4) = TS(ib1,ib4).*bexp(ib2);
+	M2(ib2,ib4) = TS(ib2,ib4).*bexp(ib2);
+	M2(ib3,ib1) = TC(ib3,ib1).*bexp(ib1);
+	M2(ib4,ib1) = TC(ib4,ib1).*bexp(ib1);
+	M2(ib3,ib2) = TC(ib3,ib2).*bexp(ib2);
+	M2(ib4,ib2) = TC(ib4,ib2).*bexp(ib2);
+
+		% S-matrix
+	SM = M2S(M1/M2);
+end
+%
+% END
+%

fmmtd_balance.m

@@ -0,0 +1,62 @@
+%{
+Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
+
+This file is part of GratingFMM.
+
+GratingFMM is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+GratingFMM is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
+%}
+%% description:
+% calculate the power balance (check the energy conservation law) in case
+% of the diffraction by 2D gratings being periodic in x and y dimensions
+%% input:
+% xno, yno: numbers of Fourier harmonics
+%  total number of Fourier harmonics is no = xno*yno;
+% V_inc: incident field amplitude matrix of size (2*no, 2)
+% V_dif: diffracted field amplitude matrix of size (2*no, 2)
+%  first index of V_inc, V_dif indicates diffraction harmonics
+%   with indices 1:no being TE orders and no+1:2*no being TM orders
+%   (0-th order index is ind_0 = (ceil(xno/2)-1)*yno+ceil(yno/2))
+%  second index of V_inc, V_dif, V_eff indicates whether the diffraction orders
+%   are in the substrate (V(:,1)) or in the superstrate (V(:,2))
+% kx0, ky0: incident plane wave wavevector x and y projections (Bloch wavevector projections)
+% kgx, kgy: wavelength-to-period ratios (grating vectors)
+% eps1, eps2: substrate and superstrate permittivities
+%% output:
+% if the incident field has propagating harmonics the function returns the
+%  normalized difference between the incident and diffractied field total
+%  power, otherwise (if the incident field is purely evanescent) it returns
+%  the total power carried by propagating diffraction orders
+%% implementation
+function [b] = fmmtd_balance(xno, yno, V_inc, V_dif, kx0, ky0, kgx, kgy, eps1, eps2)
+	no = xno*yno;
+	ib1 = 1:no; ib2 = no+1:2*no;
+
+	[kz1, kz2] = fmmtd_kxyz(xno, yno, kx0, ky0, kgx, kgy, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+
+	P_inc = sum( abs((V_inc(ib1,1)).^2).*real(kz1) + abs((V_inc(ib1,2)).^2).*real(kz2) ) ...
+				+ sum( abs((V_inc(ib2,1)).^2).*real(kz1/eps1) + abs((V_inc(ib2,2)).^2).*real(kz2/eps2) );
+	P_dif = sum( abs((V_dif(ib1,1)).^2).*real(kz1) + abs((V_dif(ib1,2)).^2).*real(kz2) ) ...
+				+ sum( abs((V_dif(ib2,1)).^2).*real(kz1/eps1) + abs((V_dif(ib2,2)).^2).*real(kz2/eps2) );
+
+	if (abs(P_inc) > 1e-15)
+		b = abs(P_dif/P_inc-1);
+	else
+		b = 0.5*P_dif;
+	end
+end
+%
+% END
+%

fmmtd_calc_T.m

@@ -0,0 +1,87 @@
+%{
+Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
+
+This file is part of GratingFMM.
+
+GratingFMM is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+GratingFMM is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
+%}
+%% description:
+% T-matrix of an interface between a photonic crystal and a homogeneous medium
+% the case of the non-collinear diffraction on 1D gratings
+%% input:
+% no: total number of Fourier harmonics
+% EV: matrix of the electric field eigenvectors for the photonic crystal
+% HV: matrix of the magnetic field eigenvectors for the photonic crystal
+% kx: row of grating vector projections in x direction
+% ky: row of grating vector projections in y direction
+% kxy: row of grating vector projections in xy plane
+% kz: row of plane wave propagation constants in the homogeneous medium
+% eps: permittivity the homogeneous medium
+%% output:
+% T: T-matrix of size 4*no by 4*no
+%% implementation:
+function [T] = fmmtd_calc_T(no, EV, HV, kx, ky, kxy, kz, eps)
+	T = zeros(4*no,4*no);
+		% pre-calculate combinations of wavevector projections
+	ikxy = 0.5./kxy;
+	ky_ikxy = ky.*ikxy;
+	kx_ikxy = kx.*ikxy;
+	ikxyz = ikxy./kz;
+	ky_ikxyz = transpose(ky.*ikxyz);
+	kx_ikxyz = transpose(kx.*ikxyz);
+	eky_ikxyz = eps*ky_ikxyz;
+	ekx_ikxyz = eps*kx_ikxyz;
+	kx_ikxy = transpose(kx_ikxy);
+	ky_ikxy = transpose(ky_ikxy);
+		% block indices
+	ib1 = 1:no;
+	ib2 = no+1:2*no;
+	ib3 = 2*no+1:3*no;
+	ib4 = 3*no+1:4*no;
+
+		% fill the T-matrix blocks
+	T(ib1,ib1) = EV(ib1,ib1).*ky_ikxy - EV(ib2,ib1).*kx_ikxy;
+	T(ib1,ib3) = T(ib1,ib1);
+	TM = HV(ib1,ib1).*kx_ikxyz + HV(ib2,ib1).*ky_ikxyz;
+	T(ib1,ib1) = T(ib1,ib1) + TM;
+	T(ib1,ib3) = T(ib1,ib3) - TM;
+	T(ib3,ib3) = T(ib1,ib1);
+	T(ib3,ib1) = T(ib1,ib3);
+
+	T(ib1,ib2) = EV(ib1,ib2).*ky_ikxy - EV(ib2,ib2).*kx_ikxy;
+	T(ib1,ib4) = T(ib1,ib2);
+	TM = HV(ib1,ib2).*kx_ikxyz + HV(ib2,ib2).*ky_ikxyz;
+	T(ib1,ib2) = T(ib1,ib2) + TM;
+	T(ib1,ib4) = T(ib1,ib4) - TM;
+	T(ib3,ib2) = T(ib1,ib4);
+	T(ib3,ib4) = T(ib1,ib2);
+
+	T(ib2,ib1) = -EV(ib1,ib1).*ekx_ikxyz - EV(ib2,ib1).*eky_ikxyz;
+	T(ib2,ib3) = T(ib2,ib1);
+	TM = HV(ib1,ib1).*ky_ikxy - HV(ib2,ib1).*kx_ikxy;
+	T(ib2,ib1) = T(ib2,ib1) + TM;
+	T(ib2,ib3) = T(ib2,ib3) - TM;
+	T(ib4,ib1) = -T(ib2,ib3);
+	T(ib4,ib3) = -T(ib2,ib1);
+
+	T(ib2,ib2) = -EV(ib1,ib2).*ekx_ikxyz - EV(ib2,ib2).*eky_ikxyz;
+	T(ib2,ib4) = T(ib2,ib2);
+	TM = HV(ib1,ib2).*ky_ikxy - HV(ib2,ib2).*kx_ikxy;
+	T(ib2,ib2) = T(ib2,ib2) + TM;
+	T(ib2,ib4) = T(ib2,ib4) - TM;
+	T(ib4,ib2) = -T(ib2,ib4);
+	T(ib4,ib4) = -T(ib2,ib2);
+end
+
+

fmmtd_efficiency.m

@@ -0,0 +1,63 @@
+%{
+Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
+
+This file is part of GratingFMM.
+
+GratingFMM is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+GratingFMM is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
+%}
+%% description:
+% calculate a matrix of diffraction efficiencies in case of the
+% diffraction by 2D gratings being periodic in x and y dimensions
+%% input:
+% xno, yno: numbers of Fourier harmonics in x and y dimensions
+% V_inc: incident field amplitude matrix of size (2*no,2)
+% V_dif: diffracted field amplitude matrix of size (2*no,2)
+% kx0, ky0: incident plane wave wavevector x and y projections (Bloch wavevector projections)
+% kgx, kgy: wavelength-to-period ratios (grating vectors)
+% eps1, eps2: substrate and superstrate permittivities
+%% output:
+% V_eff: efficiency matrix of size (2*no,2) if the if the incident field has
+%  propagating harmonics, otherwise (if the incident field is purely evanescent)
+%  the matrix of partial powers carried by each diffraction order
+% first index of V_inc, V_dif, V_eff indicates diffraction harmonics
+%   with indices 1:no being TE orders and no+1:2*no being TM orders
+%  (0-th order index is ind_0 = (ceil(xno/2)-1)*yno+ceil(yno/2))
+% second index of V_inc, V_dif, V_eff indicates whether the diffraction orders
+%  are in the substrate (V(:,1)) or in the superstrate (V(:,2))
+%% implementation
+function [V_eff] = fmmtd_efficiency(xno, yno, V_inc, V_dif, kx0, ky0, kgx, kgy, eps1, eps2)
+	no = xno*yno;
+	ib1 = 1:no; ib2 = no+1:2*no;
+
+	[kz1, kz2] = fmmtd_kxyz(xno, yno, kx0, ky0, kgx, kgy, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+
+	P_inc = sum( abs(V_inc(ib1,1).^2).*real(kz1) + abs(V_inc(ib1,2).^2).*real(kz2) ) ...
+				+ sum( abs(V_inc(ib2,1).^2).*real(kz1/eps1) + abs(V_inc(ib2,2).^2).*real(kz2/eps2) );
+	V_eff = zeros(2*no,2);
+	V_eff(ib1,1) = abs(V_dif(ib1,1).^2).*real(kz1);
+	V_eff(ib1,2) = abs(V_dif(ib1,2).^2).*real(kz2);
+	V_eff(ib2,1) = abs(V_dif(ib2,1).^2).*real(kz1/eps1);
+	V_eff(ib2,2) = abs(V_dif(ib2,2).^2).*real(kz2/eps2);
+
+	if abs(P_inc) > 1e-15
+		V_eff = (1/P_inc)*V_eff;
+	else
+		V_eff = 0.5*V_eff;
+	end
+end
+%
+% END
+%

fmmtd_kxyz.m

@@ -0,0 +1,53 @@
+%{
+Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
+
+This file is part of GratingFMM.
+
+GratingFMM is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+GratingFMM is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
+%}
+%% description:
+% calculate grating vectors for 2D grating (periodicity in x and y directions)
+% and plane wave propagation constants in homogeneous media below and 
+% above the grating
+%% input:
+% xno, yno: numbers of Fourier harmonics in X and Y dimensions
+% kx0, ky0: wavevector projections of an incident plane wave
+% kgx, kgy: wavelength-to-period ratios
+% eps1: permittivity of a medium below the grating (substrate)
+% eps2: permittivity of a medium above the grating (superstrate)
+%% output:
+% kz1: row of propagation constants in the substrate
+% kz2: row of propagation constants in the superstrate
+% kx: row of grating vectors in x-direction
+% ky: row of grating vectors in y-direction
+% kxy: row of grating vectors in xy plane
+%% implementation:
+function [kz1, kz2, kx, ky, kxy] = fmmtd_kxyz(xno, yno, kx0, ky0, kgx, kgy, eps1, eps2)
+	[kx,ky] = meshgrid(kx0 + kgx*(linspace(1,xno,xno) - ceil(xno/2)), ...
+											 ky0 + kgy*(linspace(1,yno,yno) - ceil(yno/2)));
+	kx = reshape(kx,1,[]);
+	ky = reshape(ky,1,[]);
+	kxy = kx.^2 + ky.^2;
+	
+	kz1 = sqrt(eps1 - kxy);
+	kz2 = sqrt(eps2 - kxy);
+	ind = angle(kz1) < -1e-12;
+	kz1(ind) = -kz1(ind);
+	ind = angle(kz2) < -1e-12;
+	kz2(ind) = -kz2(ind);
+
+	kxy = sqrt(kxy);
+end
+
+