Физика \ Оптика фотонных кристаллов

Оптика фотонных кристаллов: Методические указания к практическим занятиям, страница 13

В приводимой программе изначально вводятся относительные величины, так как на первых парах учет размерностей загромождает программу излишними перемеными и затрудняет понимание процесса. Даже множитель  указывается только в одном месте при расчете фазы.

При желании, можно ввести размерности самостоятельно – это технически несложная работа, однако трудоемкая и обильная ловушками. Введение размерностей позволит лучше почувствовать экспериментальную ситуацию.

Переключатель xName допускает строковые значения Friq, Wave, Angle и Space, переключающие программу в режимы расчета спектров в зависимости от частоты, длины волны, угла, или в режим расчета распределения поля по глубине среды.

Программа расчета в пакете SciLab

1.  clear all; % clears all previous variables

2.  tic; % Simulation time started

3.   

4.  % Method, variabe range

5.  xName = 'Friq'; % Friq Wave Angle Space

6.  switch xName

7.      case 'Friq'

8.          Nx = 501; Rx = [0 3]; .5+.005*[-1 1]; 

9.      case 'Wave'

10.          Nx = 201; Rx = [.1 1];

11.      case 'Space'

12.          Rx = .5+.0025; Nx = 1; Rx = [Rx Rx];

13.      case 'Angle'

14.          Nx = 501; Rx = [-85 85]; vw0 = 1;

15.  end;

16.  xg = linspace(Rx(1),Rx(2),Nx)';

17.   

18.  % Angle of incoming field

19.  Theta = 0; %Theta = Theta_degree*pi/180;

20.  % Crystall

21.  N_sublayers = 30;

22.  if ~strcmp(xName,'Space'), N_sublayers = 1; end;

23.  layer = ...  layer characteristics (n, l, grid)

24.      [1 1 N_sublayers;... air around the structure

25.      2 1/3 N_sublayers;... layer 1

26.      1 2/3 N_sublayers;... layer 2

27.      1 4/3 N_sublayers... defect layer

28.      ];

29.  % lattice 1: 123232323432323232343232323231

30.  % lattice 2: 123234323232323232323234323231

31.  latt = [1,repmat([2,3],1,4),2,4,...

32.      repmat([2,3],1,4),2,4,...

33.      repmat([2,3],1,4),2,1];

34.  % latt = [1,repmat([2,3],1,2),2,4,...

35.  %     repmat([2,3],1,8),2,4,...

36.  %     repmat([2,3],1,2),2,1];

37.   

38.  % build z and n grids

39.  zg = 0; ng = [];

40.  for r = 1:length(latt)

41.      zg = [zg zg(end)+linspace(0,layer(latt(r),2),layer(latt(r),3)+1)]; %#ok<AGROW>

42.      ng = [ng layer(latt(r),1)*ones(1,layer(latt(r),3)+1)]; %#ok<AGROW>

43.  end;

44.  zg = zg(2:end); Nz = length(zg); dzg = diff(zg);

45.  n0 = layer(latt(1),1);

46.  sinTheta = n0*sind(Theta); % Theta>0 - TEM, <0 - TMM

47.  L_opt = layer(2,1)*layer(2,2)+...

48.      layer(3,1)*layer(3,2);

49.  % (d1*n1+d2*n2); % optical dencity

50.   

51.  % Fields

52.  E = zeros(2,Nz);

53.  E(:,Nz) = [1 0]; % no incoming field from the right side

54.   

55.  % Spectr

56.  tg = zeros(1,Nx); % transmission coefficient

57.  for rx = 1:Nx,

58.      ncos = sqrt(ng.^2-sinTheta^2);

59.      switch xName

60.          case 'Friq'

61.              wv = xg(rx); % mkm^-1, w/c = 2*pi/lambda -- wave vector

62.          case 'Wave'

63.              wv = 1/xg(rx);

64.          case 'Space'

65.              wv = xg;

66.          case 'Angle'

67.              wv = vw0; sinTheta = n0*sind(xg(rx));

68.              ncos = sqrt(ng.^2-sinTheta^2);

69.      end;

70.      if sinTheta>=0,

71.          mu=ncos(2:Nz)./ncos(1:Nz-1);

72.      else

73.          mu=ncos(1:Nz-1)./ncos(2:Nz).*(ng(2:Nz)./ng(1:Nz-1)).^2;

74.      end;

75.      % Transfer matrix cycle

76.      for rz = Nz-1:-1:1,

77.          if zg(rz)==zg(rz+1),

78.              E(:,rz) = E(:,rz+1) + (mu(rz)-1)*(E(1,rz+1)-E(2,rz+1))*[.5;-.5];

79.          else

80.              Rot = exp(-i*2*pi/L_opt*ncos(rz)*wv*dzg(rz));

81.              E(:,rz) = [E(1,rz+1)/Rot;E(2,rz+1)*Rot];

82.          end;

83.      end;

84.      % Transfer matrix finished

85.      % tg(rx) = 1-abs(E(2,1)/E(1,1)).^2; % only when it is no absorbtion

86.      tg(rx) = abs(E(1,end)/E(1,1)).^2; % only for n_left = n_right

87.  end;

88.   

89.  % Sets figure to draw

90.  hf = findobj('Name',xName);

91.  if isempty(hf), hf = figure; hold on;

92.      set(hf,'Name',xName,'WindowStyle','Dock');

93.  else

94.      figure(hf); h_child = get(gca,'children');

95.      n_child = length(h_child); colors = hsv(n_child);

96.      for r=1:n_child, set(h_child(r),'Color',colors(r,:)); end;

97.  end;

98.   

99.  % DRAW

100.  if ~strcmp(xName,'Space'),

101.      plot(xg,tg,'k.-');

102.      axis ([Rx+.05*(Rx(2)-Rx(1))*[-1 1] [-.05 1.05]])

103.  else % Spatial distribution of field

104.      delete(get(gcf,'children'));

105.      plot(zg,real(ng),'g.--'); hold on;

106.      Intensity = abs(sum(E)).^2;

107.      axis ([max(zg)*[-.05 1.05] 2*[-.05 1.05]])

108.      % Intensity = real(sum(E));

109.      % axis ([max(zg)*[-.05 1.05] 2*[-.55 1.05]])

110.      plot(zg,Intensity/max(Intensity),'r-');

111.      title(['Transmission = ' num2str(tg)]);

112.  end;

113.  toc; % elapsed time

Результаты расчета и решение

а) Методом трансфер-матрицы провести расчет спектра пропускания фотонного кристалла в случае нормального падения света на образец.

На рис. 2. изображен численно построенный спектр исследуемого фотонного кристалла.  измеряется в единицах , где  – оптическая толщина. В данном случае, . В вакууме соответствующая волна будет иметь длину 200 нм, равную оптической толщине среды. Видно, что , то есть такая волна проходит через ФК без отражения. А вот волна, у которой частота в два раза меньше и длина в два раза больше, полностью отражается средой, поскольку период среды для нее равен половине длины волны, а каждый слой – четвертьволновой.

Скачать файл