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

fmm.m

@@ -0,0 +1,94 @@
+%{
+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 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)
+% 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
+% pol: polarization, either "TE" or "TM"
+%% output:
+% SM: scattering matrix of size (no,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
+% central harmonic index is ind_0 = ceil(no/2)
+% for example, an ampitude of the transmitted wave to i-th diffraction order
+%  from the substrate to the superstrate under the plane wave illumination
+%  with unit amplitude is SM(ind_0+i, ind_0, 2, 1)
+%% implementation
+function [SM] = fmm(no, kx0, kg, kh, eps1, eps2, FE, pol)
+		% wavevector projections
+	[kz1, kz2, kx] = fmm_kxz(no, kx0, 0, kg, eps1, eps2);
+		% block indices
+	ib1 = 1:no;
+	ib2 = no+1:2*no;
+	  % solve the eigenvalue problem:
+	ME = toeplitz(FE(no:2*no-1,1),FE(no:-1:1,1)); % permittivity Toeplitz matrix	
+	if (strcmp(pol,'TE')) % if TE polarization
+		ME(1:no+1:end) = ME(1:no+1:end) - (kx.^2);
+		[EV,MB] = eig(ME);
+		beta = transpose(sqrt(diag(MB)));
+			% check the branch of the square root
+		ind = angle(beta) < -1e-7;
+		beta(ind) = -beta(ind);
+			% eigen vector of the magnetic field
+		HV = -EV.*beta; % Hx
+	else % if TM polarization
+		MU = eye(no) / toeplitz(FE(no:2*no-1,2),FE(no:-1:1,2)); % inverce permittivity Toeplitz matrix
+		ME = -(diag(kx) / ME).*kx;
+		ME(1:no+1:end) = ME(1:no+1:end) + 1;
+		[EV,MB] = eig(ME*MU);
+		beta = transpose(sqrt(diag(MB)));
+			% check the branch of the square root
+		ind = angle(beta) < -1e-7;
+		beta(ind) = -beta(ind);
+			% eigen vector of the magnetic field
+		HV = (MU*EV)./beta;
+	end
+
+	bexp = exp((1i*kh)*beta);
+	  % apply the boundary conditions:
+	TS = fmm_calc_T(no,EV,HV,kz1,eps1,pol); % susbtrate-grating T-matrix
+	TC = fmm_calc_T(no,EV,HV,kz2,eps2,pol); % grating-cover T-matrix
+		% combine T-matrices
+	M1 = zeros(2*no,2*no);
+	M2 = zeros(2*no,2*no);
+	M1(ib1,ib1) = TS(ib2,ib1);
+	M1(ib1,ib2) = TS(ib2,ib2).*bexp;
+	M1(ib2,ib1) = TC(ib1,ib1).*bexp;
+	M1(ib2,ib2) = TC(ib1,ib2);
+	M2(ib1,ib1) = TS(ib1,ib1);
+	M2(ib1,ib2) = TS(ib1,ib2).*bexp;
+	M2(ib2,ib1) = TC(ib2,ib1).*bexp;
+	M2(ib2,ib2) = TC(ib2,ib2);
+
+	SM = M2S(M1/M2);
+end
+
+
+
+

fmm_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 collinear diffraction by 1D gratings
+% the returned value should be close to zero for pure dielectric structures
+% with no losses
+%% input:
+% no: number of Fourier harmonics
+% V_inc: incident field amplitude matrix of size (no, 2)
+% V_dif: diffracted field amplitude matrix of size (no, 2)
+%  first index of V_inc, V_dif indicates diffraction harmonics
+%   (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)
+% kg: wavelength-to-period ratio (grating vector)
+% eps1, eps2: substrate and superstrate permittivities
+% pol: polarization (either "TE" or "TM")
+%% 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] = fmm_balance(no, V_inc, V_dif, kx0, kg, eps1, eps2, pol)
+	[kz1, kz2] = fmm_kxz(no, kx0, 0, kg, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+
+	if (strcmp(pol,'TM'))
+		kz1 = kz1/eps1;
+		kz2 = kz2/eps2;
+	end
+	P_inc = sum( abs(V_inc(:,1).^2).*real(kz1) + abs(V_inc(:,2).^2).*real(kz2) );
+	P_dif = sum( abs(V_dif(:,1).^2).*real(kz1) + abs(V_dif(:,2).^2).*real(kz2) );
+
+	if (abs(P_inc) > 1e-15)
+		b = abs(P_dif/P_inc-1);
+	else
+		b = 0.5*P_dif;
+	end
+end
+%
+% END
+%

fmm_calc_T.m

@@ -0,0 +1,52 @@
+%{
+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 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
+% kz: row of plane wave propagation constants in the homogeneous medium
+% eps: permittivity the homogeneous medium
+% pol: polarization
+%% output:
+% T: T-matrix of size 4*no by 4*no
+%% implementation:
+function [T] = fmm_calc_T(no, EV, HV, kz, eps, pol)
+    T = zeros(2*no,2*no);
+		ib1 = [1:no];
+		ib2 = [no+1:2*no];
+    if strcmp(pol,'TE')
+			ikz = transpose(0.5./kz);
+			T(ib1,ib1) = HV.*ikz;
+			T(ib2,ib1) = 0.5*EV - T(ib1,ib1);
+			T(ib1,ib1) = T(ib1,ib1) + 0.5*EV;
+			T(ib1,ib2) = T(ib2,ib1);
+			T(ib2,ib2) = T(ib1,ib1);
+    elseif strcmp(pol,'TM')
+			eikz = transpose((0.5*eps)./kz);
+			T(ib1,ib1) = EV.*eikz;
+			T(ib2,ib1) = 0.5*HV - T(ib1,ib1);
+			T(ib1,ib1) = T(ib1,ib1) + 0.5*HV;
+			T(ib1,ib2) = -T(ib2,ib1);
+			T(ib2,ib2) = -T(ib1,ib1);
+    end
+end
+

fmm_efficiency.m

@@ -0,0 +1,59 @@
+%{
+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
+% collinear diffraction by 1D gratings
+%% input:
+% no: number of Fourier harmonics
+% V_inc: incident field amplitude matrix of size (no,2)
+% V_dif: diffracted field amplitude matrix of size (no,2)
+% kx0: incident plane wave wavevector x-projection (Bloch wavevector)
+% kg: wavelength-to-period ratio (grating vector)
+% eps1, eps2: substrate and superstrate permittivities
+% pol: polarization (either "TE" or "TM")
+%% output:
+% V_eff: efficiency matrix of size (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
+%  (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] = fmm_efficiency(no, V_inc, V_dif, kx0, kg, eps1, eps2, pol)
+	[kz1, kz2] = fmm_kxz(no, kx0, 0, kg, eps1, eps2);
+	kz1 = transpose(kz1);
+	kz2 = transpose(kz2);
+	if (strcmp(pol,'TM'))
+		kz1 = kz1/eps1;
+		kz2 = kz2/eps2;
+	end
+
+	V_eff = zeros(no,2);
+	P_inc = sum( abs(V_inc(:,1).^2).*real(kz1) + abs(V_inc(:,2).^2).*real(kz2) );
+	V_eff(:,1) = abs(V_dif(:,1).^2).*real(kz1);
+	V_eff(:,2) = abs(V_dif(:,2).^2).*real(kz2);
+
+	if (abs(P_inc) > 1e-15)
+		V_eff = V_eff/P_inc;
+	end
+end
+%
+% END
+%

fmm_kxz.m

@@ -0,0 +1,51 @@
+%{
+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 1D grating (periodicity in x direction)
+% and plane wave propagation constants in homogeneous media below and 
+% above the grating (substrate and superstrate)
+%% input:
+% no: number of Fourier harmonics
+% kx0, ky0: wavevector projections of an incident plane wave
+% kg: wavelength-to-period ratio
+% 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
+% kxy: row of grating vectors in xy plane
+%% implementation:
+function [kz1, kz2, kx, kxy] = fmm_kxz(no, kx0, ky0, kg, eps1, eps2)
+	ind = linspace(1,no,no);
+	kx = kx0 + kg*(ind - ceil(no/2));
+	kxy = kx.^2 + ky0^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
+%
+% end of fmm_kxz
+%