%%  SET UP PARAMETER

mesh.name ='Mesh_Kalscheuer2008_1ano_v2';% 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     = [2  4];      %  Domain ID of air
P.sig_air   = 1e-8;        % Conductivity of air S/m
P.sigp      = 1./[10000];  % Background (residual domain) conductivity
P.sig_anoID = [5];         % Domain number of the abnormal body
P.sigax    = 1./[10000];   % Conductivity of the abnormal body 
P.sigay    = 1./[10000] ;
P.sigaz    = 1./[10000];

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   = [aS];     
P.aD   = [aD];
P.aL   =  [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    = [];        %: anomalous domain ID for comsol mesh,
P.mura        = [];
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      = 1;
P.eps_anoID  = [5];        % Domain ID of the abnormal body    
P.epsax      = [epsrx];   % dielectric permittivity of the abnormal body
P.epsay      = [epsry];   
P.epsaz      = [epsrz];
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);