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

calc_emn_bin.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:
+% Fourier matrix of the permittivity and the inverse permittivity of a
+% 1D binary grating
+%% input:
+% no - number of Fourier harmonics
+% alps: relative widths of the grating period elements (should meet the conditions
+%   0 < alps(i) < 1 and sum(alps) == 1)
+% poss: relative positions of the centers of the grating period elements
+%   (should not overlap)
+% eps: permittivities of the grating period elements
+%% output:
+% FE: Fourier matrix of size (2*no-1,2), which contains
+% Fourier components of the permittivity (FE(:,1)) and the inverse
+% permittivity (FE(:,2))
+%% implementation
+function [FE] = calc_emn_bin(no, alps, poss, eps)
+	if (length(alps) ~= length(poss)) || (length(alps) ~= length(eps)) || (abs(sum(alps)-1) > 1e-12)
+		error("calc_emn_bin input error");
+	end
+
+	FE = zeros(2*no-1,2);
+	ind = transpose(linspace(1,no-1,no-1));
+
+	for k = 1:length(alps)
+		ifun = sin(ind*pi*alps(k))./(ind*pi);
+		te = exp(-(2*pi*1i*poss(k))*ind);
+			% zero harmonics:
+		FE(no,1) = FE(no,1) + eps(k)*alps(k);
+		FE(no,2) = FE(no,2) + alps(k)/eps(k);
+			% non-zero harmonics:
+		tmp = eps(k)*ifun;
+		FE(no+1:2*no-1,1) = FE(no+1:2*no-1,1) + tmp.*te;
+		FE(no-1:-1:1,1) = FE(no-1:-1:1,1) + tmp.*conj(te);
+		tmp = (1/eps(k))*ifun;
+		FE(no+1:2*no-1,2) = FE(no+1:2*no-1,2) + tmp.*te;
+		FE(no-1:-1:1,2) = FE(no-1:-1:1,2) + tmp.*conj(te);
+	end
+end
+%
+% end of calc_emn_bin
+%

calc_emn_lam.m

@@ -0,0 +1,49 @@
+%{
+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:
+% Fourier matrix of the permittivity and the inverse permittivity of a
+% 1D lamellar grating
+%% input:
+% no - number of Fourier harmonics
+% alp - lamellae width relative to the grating period (0 < alp < 1)
+% eps1 - lamellae permittvity
+% eps2 - surrounding medium permittvity
+%% output:
+% FE: Fourier matrix of size (2*no-1,2), which contains
+% Fourier components of the permittivity (FE(:,1)) and the inverse
+% permittivity (FE(:,2))
+%% implementation
+function [FE] = calc_emn_lam(no, alp, eps1, eps2)
+	FE = zeros(2*no-1,2);
+	te1 = eps1 - eps2;
+	te2 = 1/eps1 - 1/eps2;
+		% zero harmonic
+	FE(no,1) = eps1*alp + eps2*(1-alp);
+	FE(no,2) = alp/eps1 + (1-alp)/eps2;
+		% non-zero harmonics:
+	ind = linspace(1,no-1,no-1);
+	ifun = transpose(sin(ind*pi*alp)./(ind*pi));
+	FE(no+1:2*no-1,1) = te1*ifun;
+	FE(no+1:2*no-1,2) = te2*ifun;
+	FE(no-1:-1:1,1) = FE(no+1:2*no-1,1);
+	FE(no-1:-1:1,2) = FE(no+1:2*no-1,2);
+end
+%
+% end of calc_emn_lam
+%

calc_emntd_bin.m

@@ -0,0 +1,76 @@
+%{
+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 permittivity Fourier matrix of a 2D binary grating
+% being periodic in x and y dimensions of the 3D Cartesian coordinates
+%% input:
+% xno, yno: numbers of Fourier harmonics
+% cx, cy: rows of centers of a 2D rectangular mesh filling the grating period along
+%   x and y dimensions normalized by the period (each value should be between -0.5 and 0.5)
+% dx, dy: rows of widths of a 2D rectangular mesh elements along
+%   x and y dimensions normalized by the period (each value should be between 0 and 1)
+% eps: row of permittivities for each mesh element (length(eps) should be
+%   equal to length(cx)*length(cy))
+%% output:
+% FE: cell array containing two Fourier matrices of the permittivity and
+% inverse permittivity
+%% implementation:
+function [FE] = calc_emntd_bin(xno, yno, cx, cy, dx, dy, eps)
+	nx = length(cx);
+	ny = length(cy);
+	if (length(cx)~=length(dx)) || (length(cy)~=length(dy)) || (length(eps)~=(nx*ny))
+		error("incorrect binary grating definition");
+	end
+
+	FE = cellmat(1,2,2*yno-1,2*xno-1);
+
+	[CX,CY] = meshgrid(cx,cy);
+	[DX,DY] = meshgrid(dx,dy);
+
+	ix = linspace(1,xno-1,xno-1);
+	iy = linspace(1,yno-1,yno-1);
+	[IX,IY] = meshgrid(ix,iy);
+	
+	for ip = 1:nx*ny
+		fx = (sin(ix*pi*DX(ip))./(pi*ix)).*exp((-2*pi*1i*CX(ip))*ix);
+		fy = (sin(iy*pi*DY(ip))./(pi*iy)).*exp((-2*pi*1i*CY(ip))*iy);
+		FX = (sin(IX*pi*DX(ip))./(pi*IX)).*exp((-2*pi*1i*CX(ip))*IX);
+		FY = (sin(IY*pi*DY(ip))./(pi*IY)).*exp((-2*pi*1i*CY(ip))*IY);
+
+		M = zeros(2*yno-1,2*xno-1);
+
+		M(yno+1:2*yno-1,xno) = DX(ip)*fy;
+		M(yno-1:-1:1,xno) = conj(M(yno+1:2*yno-1,xno));
+		M(yno,xno+1:2*xno-1) = DY(ip)*fx;
+		M(yno,xno-1:-1:1) = conj(M(yno,xno+1:2*xno-1));
+
+		M(yno+1:2*yno-1,xno+1:2*xno-1) = FX.*FY;
+		M(yno+1:2*yno-1,xno-1:-1:1) = conj(FX).*FY;
+		M(yno-1:-1:1,xno+1:2*xno-1) = FX.*conj(FY);
+		M(yno-1:-1:1,xno-1:-1:1) = conj(FX.*FY);
+
+		FE{1,1} = FE{1,1} + eps(ip)*M;
+		FE{1,2} = FE{1,2} + (1/eps(ip))*M;
+		FE{1,1}(yno,xno) = FE{1,1}(yno,xno) + DX(ip)*DY(ip)*eps(ip);
+		FE{1,2}(yno,xno) = FE{1,2}(yno,xno) + DX(ip)*DY(ip)/eps(ip);
+	end
+end
+%
+% end of calc_emntd_cyl
+%

calc_emntd_cyl.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 permittivity Fourier matrix of a 2D grating with cylindrical
+% pithches being periodic in x and y dimensions of the 3D Cartesian coordinates
+%% input:
+% xno, yno: numbers of Fourier harmonics
+% rpx, rpy: radius-to-period ratios
+% eps_c: cylinder permittivity
+% eps_m: permittivity of a medium which surrounds the cylinder
+%% output:
+% FE: cell array containing two Fourier matrices of the permittivity and
+% inverse permittivity
+%% implementation:
+function FE = calc_emntd_cyl(xno, yno, rpx, rpy, eps_c, eps_m)
+	FE = cellmat(1,2,2*yno-1,2*xno-1);
+
+	ix = linspace(1,xno-1,xno-1);
+	iy = linspace(1,yno-1,yno-1);
+	[IX,IY] = meshgrid(ix,iy);
+	
+	fx = rpy*besselj(1,(2*pi*rpx)*ix) ./ ix;
+	fy = rpx*besselj(1,(2*pi*rpy)*iy) ./ iy;
+	FXY = sqrt((rpx*IX).^2 + (rpy*IY).^2);
+	FXY = (rpx*rpy)*besselj(1,(2*pi)*FXY) ./ FXY;
+
+	M = zeros(2*yno-1,2*xno-1);
+
+	M(yno,xno) = pi*rpx*rpy;
+
+	M(yno+1:2*yno-1,xno) = fy;
+	M(yno-1:-1:1,xno) = M(yno+1:2*yno-1,xno);
+	M(yno,xno+1:2*xno-1) = fx;
+	M(yno,xno-1:-1:1) = M(yno,xno+1:2*xno-1);
+
+	M(yno+1:2*yno-1,xno+1:2*xno-1) = FXY;
+	M(yno+1:2*yno-1,xno-1:-1:1) = FXY;
+	M(yno-1:-1:1,xno+1:2*xno-1) = FXY;
+	M(yno-1:-1:1,xno-1:-1:1) = FXY;
+
+	FE{1,1} = FE{1,1} + (eps_c - eps_m)*M;
+	FE{1,2} = FE{1,2} + (1/eps_c - 1/eps_m)*M;
+	FE{1,1}(yno,xno) = FE{1,1}(yno,xno) + eps_m;
+	FE{1,2}(yno,xno) = FE{1,2}(yno,xno) + 1/eps_m;
+end
+%
+% end of calc_emntd_cyl
+%

calc_emntd_lam.m

@@ -0,0 +1,65 @@
+%{
+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 permittivity Fourier matrix of a 2D lamellar grating
+% being periodic in x and y dimensions of the 3D Cartesian coordinates
+%% input:
+% xno, yno: numbers of Fourier harmonics
+% ax, ay: aspect ratios, should be between 0 and 1
+% eps1: pitch permittivity (occupies ax*ay volume fraction of the period)
+% eps2: permittivity of a medium which surrounds the pitch
+%   (occupies (1-ax)*(1-ay) volume fraction of the period)
+%% output:
+% FE: cell array containing two Fourier matrices of the permittivity and
+% inverse permittivity
+%% implementation:
+function [FE] = calc_emntd_lam(xno, yno, ax, ay, eps1, eps2)
+	FB = zeros(2*yno-1,2*xno-1);
+	FE = cell(1,2);
+	FE{1,1} = eps2 + FB;
+	FE{1,2} = 1/eps2 + FB;
+
+	ix = linspace(1,xno-1,xno-1);
+	iy = linspace(1,yno-1,yno-1);
+	[IX,IY] = meshgrid(ix,iy);
+		
+	fx = sin(ix*pi*ax)./(pi*ix);
+	fy = sin(iy*pi*ay)./(pi*iy);
+	FX = sin(IX*pi*ax)./(pi*IX);
+	FY = sin(IY*pi*ay)./(pi*IY);
+
+	FB(yno+1:2*yno-1,xno) = ax*fy;
+	FB(yno-1:-1:1,xno) = FB(yno+1:2*yno-1,xno);
+
+	FB(yno,xno+1:2*xno-1) = ay*fx;
+	FB(yno,xno-1:-1:1) = FB(yno,xno+1:2*xno-1);
+
+	FB(yno+1:2*yno-1,xno+1:2*xno-1) = FX.*FY;
+	FB(yno+1:2*yno-1,xno-1:-1:1) = FB(yno+1:2*yno-1,xno+1:2*xno-1);
+	FB(yno-1:-1:1,xno+1:2*xno-1) = FB(yno+1:2*yno-1,xno+1:2*xno-1);
+	FB(yno-1:-1:1,xno-1:-1:1) = FB(yno+1:2*yno-1,xno+1:2*xno-1);
+
+	FE{1,1} = (eps1 - eps2)*FB;
+	FE{1,2} = (1/eps1 - 1/eps2)*FB;
+	FE{1,1}(yno,xno) = eps2 + (eps1 - eps2)*ax*ay;
+	FE{1,2}(yno,xno) = 1/eps2 + (1/eps1 - 1/eps2)*ax*ay;
+end
+%
+% end of calc_emntd_lam
+%