% 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=block3('360e-9','1400e-9','220e-9','base','center','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');
g2=block3('190e-9','990e-9','220e-9','base','center','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1,g2};
s.name={'BLK1','BLK2'};
s.tags={'g1','g2'};

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 = 'ElectromagneticWaves';
appl.module = 'RF';
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_rfw';
clear prop
prop.analysis='eigen';
appl.prop = prop;
clear bnd
bnd.type = {'E0','cont','periodic'};
bnd.index = {0,0,1};
bnd.kper = {{0;0;0},{0;0;0},{'kx';0;0}};
bnd.pertype = {'sym','sym','floque'};
bnd.ind = [3,1,1,1,1,2,2,1,1,2,2,3];
appl.bnd = bnd;
clear equ
equ.n = {3.493,1};
equ.shape = [1;2];
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;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);
kk1=0;
n=1;
bands=4;
for kk1=0:0.02:1
% Global expressions
fem.globalexpr = {'a','360E-9[m]', ...
  'b','pi/a', ...
  'k1',kk1, ...
  'kx','(k1+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',{'tExEyEz10','tExEyEz20','tExEyEz21'}, ...
               'outcomp',{'tExEyEz10','tExEyEz20','tExEyEz21'}, ...
               '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');