close all;
clear all;
clc

%addpath('./Functions');
addpath('../MatlabFemBuild');

eps0=8.854187817e-12;
mu0=1.256637061e-6;


SymOrder=10;
Resolution=20;

%Geo-Setup
rad=0.667e-6;   %radius of particle
lambda=1.55e-6; %excite wavelength
epsR=10;        %epsilon_r of particle
glob_epsR=1;    %surrounding eps_r
optDist=2e-9;   %boundary meshing distance

mesh=lambda/Resolution/sqrt(glob_epsR);
submesh=lambda/Resolution/sqrt(epsR);
Np=round((pi/SymOrder*rad/optDist));

sizeSimBox=3e-6;    %size of simulation area (sim box)
sizePML=lambda/2;   %size of pml

sizeSB=sizeSimBox+2*sizePML;

%init particle
alpha=linspace(0,pi/SymOrder,Np);
x=-cos(alpha)*rad;
y=sin(alpha)*rad;

    %areaSize=polyarea(x,y);
    
%Sweep-Setup
Nlbd=51;
lambdaBW=0.02e-6;  %simulation bandwidth
fitBW=0.002e-6;  %simulation bandwidth
lambdaSw=linspace(lambda-lambdaBW/2,lambda+lambdaBW/2,Nlbd);
Wl=zeros(2,Nlbd);
    
   
%do Simulation
disp(['fem-Simulation']);
args={lambda,glob_epsR,epsR,sizeSB,sizePML,mesh,submesh,SymOrder};
femHandle=@femCreate_new_v1;
femSweepHandle=@femFreqSweep_new_v1;

fem = femHandle(x,y,args{:},1);

    
postplot(fem, ...
         'tridata',{'normH_rfwh','cont','internal','unit','A/m'}, ...
         'trimap','jet(1024)', ...
         'title','Surface: Magnetic field, norm [A/m]', ...
         'sdl',{[3,5]});
xlim([-sizeSimBox/2 sizeSimBox/2]);
ylim([0 sizeSimBox/2]);
axis equal
pause(0.1);

[a,b,c,e]=FemLorentzFit(x,y,lambda,fitBW,femHandle,femSweepHandle,args);

ParticleSD = GetDrawSubDomains(fem,'Particle');
W = femSweepHandle(fem,lambdaSw,ParticleSD);

figure
plot(lambdaSw,W/max(W));
hold on;
grid on;
plot(lambdaSw,a./((lambdaSw-b).^2+c)*(c/a),'r');
xlim([min(lambdaSw) max(lambdaSw)]);