% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.4';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 248;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2007/10/10 16:07:51 $';
fem.version = vrsn;

% Geometry
g1=rect2(0.8,0.4,'base','center','pos',[0,0]);
g2=rect2('360E-9','1400E-9','base','center','pos',{'0','0'},'rot','0');
g3=rect2(8.0E-7,4.0E-7,'base','center','pos',[0,0]);
g4=rect2('190E-9','990E-9','base','center','pos',{'0','0'},'rot','0');

% Analyzed geometry
clear s
s.objs={g2,g4};
s.name={'R1','E1'};
s.tags={'g2','g4'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'InPlaneWaves';
appl.module = 'RF';
appl.assignsuffix = '_rfwe';
clear prop
prop.analysis='eigen';
appl.prop = prop;
clear bnd
bnd.pertype = {'sym','sym','floque'};
bnd.type = {'periodic','cont','periodic'};
bnd.index = {2,0,0};
bnd.kper = {{0;0},{0;0},{'kx';'ky'}};
bnd.ind = [3,1,1,2,2,2,2,3];
appl.bnd = bnd;
clear equ
equ.n = {{1;1;3.493},1};
equ.matparams = 'n';
equ.ind = [1,2];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;
kk1=0;
kk2=0;
kk=0;
n=1;
bands=6;
for kk1=0:0.02:1
% Global expressions
fem.globalexpr = {'a','360E-9[m]', ...
  'b','pi/a', ...
  'k1',kk1,...
  'k2',kk2,...
  'k',kk,...
  'kx','(k1+k)*b', ...
  'ky','(k2+k)*b', ...
  };
% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femeig(fem, ...
               'complexfun','on', ...
               'conjugate','on', ...
               'solcomp',{'Ez'}, ...
               'outcomp',{'Ez'}, ...
               'neigs',bands);

% Save current fem structure for restart purposes
fem0=fem;
i=1;
while i<=bands
f(i,n)=-imag(fem.sol.lambda(i)/pi/2);
i=i+1;
end
n=n+1;
end
i=1;
while i<=bands
    plot(0:0.02:1,f(i,:)/1e12);hold on; grid on;
i=i+1;
end
hold on;
plot(0:0.5*250*24/1e4:250*24/1e4,0:250*0.5:250,'r');