%%  SET PARAMETER
addpath(genpath('D:\em3d_MT'));

mesh.name  = 'mesh_layeredmodel'; % Contains grid information for the model, output by GeoMesh
[mesh.node2coord,mesh.elem2node,mesh.entity] = load_comsol(mesh.name);
% load_mesh;

%% Set the conductivity
P.airID     = [4];              % Domain ID of air
P.sig_air   = 1e-8;             % Conductivity of air S/m
P.sigp      = 1./[2000];        % Conductivity of background (residual domain) 
P.sig_anoID = [2 1];            % Domain ID of the abnormal body
P.sigax     = 1./[20000 20000]; % Conductivity of the abnormal body  
P.sigay     = 1./[20000 20000];
P.sigaz     = 1./[20000 20000];

P.sig = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.sigax,P.sigay,P.sigaz,P.sig_air,P.sigp) ;

% Set the Euler rotation Angle of the abnormal body 
P.aS   = [0 0];     
P.aD   = [45 45];
P.aL   = [0 0];
P.AniAngle = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.aS,P.aD,P.aL,0,0) ;
P.sig  = get_sig_ani(P.sig(:,1),P.sig(:,2),P.sig(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));

%% Set magnetic permitivity
P.mu0 = 4*pi*1e-7 ;  % V·s/ (A·m)
P.mu_anoID   = [];  % Domain ID of the abnormal body
P.mura       = [];  % magnetic permitivity of the abnormal body
P.mur = Set_material_3Ani(mesh.entity,P.mu_anoID,P.airID,P.mura,P.mura,P.mura,1,1) ;
P.mur = P.mur(:,1);
P.mu  = P.mur.*P.mu0 ;

%% Set dielectric permittivity
P.epsr_air   = 1;    
P.epsrp      = 8;       % dielectric permittivity of background (residual domain) 
P.eps_anoID  = [2 1];   % Domain ID of the abnormal body
P.epsax      = [8 8];   % dielectric permittivity of the abnormal body
P.epsay      = [5 5];
P.epsaz      = [6 6];
P.epsr = Set_material_3Ani(mesh.entity,P.eps_anoID,P.airID,P.epsax,P.epsay,P.epsaz,P.epsr_air,P.epsrp) ;
P.epsr = get_sig_ani(P.epsr(:,1),P.epsr(:,2),P.epsr(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));
P.eps0 = 8.854187817620389*10^-12 ;
P.eps  = P.epsr *P.eps0 ;

%% Set frequency
freq  = logspace(4,6,11);
omega = 2*pi*freq;
nfreq = size(freq,2);