demo_fmmtd.m 2.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970
  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. %% demonstration script for the 2D grating Fourier Modal Method calculations
  16. clc;
  17. %format long;
  18. %% initialization
  19. wl = 1; % wavelength in micrometers
  20. phi = 0;
  21. theta = 0.001; % for normal incidence take a small non-zero value
  22. kx0 = sin(theta*pi/180)*cos(phi*pi/180); % incidence wavevector horizontal projection (dimensionless) (Bloch wavevector)
  23. ky0 = sin(theta*pi/180)*sin(phi*pi/180); % incidence wavevector horizontal projection (dimensionless) (Bloch wavevector)
  24. gpx = 0.72; % grating period in x dimension
  25. gpy = 0.72; % grating period it y dimension
  26. gh = 0.5; % grating depth
  27. wv = 2*pi/wl; % wavevector
  28. % dimensionless variables
  29. kgx = wl/gpx;
  30. kgy = wl/gpy;
  31. kh = wv*gh;
  32. xno = 15; % number of Fourier modes
  33. yno = 15; % number of Fourier modes
  34. no = xno*yno;
  35. ixy = (ceil(xno/2)-1)*yno+ceil(yno/2);
  36. eps_sub = 1.5; % substrate permittivity
  37. eps_gr = 3.17^2; % get_epsAu_Drude(wl); %grating permittivity
  38. eps_sup = 1; % superstrate permittivity
  39. %% S-matrix calculation
  40. % calculate Fourier image matrix of the dielectric permittivity function
  41. % for a 2D lamellar grating with filling factors 0.5,0.5
  42. FE = calc_emntd_lam(xno,yno,0.5,0.5,eps_gr,eps_sup);
  43. % scattering matrix of the grating
  44. SM = fmmtd(xno,yno,kx0,ky0,kgx,kgy,kh,eps_sub,eps_sup,FE);
  45. %% diffraction of a plane wave example
  46. % incident field amplitude vector
  47. V_inc = zeros(2*no,2);
  48. V_inc((ceil(xno/2)-1)*yno+ceil(yno/2),2) = 1; % TE polarized plane wave (0-th harmonic) coming from the superstrate
  49. % define diffracted amplitude vector
  50. V_dif = zeros(2*no,2);
  51. % calculate diffraction vector
  52. V_dif(:,1) = SM(:,:,1,1)*V_inc(:,1) + SM(:,:,1,2)*V_inc(:,2); % diffraction to the substrate
  53. V_dif(:,2) = SM(:,:,2,1)*V_inc(:,1) + SM(:,:,2,2)*V_inc(:,2); % diffraction to the superstrate
  54. V_eff = fmmtd_efficiency(xno,yno,V_inc,V_dif,kx0,ky0,kgx,kgy,eps_sub,eps_sup);
  55. disp("efficiency:");
  56. disp(V_eff(ixy,1)); % 0th order power transmission coefficient
  57. disp(V_eff(ixy,2)); % 0th order power reflection coefficient
  58. % calculate the power balance
  59. b = fmmtd_balance(xno,yno,V_inc,V_dif,kx0,ky0,kgx,kgy,eps_sub,eps_sup);
  60. disp("balance:");
  61. disp(b);