flclear fem
%clc 
close all
clear all
tStart=tic; 

%Geometry

%Create structure array for geometry dimensions
gd = struct('length_tot',10e-3,...    %sim raum
    'w_tot',30e-3,...
    'h_tot',20e-3,...
    'h_sbstr',0.64e-3,...                   %dielectric layer  thickness (substrate)
    'w_mstrip',5e-3);                       %microstrip width thickness


%air
g1=block3(gd.length_tot,gd.w_tot,gd.h_tot-gd.h_sbstr,'base','center','pos',[0,0, (gd.h_tot)/2]);
%dielectric
g2=block3(gd.length_tot,gd.w_tot,gd.h_sbstr,'base','center','pos',[0 0 0]);
%microstrip
g3 = face3(rect2(gd.length_tot,gd.w_mstrip,'base','center','pos',[0 gd.h_sbstr/2]));
g4 =  move(g3,[0, 0 ,gd.h_sbstr/2]);


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

scrsz = get(0,'ScreenSize');

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')

% Constants
fem.const = {'k','0'}; %6.3

%Create structure array for output
out_vars = struct('lambda',[],'omega',[],'damp',[],'Qfact',[],'nu',[]);

freq0 = 9.84e7; %Eigref Linearization point


%% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'hmaxedg',[9,1e-3,10,1e-3,12,1e-3,23,1e-3], ...
                  'hmaxfac',[8,1e-3], ...
                  'point',[], ...
                  'edge',[], ...
                  'face',[1,4], ...
                  'subdomain',[]);
              
% Copy boundary mesh
fem.mesh=meshcopy(fem,'source',1,'target',12,'mcase',0);
fem.mesh=meshcopy(fem,'source',4,'target',13,'mcase',0);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'hmaxedg',[9,1e-3,10,1e-3,12,1e-3,23,1e-3], ...
                  'hmaxfac',[8,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)



k_x=linspace(0,2*pi/(2*gd.length_tot),100); %k-vector for floquet periodicity

for j_loop=1:length(k_x);
    
    % 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 = {'periodic','E0','cont'};
    bnd.kper = {{k_x(j_loop);0;0},{0;0;0},{0;0;0}};
    bnd.pertype = {'floque','sym','sym'};
    bnd.ind = [1,2,2,1,2,3,2,2,3,2,2,1,1];
    appl.bnd = bnd;

    clear equ
    equ.epsilonr = {10,1};
    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);
    fem.xmesh = meshextend(fem);

    %DOFs limitation
    if (xmeshinfo(fem,'out','ndofs') > 100000)
        disp(['Number of DOFs : ', num2str(xmeshinfo(fem,'out','ndofs'))]);
        disp('Exittinng, number of DOFs is too much...')
        return
    end 

%% Solve problem
    n_fine=4; %eigenfrequency fine search parameter;
    
    for i_loop=1:1:n_fine;
        fem.sol = femeig(fem,'neigs',1,'shift',freq0(i_loop),'eigref',sprintf('-i*2*pi*%g',freq0(i_loop)));
        out_vars.lambda=fem.sol.lambda;
        out_vars.omega=imag(-out_vars.lambda);           %Angular frequency
        out_vars.damp=real(out_vars.lambda);             %Damping in time
        out_vars.Qfact=0.5*out_vars.omega./abs(out_vars.damp);    %Quality factor
        out_vars.nu=out_vars.omega/(2*pi);               %Frequency
        freq0(i_loop+1) = out_vars.nu;
    end
  
    disp(['Loop: ',num2str(j_loop), '	k: ',num2str(k_x(j_loop)),'	Frq:',num2str(freq0, '%5.2e     ') ]);
    if freq0(n_fine) < 0
        freq0=9.84e7;
    else
        freq0=freq0(n_fine); %for next eigenfrequency search
    end   

    loop_out_vars.solutions(j_loop)=out_vars;
end


%%

figure('Name','Solution Plot 1','Position',[2*scrsz(3)/3 scrsz(4)*2/3-30 ...
scrsz(3)/3 scrsz(4)/3],'MenuBar','none');
%extacting the results from the solutions for plotting
for i=1:length(k_x)
    freq_pl(i)=loop_out_vars.solutions(i).nu;
    qfact_pl(i)=loop_out_vars.solutions(i).Qfact;
    hold on
end

subplot(211)
plot(k_x(1,1:10),freq_pl(1,1:10));
xlabel('k[rad/m]');
ylabel('Frequency [Hz]');

subplot(212)
plot(k_x(1,1:10),qfact_pl(1,1:10));
xlabel('k[rad/m]');
ylabel('Quality Factor');


%% Plot solution
figure('Name','Solution Plot 2','Position',[0*scrsz(3)/3, ...
    scrsz(4)*1/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')  
tElapsed=toc(tStart);
disp('----------DONE-------------');
disp(['Elepsed Time :',num2str(tElapsed/60), ' min ' ]);