%Physical Constants
mu0      = 4*pi*1e-7;  % Free space permiability H/m
epsilon0 = 8.85418782; % Free space permittivity F/m
eta0    =  sqrt(mu0/epsilon0); %Free space impedance

frequency = 433e6;  
C = 299792458; % Velocity of light in free space m/s
beta = 2*pi*frequency / C; % Wave number

%dimensions: A square cavity is considered 
xa = -1:0.01:1;
ya = xa;
za = xa;
x_src = 0;
y_src = 0;
z_src = 0;
[xx,yy,zz] = meshgrid(xa,ya,za);
%Find the distance from the source
Rsq        = (xx-x_src).^2+(yy-y_src).^2+(zz-z_src).^2;
R          = sqrt(Rsq);

% At the point of singularity i.e, at R=0, set it to a small value
ind        = find(R(:)==0);
[i,j,k]    = ind2sub(size(xx),ind);
Rsq(i,j,k) = 1e-6;
R(i,j,k)   = 1e-3;

%Evaluate cosine and sine of angles in spherical coordinates
costh      = zz./sqrt(R); %cos(theta)
sinth      = sqrt((xx-x_src).^2+(yy-y_src).^2)./R;
phi        = atan2(yy,xx);
cosphi     = cos(phi);
sinphi     = sin(phi);



%Safety check to limit cos(phi) and sin(phi) with in [-1,1]
cosphi     = max(-1,cosphi);
cosphi     = min(cosphi,1);
sinphi     = max(-1,sinphi);
sinphi     = min(1,sinphi);
costh      = max(-1,costh);
costh      = min(1,costh);
sinth      = min(1,sinth);

%Dipole moment P = I*dl
P = 0.01; % in A.m

%Field calculations
f_r       = (eta0*P/2/pi)* costh./R.^2 .*(1-1j/beta./R).*exp(-1j*beta*R);
f_th      = (1j*eta0*beta*P/4/pi)* sinth./R .*(1-1j/beta./R-1/beta^2./Rsq).*exp(-1j*beta*R); 

%Use Spherical to cartesian transformation 
fx   = f_r.*sinth.*cosphi + f_th.*costh.*cosphi;
fy   = f_r.*sinth.*sinphi + f_th.*costh.*sinphi;
fz   = f_r.*costh - f_th.*sinth;
figure(1)
mesh(xa,ya,abs(fz(:,:,k)));

%Load the comsol exported file
t1= load('Ez_pml.txt');
x1= t1(:,1);
y1= t1(:,2);
fz_comsol= t1(:,4)+1j*t1(:,5);

F = scatteredInterpolant(x1,y1,abs(fz_comsol),'linear','nearest');
x = linspace(min(x1),max(x1),100);
y = linspace(min(y1),max(y1),100);
[X,Y] = meshgrid(x,y);
Z = F(X,Y);
figure(2)
mesh(X,Y,Z)