% This Document is created by Hakan Cakmak, it is for test purpose and not
% realy working if you find a better result please contact me
% cakmak.hakan@yahoo.com

% flclear fem
% 
%clc 
close all
clear all

length_tot=20e-3;    %sim raum
w_tot=20e-3;         %sim raum
h_tot=10e-3;         %sim raum
H=0.64e-3;          %dielectric layer  thickness (substrate)
w=5e-3;             %microstrip width thickness
%Geometry

pml_thicknes=0.003;

%sim raum
g1=block3(length_tot,w_tot,h_tot,'base','center','pos',[0,0, (h_tot+H)/2]);

%dielectric
g2=block3(length_tot,w_tot,H,'base','center','pos',[0 0 0]);

%microstrip
g3 = face3(rect2(length_tot,w,'base','center','pos',[0 H/2]));
g4 =  move(g3,[0, 0 ,H/2]);

%PML
g5=block3(length_tot,w_tot,pml_thicknes,'base','center','pos',[0,0, (pml_thicknes+H)/2+h_tot]);





clear f 
f.objs={g1,g2,g4,g5};
f.name={'g1','g2','g4','g5'};
f.tags={'g1','g2','g4','g5'};
fem.draw=struct('f',f);
fem.geom=geomcsg(fem);






scrsz = get(0,'ScreenSize');%1,1,1600,1200%
figure('Name','Simulation Plot','Position',[1 scrsz(4)*2/3-30 ...
    scrsz(3)/3 scrsz(4)/3],'MenuBar','none');
hold on

geomplot(g1,'facemode','on','facelabels','on','transparency',0.4,'camlight','on','renderer','opengl','axisvisible','off') 
geomplot(g2,'facemode','on','facelabels','on','transparency',0.4,'camlight','on','renderer','opengl','axisvisible','off') 
geomplot(g4,'facemode','on','facelabels','on','transparency',0.4,'camlight','on','renderer','opengl','axisvisible','off')
geomplot(g5,'facemode','on','facelabels','on','transparency',0.4,'camlight','on','renderer','opengl','axisvisible','off') 

pause(1)
% Constants
fem.const = {'k','6.3'};

%return
% % Initialize mesh
% fem.mesh=meshinit(fem, ...
%                   'hauto',7, ...
%                   'hmaxedg',[12,1e-3,13,1e-3,15,1e-3,30,1e-3], ...
%                   'hmaxfac',[11,1e-3], ...
%                   'point',[], ...
%                   'edge',[12,13,15,30], ...
%                   'face',[], ...
%                   'subdomain',[]);
% 
% % Initialize mesh
% fem.mesh=meshinit(fem, ...
%                   'hauto',7, ...
%                   'hmaxedg',[12,1e-3,13,1e-3,15,1e-3,30,1e-3], ...
%                   'hmaxfac',[11,1e-3], ...
%                   'point',[], ...
%                   'edge',[], ...
%                   'face',[11], ...
%                   'subdomain',[], ...
%                   'meshstart',fem.mesh);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'hmaxedg',[12,1e-3,13,1e-3,15,1e-3,30,1e-3], ...
                  'hmaxfac',[11,1e-3], ...
                  'point',[], ...
                  'edge',[], ...
                  'face',[1,4,7], ...
                  'subdomain',[]);
              
% Copy boundary mesh
fem.mesh=meshcopy(fem,'source',1,'target',16,'mcase',0);
fem.mesh=meshcopy(fem,'source',4,'target',17,'mcase',0);
fem.mesh=meshcopy(fem,'source',7,'target',18,'mcase',0);
              


% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'hmaxedg',[12,1e-3,13,1e-3,15,1e-3,30,1e-3], ...
                  'hmaxfac',[11,1e-3], ...
                  'point','auto', ...
                  'edge','auto', ...
                  'face','auto', ...
                  'subdomain','auto', ...
                  'meshstart',fem.mesh);
              
%mesh plot
figure('Name','Mesh Plot','Position',[scrsz(3)/3 scrsz(4)*2/3-30 ...
    scrsz(3)/3 scrsz(4)/3],'MenuBar','none');
meshplot(fem,'transparency',0.4), axis equal

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'ElectromagneticWaves';
appl.module = 'RF';
appl.gporder = 4;
appl.cporder = 2;
appl.border = 'on';
appl.assignsuffix = '_rfw';
clear prop
prop.analysis='eigen';
prop.divcond='on';
appl.prop = prop;

clear bnd
bnd.type = {'E0','cont','periodic','SC','H0'};
bnd.kper = {{0;0;0},{0;0;0},{'k';0;0},{0;0;0},{0; ...
  0;0}};
bnd.pertype = {'sym','sym','floque','sym','sym'};
bnd.ind = [3,5,1,3,5,2,3,5,2,4,1,2,5,5,5,3,3,3];
appl.bnd = bnd;

clear equ
equ.coordOn = {{0;0;0},{0;0;0},{0;0;1}};
equ.epsilonr = {10,1,1};
equ.Stype = {'none','none','coord'};
equ.ind = [1,2,3];
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);

return
fem.xmesh = meshextend(fem);
fem.sol = femeig(fem,'neigs',1,'shift',1);

% Neigs:Number of eigenvalues sought
% Out:Output variables
% Shift:Eigenvalue search location
% Eigname:Name of eigenvalue variable
% Eigref:Linearization point for the eigenvalue
% Etol:Eigenvalue tolerance
% In:Input matrices
% Krylovdim:Dimension of Krylov space

% figure('Name','Solution Plot','Position',[2*scrsz(3)/3 scrsz(4)*2/3-30 ...
%     scrsz(3)/3 scrsz(4)/3],'MenuBar','none');
% postplot(fem, ...
%          'slicedata',{'Ez','cont','internal','unit','V/m'}, ...
%          'slicexspacing',1, ...
%          'sliceyspacing',1, ...
%          'slicezspacing',[1e-3,0.1e-3], ...
%          'slicemap','Rainbow','transparency',0.2,'renderer','opengl')     
%