123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687 |
- %{
- 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
|