%  Cloak5 - PML's all around - no target
%======================================================================================================

nEl = 20;
bSolve = 0;

flclear fem
tic

%  Material Parameters
rho0 = 1;
kappa0 = 1;

%  Frequency Parameters
lambda = 1;
k = 2*pi/lambda;

%  PML Parameters
hp   = 1.5*lambda;
pml  = 2;
pmlx = '1/(1+pml*j*min((x/hp)^2,((x-hx)/hp)^2))';
pmly = '1/(1+pml*j*min((y/hp)^2,((y-hy)/hp)^2))';

c0  = 1;  c1 = 10;
R1  = 2/3;  R2 = 2*R1;
r   = 'sqrt((x-x0)^2+(y-y0)^2)';
cs  = '(x-x0)/r';
sn  = '(y-y0)/r';
pr  = '(r-R1)/r';
pt  = 'r/(r-R1)';  
p11 = 'pr*cs^2+pt*sn^2';
p12 = 'cs*sn*(pr-pt)';
p22 = 'pr*sn^2+pt*cs^2';
kappa = '((R2-R1)/R2)^2*r/(r-R1)*(1-0.0*j)';
%kappa = kappa0;

%  Plate Parameters
hx  = 8;                                                %  Width
hy  = 16/3;                                                %  Height

%  Targets                                               
tgtx = hx/2;
tgty = hy/2;                                              
tgtr = R1+1.e-6;                                              

clkx = tgtx;
clky = tgty;
clkr = R2;

%  Target Region
gPlate = rect2( hx, hy, 'base','corner', 'pos',[0,0], 'rot','0' );

%  Add Targets + Cloak
clkThick = clkr-tgtr;
thickInc = clkThick/5;
gTgt = circ2( tgtr, 'base','center', 'pos',[tgtx,tgty], 'rot','0' );
ii = 1;
rad = tgtr + ii*thickInc;
gLay1 = circ2( rad, 'base','center', 'pos',[clkx,clky], 'rot','0' );
ii = 2;
rad = tgtr + ii*thickInc;
gLay2 = circ2( rad, 'base','center', 'pos',[clkx,clky], 'rot','0' );
ii = 3;
rad = tgtr + ii*thickInc;
gLay3 = circ2( rad, 'base','center', 'pos',[clkx,clky], 'rot','0' );
ii = 4;
rad = tgtr + ii*thickInc;
gLay4 = circ2( rad, 'base','center', 'pos',[clkx,clky], 'rot','0' );
ii = 5;
rad = tgtr + ii*thickInc;
gClk = circ2( rad, 'base','center', 'pos',[clkx,clky], 'rot','0' );
gPlate = geomcomp( {gPlate,gClk},'ns',{'gPlate','gClk'}, 'sf','gPlate-gClk', 'edge','none' );
gClk = geomcomp( {gClk,gLay4},'ns',{'gClk','gLay4'}, 'sf','gClk-gLay4', 'edge','none' );
gLay4 = geomcomp( {gLay4,gLay3},'ns',{'gLay4','gLay3'}, 'sf','gLay4-gLay3', 'edge','none' );
gLay3 = geomcomp( {gLay3,gLay2},'ns',{'gLay3','gLay2'}, 'sf','gLay3-gLay2', 'edge','none' );
gLay2 = geomcomp( {gLay2,gLay1},'ns',{'gLay2','gLay1'}, 'sf','gLay2-gLay1', 'edge','none' );
gLay1 = geomcomp( {gLay1,gTgt},'ns',{'gLay1','gTgt'}, 'sf','gLay1-gTgt', 'edge','none' );

%%  Edge PMLs
gpl = rect2( hp, hy, 'base','corner', 'pos',[-hp, 0], 'rot','0' );     %  Left
gpr = rect2( hp, hy, 'base','corner', 'pos',[ hx, 0], 'rot','0' );     %  Right
gpb = rect2( hx, hp, 'base','corner', 'pos',[0, -hp], 'rot','0' );     %  Top
gpt = rect2( hx, hp, 'base','corner', 'pos',[0,  hy], 'rot','0' );     %  Bottom
%
%%  Corner PMLs
gpbl = rect2( hp, hp, 'base','corner', 'pos',[-hp, -hp], 'rot','0' );  %  Bottom Left
gptl = rect2( hp, hp, 'base','corner', 'pos',[-hp,  hy], 'rot','0' );  %  Top Left 
gpbr = rect2( hp, hp, 'base','corner', 'pos',[ hx, -hp], 'rot','0' );  %  Bottom Right
gptr = rect2( hp, hp, 'base','corner', 'pos',[ hx,  hy], 'rot','0' );  %  Top Right

% Analyzed geometry
clear s

%  Subdomain PML
s.objs = {gPlate, gClk, gLay4, gLay3, gLay2, gLay1, gpl,gpr,gpb,gpt, gpbl,gpbr,gptl,gptr};
SubPml = [4,5, 5, 5, 5, 5, 2,2,3,3, 1,1,1,1]; 
%
fem.draw=struct('s',s);
[g,st,ft] = geomcsg(fem);
[SubInd,s0] = find(st);                                                 %  SubInd gives subdomain indices
fem.geom=g;

SubAll = ones( 1, size(st,1) );                                         %  Used to assign values to all 
BndAll = ones( 1, size(ft,1) );                                         %  subdomains, boundaries, and points

%  Equation:  Wave Equation - Weak Form
pdeStart = '1/dx_pml/dy_pml*(';                                                                        %  Jacobian
pde1 = '-kappa*wx_test*dx_pml*(p11*wx*dx_pml+p12*wy*dy_pml)';
pde2 = '-kappa*wy_test*dy_pml*(p12*wx*dx_pml+p22*wy*dy_pml)';                               
pde3 = '+k^2*w*w_test';                                                 %  Time Dependent Terms
pdeEnd = ')';
pde = strcat(pdeStart,pde1,pde2,pde3,pdeEnd);

%  Create Mesh
%hmax = lambda / nEl;
%fem.mesh=meshinit( fem, 'hmax',[hmax] );

%  Weak PDE Application Mode 
clear appl
appl.mode.class = 'FlPDEW';
appl.dim = {'w','w_t'};
appl.border = 'on';
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_w';

%  Differential Equation - Weak Form - The same PDE applies to all subdomains - dx_pml and dy_pml are subdomain dependent
clear equ
equ.weak = pde;
equ.dweak = 0;
equ.ind = SubAll;
appl.equ = equ;

%  Subdomain settings - PML Variables
clear equ
equ.ind( SubInd ) = SubPml;
equ.dim = {'w'};
equ.expr = { 'dx_pml',{pmlx,pmlx,1,1,1}, 'dy_pml',{pmly,1,pmly,1,1}, 'kappa',{c0,c0,c0,c0,kappa}, ...
             'p11',{1,1,1,1,p11}, 'p12',{0,0,0,0,p12}, 'p22',{1,1,1,1,p22} };
fem.equ = equ;

%  Boundary Settings - Simply Support Edges Along the Boundary
clear bnd
bnd.constr = {'0','w','w-1'};
bnd.ind = BndAll;
bnd.ind( CheckBound(g) ) = 2;
%bnd.ind( [27,28,30,31] ) = 2;
bnd.ind( [35,36,42,43] ) = 2;
bnd.ind( [8,10,12] ) = 3;
appl.bnd = bnd;

%  Give Comsol Relevant FEM Constants
fem.const = {'R1',R1, 'R2',R2, 'k',k, 'pml',pml, 'hp',hp, 'hx',hx, 'hy',hy, 'x0',tgtx, 'y0',tgty };
fem.globalexpr = {'r',r, 'cs',cs, 'sn',sn, 'pr',pr, 'pt',pt };

%  Final Setup - Units, Comsol Stuff
clear units;
units.basesystem = 'SI';
fem.units = units;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
fem=multiphysics(fem);