%   Example - e08 Radiation force on the solid plate
%   as a function of frequency

clear all
setpath

%   Define scatterer geometry - plate
S		= .1;               %   Size
[P,t]	= g_plate(S,S,20,20,0);
P(2,:)  = P(1,:);
P(1,:)  = 0;
geom    = rwgm(P,t,4);      % basis functions

%   Define constants and the incident field
const.epsilon       = 8.854e-012;
const.mu            = 1.25664e-006;
const.c             = 1/sqrt(const.epsilon*const.mu);
const.eta           = sqrt(const.mu/const.epsilon);

const.dir           = [1; 0; 0];    %Direction
const.pol           = [0; 0; 1];    %Polarization

%   Define observation points in the near field
Points  =   geom.P;    

%   Define unperturbed force
Force0  =   S*S/(2*const.eta)/const.c*2; %double reflection

%   Define frequency and solve the problem
frequency       = 3e9;
I                   = solver(geom,const,frequency);
[Poynting,Ev,Hv]    = field(const,geom,I,frequency,Points,1);
Ep                  = (Ev(:,geom.t(1,:))+Ev(:,geom.t(2,:))+Ev(:,geom.t(3,:)))/3;
Hp                  = (Hv(:,geom.t(1,:))+Hv(:,geom.t(2,:))+Hv(:,geom.t(3,:)))/3;

for m =1:geom.TrianglesTotal
    Ep(:,m) = (Ev(:,geom.t(1,m)) + Ev(:,geom.t(2,m)) + Ev(:,geom.t(3,m)))/3;
    Hp(:,m) = (Hv(:,geom.t(1,m)) + Hv(:,geom.t(2,m)) + Hv(:,geom.t(3,m)))/3;
end
    
[ForceB, dummy]     = force (geom,const,frequency,I,Ep,Hp,0);
[ForceT, dummy]     = force (geom,const,frequency,I,Ep,Hp,1);
for i = 1:3
        FTB(i)            = sum(ForceB(i,:).*geom.AreaF)/Force0
        FTT(i)            = sum(ForceT(i,:).*geom.AreaF)/Force0
end