calc_emn_bin.m 2.1 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859
  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. % Fourier matrix of the permittivity and the inverse permittivity of a
  17. % 1D binary grating
  18. %% input:
  19. % no - number of Fourier harmonics
  20. % alps: relative widths of the grating period elements (should meet the conditions
  21. % 0 < alps(i) < 1 and sum(alps) == 1)
  22. % poss: relative positions of the centers of the grating period elements
  23. % (should not overlap)
  24. % eps: permittivities of the grating period elements
  25. %% output:
  26. % FE: Fourier matrix of size (2*no-1,2), which contains
  27. % Fourier components of the permittivity (FE(:,1)) and the inverse
  28. % permittivity (FE(:,2))
  29. %% implementation
  30. function [FE] = calc_emn_bin(no, alps, poss, eps)
  31. if (length(alps) ~= length(poss)) || (length(alps) ~= length(eps)) || (abs(sum(alps)-1) > 1e-12)
  32. error("calc_emn_bin input error");
  33. end
  34. FE = zeros(2*no-1,2);
  35. ind = transpose(linspace(1,no-1,no-1));
  36. for k = 1:length(alps)
  37. ifun = sin(ind*pi*alps(k))./(ind*pi);
  38. te = exp(-(2*pi*1i*poss(k))*ind);
  39. % zero harmonics:
  40. FE(no,1) = FE(no,1) + eps(k)*alps(k);
  41. FE(no,2) = FE(no,2) + alps(k)/eps(k);
  42. % non-zero harmonics:
  43. tmp = eps(k)*ifun;
  44. FE(no+1:2*no-1,1) = FE(no+1:2*no-1,1) + tmp.*te;
  45. FE(no-1:-1:1,1) = FE(no-1:-1:1,1) + tmp.*conj(te);
  46. tmp = (1/eps(k))*ifun;
  47. FE(no+1:2*no-1,2) = FE(no+1:2*no-1,2) + tmp.*te;
  48. FE(no-1:-1:1,2) = FE(no-1:-1:1,2) + tmp.*conj(te);
  49. end
  50. end
  51. %
  52. % end of calc_emn_bin
  53. %