fmmtd_calc_T.m 3.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687
  1. %{
  2. Copyright © 2020 Alexey A. Shcherbakov. All rights reserved.
  3. This file is part of GratingFMM.
  4. GratingFMM is free software: you can redistribute it and/or modify
  5. it under the terms of the GNU General Public License as published by
  6. the Free Software Foundation, either version 2 of the License, or
  7. (at your option) any later version.
  8. GratingFMM is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11. GNU General Public License for more details.
  12. You should have received a copy of the GNU General Public License
  13. along with GratingFMM. If not, see <https://www.gnu.org/licenses/>.
  14. %}
  15. %% description:
  16. % T-matrix of an interface between a photonic crystal and a homogeneous medium
  17. % the case of the non-collinear diffraction on 1D gratings
  18. %% input:
  19. % no: total number of Fourier harmonics
  20. % EV: matrix of the electric field eigenvectors for the photonic crystal
  21. % HV: matrix of the magnetic field eigenvectors for the photonic crystal
  22. % kx: row of grating vector projections in x direction
  23. % ky: row of grating vector projections in y direction
  24. % kxy: row of grating vector projections in xy plane
  25. % kz: row of plane wave propagation constants in the homogeneous medium
  26. % eps: permittivity the homogeneous medium
  27. %% output:
  28. % T: T-matrix of size 4*no by 4*no
  29. %% implementation:
  30. function [T] = fmmtd_calc_T(no, EV, HV, kx, ky, kxy, kz, eps)
  31. T = zeros(4*no,4*no);
  32. % pre-calculate combinations of wavevector projections
  33. ikxy = 0.5./kxy;
  34. ky_ikxy = ky.*ikxy;
  35. kx_ikxy = kx.*ikxy;
  36. ikxyz = ikxy./kz;
  37. ky_ikxyz = transpose(ky.*ikxyz);
  38. kx_ikxyz = transpose(kx.*ikxyz);
  39. eky_ikxyz = eps*ky_ikxyz;
  40. ekx_ikxyz = eps*kx_ikxyz;
  41. kx_ikxy = transpose(kx_ikxy);
  42. ky_ikxy = transpose(ky_ikxy);
  43. % block indices
  44. ib1 = 1:no;
  45. ib2 = no+1:2*no;
  46. ib3 = 2*no+1:3*no;
  47. ib4 = 3*no+1:4*no;
  48. % fill the T-matrix blocks
  49. T(ib1,ib1) = EV(ib1,ib1).*ky_ikxy - EV(ib2,ib1).*kx_ikxy;
  50. T(ib1,ib3) = T(ib1,ib1);
  51. TM = HV(ib1,ib1).*kx_ikxyz + HV(ib2,ib1).*ky_ikxyz;
  52. T(ib1,ib1) = T(ib1,ib1) + TM;
  53. T(ib1,ib3) = T(ib1,ib3) - TM;
  54. T(ib3,ib3) = T(ib1,ib1);
  55. T(ib3,ib1) = T(ib1,ib3);
  56. T(ib1,ib2) = EV(ib1,ib2).*ky_ikxy - EV(ib2,ib2).*kx_ikxy;
  57. T(ib1,ib4) = T(ib1,ib2);
  58. TM = HV(ib1,ib2).*kx_ikxyz + HV(ib2,ib2).*ky_ikxyz;
  59. T(ib1,ib2) = T(ib1,ib2) + TM;
  60. T(ib1,ib4) = T(ib1,ib4) - TM;
  61. T(ib3,ib2) = T(ib1,ib4);
  62. T(ib3,ib4) = T(ib1,ib2);
  63. T(ib2,ib1) = -EV(ib1,ib1).*ekx_ikxyz - EV(ib2,ib1).*eky_ikxyz;
  64. T(ib2,ib3) = T(ib2,ib1);
  65. TM = HV(ib1,ib1).*ky_ikxy - HV(ib2,ib1).*kx_ikxy;
  66. T(ib2,ib1) = T(ib2,ib1) + TM;
  67. T(ib2,ib3) = T(ib2,ib3) - TM;
  68. T(ib4,ib1) = -T(ib2,ib3);
  69. T(ib4,ib3) = -T(ib2,ib1);
  70. T(ib2,ib2) = -EV(ib1,ib2).*ekx_ikxyz - EV(ib2,ib2).*eky_ikxyz;
  71. T(ib2,ib4) = T(ib2,ib2);
  72. TM = HV(ib1,ib2).*ky_ikxy - HV(ib2,ib2).*kx_ikxy;
  73. T(ib2,ib2) = T(ib2,ib2) + TM;
  74. T(ib2,ib4) = T(ib2,ib4) - TM;
  75. T(ib4,ib2) = -T(ib2,ib4);
  76. T(ib4,ib4) = -T(ib2,ib2);
  77. end