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

fmmnc.m

@@ -0,0 +1,119 @@
+%{
+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 non-collinear diffraction by 1D gratings being periodic in x
+%  dimension of the Cartesian coordinates
+%% input:
+% no: number of Fourier harmonics
+% kx0: incident plane wave wavevector x-projection (Bloch wavevector)
+% ky0: incident plane wave wavevector y-projection
+%  (ky0=0 corresponds to the collinear diffraction)
+% kg: wavelength-to-period ratio (grating vector)
+% 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)
+% 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 if 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(no/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] = fmmnc(no, kx0, ky0, kg, kh, eps1, eps2, FE)
+		% block indices
+	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, kxy] = fmm_kxz(no, kx0, ky0, kg, eps1, eps2);
+
+		% permittivity Toeplitz matrices
+	ME = toeplitz(FE(no:2*no-1,1),FE(no:-1:1,1)); % permittivity Toeplitz matrix
+	MU = toeplitz(FE(no:2*no-1,2),FE(no:-1:1,2)); % inverse permittivity Toeplitz matrix
+
+		% initialize the eigenvectors
+	EV = zeros(2*no,2*no);
+	HV = zeros(2*no,2*no);
+
+		% matrix for the electric field
+	IMU = eye(no)/MU;
+	TM = ME\((kx').*IMU);
+	EV(ib1,ib1) = IMU - (kx').*TM - (ky0^2)*eye(no);
+	EV(ib2,ib1) = ky0*(diag(kx) - TM);
+	EV(ib2,ib2) = ME - diag(kxy.*kxy);
+
+		% solve the eigenvalue problem for the electric field
+	[EV,MB] = eig(EV);
+	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 eigen vectors
+	HV(ib2,ib2) = diag(ky0*kx);
+	HV(ib1,ib1) = -HV(ib2,ib2);
+	HV(ib2,ib1) = IMU - (ky0^2)*eye(no);
+	HV(ib1,ib2) = diag(kx.*kx) - ME;
+	HV = (HV*EV).*(1./beta);
+
+	bexp = exp((1i*kh)*beta);
+
+		% apply the boundary conditions:
+		% calculate T-matrices
+	TS = fmmnc_calc_T(no,EV,HV,kx,ky0,kxy,kz1,eps1); % susbtrate-grating T-matrix
+	TC = fmmnc_calc_T(no,EV,HV,kx,ky0,kxy,kz2,eps2); % grating-cover T-matrix
+		% initialization
+	M1 = zeros(4*no,4*no);
+	M2 = zeros(4*no,4*no);
+		% combine T-matrices
+	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);
+
+		% attain S-matrix
+	SM = M2S(M1/M2);
+end

fmmnc_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 non-collinear diffraction by 1D gratings being periodic in x dimension
+%% input:
+% no: number of Fourier harmonics
+% 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(no/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: incident plane wave wavevector x-projection (Bloch wavevector)
+% ky0: incident plane wave wavevector y-projection
+% kg: wavelength-to-period ratio (grating vector)
+% 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] = fmmnc_balance(no, V_inc, V_dif, kx0, ky0, kg, eps1, eps2)
+	[kz1, kz2] = fmm_kxz(no, kx0, ky0, kg, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+
+	ib1 = 1:no;
+	ib2 = no+1:2*no;
+
+	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
+%

fmmnc_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: 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
+% ky0: Bloch wavevector y-projection
+% 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] = fmmnc_calc_T(no, EV, HV, kx, ky0, kxy, kz, eps)
+	T = zeros(4*no,4*no);
+
+		% pre-calculate combinations of wavevector projections
+	ikxy = 0.5./kxy;
+	ky_ikxy = ky0*ikxy;
+	kx_ikxy = kx.*ikxy;
+	ikxyz = ikxy./kz;
+	ky_ikxyz = transpose(ky0*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

fmmnc_efficiency.m

@@ -0,0 +1,64 @@
+%{
+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
+% non-collinear diffraction by 1D gratings being periodic in x dimension
+%% input:
+% no: number of Fourier harmonics
+% V_inc: incident field amplitude matrix of size (2*no,2)
+% V_dif: diffracted field amplitude matrix of size (2*no,2)
+% kx0: incident plane wave wavevector x-projection (Bloch wavevector)
+% ky0: incident plane wave wavevector y-projection
+% kg: wavelength-to-period ratio (grating vector)
+% 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(no/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] = fmmnc_efficiency(no, V_inc, V_dif, kx0, ky0, kg, eps1, eps2)
+	[kz1, kz2] = fmm_kxz(no, kx0, ky0, kg, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+		
+	ib1 = 1:no;
+	ib2 = no+1:2*no;
+
+	V_eff = zeros(2*no,2);
+	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(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 = V_eff/P_inc;
+	else
+		V_eff = 0.5*V_eff;
+	end
+end
+%
+% END
+%