% that is general plasmons calculations, also there is comparision with
% rigorous theory

clear all


import com.comsol.model.*
import com.comsol.model.util.*
global model;
model = ModelUtil.create(['Model' num2str(fix(100*rand))]);
model.name('plasmonic waveguide mode analysis');

disp('started')

geom1 = model.geom.create('geom1', 2);
model.geom('geom1').lengthUnit([native2unicode(hex2dec('00b5'), 'Cp1252') 'm']);  % set length unit in um

MyGeometry=geom1; DoSetMaterials=true;DoSetMesh=true;MyName00='A';
MyGeometryName='geom1';
%%%%%%%%%%%%%% geometry

newMaterial=false;
%PlasmonGeometry1  %metal strip with dielectric inset
%PlasmonGeometry3  %metal strip  on double dielectric 
PlasmonGeometry3a;newMaterial=true;    %metalic strip on single dielectric
%PlasmonGeometry3b   %metalic strip on triple dielectric
%PlasmonGeometry4     % metalic strp on dielectric with inset. Difference with Geometry 1?

disp('geometry OK')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
geom1.run;

% set materials
if newMaterial
MyNewMaterial=setMaterial(MyMaterial,MyGeometryName,2 );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   create mesh



% create mesh
%  create mesh
mesh2=model.mesh.create('mesh2', 'geom1');
meshSize2=mesh2.feature('size');
meshSize2.set('hauto', '4');
mesh2.feature.create('ftri1', 'FreeTri');

%% mesh everywere else
   model.mesh('mesh2').feature('size').set('custom', 'on');
model.mesh('mesh2').feature('size').set('hmax', ['lambda/' num2str(nSubstrate) '/4']);
model.mesh('mesh2').feature('size').set('hmin', ['lambda/' num2str(nSubstrate) '/8']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

modeSearchNumber=1;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set aproximate refractive index
nAprox=nSiO2(lambda);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  physics for Hp
model.physics.create('emw', 'ElectromagneticWaves', 'geom1');

model.study.create('std1');
model.study('std1').feature.create('mode', 'ModeAnalysis');

model.physics('emw').feature('wee1').set('DisplacementFieldModel', 1, 'RefractiveIndex');
model.physics('emw').feature('wee1').set('n_mat', 1, 'from_mat');

% set all boundary scatering boundary conditions
model.physics('emw').feature.create('sctr1', 'Scattering', 1);
model.physics('emw').feature('sctr1').selection.all;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.study('std1').feature('mode').set('geomselection', 'geom1');
model.study('std1').feature('mode').set('physselection', 'emw');
model.study('std1').feature('mode').set('shift', num2str(nAprox));  %  aproximate mode refractive index
model.study('std1').feature('mode').set('modeFreq', ['c_const/' num2str(lambda) '[um]']);
model.study('std1').feature('mode').set('neigs', num2str(modeSearchNumber));    % number of modes to search



model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'mode');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('e1', 'Eigenvalue');
model.sol('sol1').feature('e1').set('shift',  num2str(nAprox));   %  aproximate mode refractive index
model.sol('sol1').feature('e1').set('neigs', modeSearchNumber);                            %  number of mode to search
model.sol('sol1').feature('e1').set('transform', 'effective_mode_index');
model.sol('sol1').feature('e1').set('control', 'mode');
model.sol('sol1').feature('e1').feature('aDef').set('complexfun', true);
model.sol('sol1').attach('std1');

%ModelUtil.showProgress(true)
model.sol('sol1').runAll;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% find which solutions are realistic
n_eff_all = mphglobal(model,'emw.neff','Complexout','on');
jMode=0;
for j=1:length(n_eff_all)
    if real(n_eff_all(j))>0       
        jMode=jMode+1;
        mode(jMode).neffHp=n_eff_all(j);
        mode(jMode).solutionN=j;
    end
end
disp([' mode total number:  ' num2str(jMode)])

%





for j=1:jMode


    
    
  
    mode(j).neff=mode(j).neffHp;
    mode(j).LossHp=kdb(-imag(mode(j).neffHp),lambda)/1E4; %dB/um

    mode(j).PropagationDistance=-lambda/2/pi/imag(mode(j).neffHp); % in um 1/e propagation distance
  


end


disp('calculations OK')


n_eff_all = mphglobal(model,'emw.neff','Complexout','on');
jMode=0;
for j=1:length(n_eff_all)
    %if real(n_eff_all(j))>nAlGaAs(xClad,lambda)      
        jMode=jMode+1;
        mode(jMode).neff=real(n_eff_all(j));
        mode(jMode).loss=kdb(imag(n_eff_all(j)),lambda)/10;  %loss in dB/mm
       
        mode(jMode).n=n_eff_all(j);
        mode(jMode).solutionN=j;
   % end
end




%%%%%% find whether the mode TE or TM
jTM=0;jTE=0;jTEM=0;


for j=1:jMode
    
mode(j).TEintens=mphint(model,'abs(emw.Ex)^2','Solnum',num2str(mode(j).solutionN));  % integration over all domains
mode(j).TMintens=mphint(model,'abs(emw.Ey)^2','Solnum',num2str(mode(j).solutionN));
if mode(j).TEintens>3*mode(j).TMintens 
    mode(j).Polar='TE';
    jTE=jTE+1;TEmod(jTE,1)=j;TEmod(jTE,2)=mode(j).neff;
    
elseif mode(j).TMintens>3*mode(j).TEintens
    mode(j).Polar='TM';
    jTM=jTM+1;TMmod(jTM,1)=j;TMmod(jTM,2)=mode(j).neff;
else
    mode(j).Polar='TEM';
    jTEM=jTEM+1;TEMmod(jTEM,1)=j;TEMmod(jTEM,2)=mode(j).neff;
end

end

% put each mode in order, mode with highest refractive index first
mes='';
if jTE>0 TEmod=sortrows(TEmod,-2); mes=['number of TE modes ' num2str(size(TEmod,1))];end
if jTE>0 TMmod=sortrows(TMmod,-2);mes=[mes ' of TM modes ' num2str(size(TMmod,1))]; end
 if jTEM>0 TEMmod=sortrows(TEMmod,-2);mes=[mes ' of TEM modes ' num2str(size(TEMmod,1))]; end

disp(mes )



jplot=3; 
if jTE>0
for j=1:size(TEmod,1)
    jplot=jplot+1;
dataName=['pg' num2str(jplot)];
model.result.create(dataName, 'PlotGroup2D');
model.result(dataName).feature.create('surf1', 'Surface');
model.result(dataName).feature('surf1').set('descr', ['Electric field TE' num2str(j-1) ' mode, x component loss' num2str(mode(TEmod(j,1)).loss) ' dB/mm']);
model.result(dataName).set('title', [' TE' num2str(j-1) ' mode, abs(Ex)  n= ' num2str(mode(TEmod(j,1)).neff) 'loss: ' num2str(mode(TEmod(j,1)).loss) ' dB/mm']);
model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ex)');
model.result(dataName).name(['TE' num2str(j-1)]);
model.result(dataName).set('solnum', num2str(TEmod(j,1)));
%model.result('dataName').set('data', 'dset1');
model.result(dataName).run;


end 
end

if jTM>0
for j=1:size(TMmod,1)
    jplot=jplot+1;
dataName=['pg' num2str(jplot)];
model.result.create(dataName, 'PlotGroup2D');
model.result(dataName).feature.create('surf1', 'Surface');
model.result(dataName).feature('surf1').set('descr', ['Electric field TE' num2str(j-1) ' mode, x component']);
MyTitle= [' TM' num2str(j-1) ' mode, abs(Ey)  n= ' num2str(mode(TMmod(j,1)).neff) 'loss: ' num2str(mode(TMmod(j,1)).loss) ' dB/mm'];
model.result(dataName).set('title', MyTitle);
model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ey)');
model.result(dataName).name(['TM' num2str(j-1)]);
model.result(dataName).set('solnum', num2str(TMmod(j,1)));
%model.result(dataName).set('data', 'dset1');

myMaxV= mphinterp(model,{'abs(emw.Ey)'},'coord',[0; -0.1],'unit',{'V/m'},'solnum', num2str(TMmod(j,1)));   % data at center
model.result(dataName).feature('surf1').set('rangecoloractive', 'on');
model.result(dataName).feature('surf1').set('rangecolormin', '0');
model.result(dataName).feature('surf1').set('rangecolormax', num2str(myMaxV*1.2));


end 
end

disp('OK')
return