%EHKFUN Returns the fields for the given plane wave in a homogeneous
%material; 
%     eps = relative permittivity
%     mu=relative permeability
%     omg= angular frequency
%     lx,ly= waveguide size
%     nx, ny = mode numbers
%     aE,aH = n

% c0const=3e8; mu0const=4e-7*pi; 

lam0=1; 
eps=1; mu=1; 
lx=1; 
ly=1; 

%setup plotting field
xVec=linspace(0,lx,15); 
yVec=linspace(0,ly,15); 
[xMesh,yMesh]=meshgrid(xVec,yVec); 

xVec=linspace(0,lx,50); 
yVec=linspace(0,ly,50); 
[xMeshZ,yMeshZ]=meshgrid(xVec,yVec); 

nyLst=(1:40); nxLst=0+0*nyLst; 
aHLst=0.1*ly*((-1).^nyLst-1); aELst=0*aHLst; 
% aELst=0.1*ly*((-1).^nyLst-1)./nyLst; aHLst=0*aELst; 

% nxLst=0; nyLst=1; 
% aHLst=1; aELst=0; 

for z0=(0:0.02:3)
ExMesh=0*xMesh; 
EyMesh=0*xMesh; 
EzMesh=0*xMeshZ; 


for imd=1:length(nxLst)
%modal parameters
nx=nxLst(imd); ny=nyLst(imd); 
aE=aELst(imd); aH=aHLst(imd); 


qx=pi*nx/lx; qy=pi*ny/ly; 
kz=sqrt(eps*mu*(2*pi/lam0)^2-qx^2-qy^2); %calculate kz
% enforce sign of kz
kz(imag(kz)<0)=-kz(imag(kz)<0);

% calculate the fields
Ez=aE*sin(qx*xMeshZ).*sin(qy*yMeshZ); 
Ex=(kz*lam0/2/pi/eps*aE*qx-aH*qy)*cos(qx*xMesh).*sin(qy*yMesh)/(qx^2+qy^2);
Ey=(kz*lam0/2/pi/eps*aE*qy+aH*qx)*sin(qx*xMesh).*cos(qy*yMesh)/(qx^2+qy^2);

ExMesh=ExMesh+Ex*exp(1i*kz*z0); EyMesh=EyMesh+Ey*exp(1i*kz*z0); EzMesh=EzMesh+Ez*exp(1i*kz*z0); 
end 

ExMesh=real(ExMesh); EyMesh=real(EyMesh); EzMesh=real(EzMesh);
zmax=0.25; EzMesh(1)=zmax; EzMesh(2)=-zmax; 
EzMesh(EzMesh>zmax)=zmax; EzMesh(EzMesh<-zmax)=-zmax; 

figure(1)
clf
surf(xMeshZ,yMeshZ,EzMesh,'EdgeColor','none')
set(gca,'FontSize',18)
view(2)
hold on 
quiver3(xMesh,yMesh,max(abs(EzMesh(:)))+0*ExMesh,ExMesh,EyMesh,0*ExMesh, 0,...
    'color',0*[1 1 1])
hold off
zlim([-zmax zmax])
colorbar


xlim([0 lx]); ylim([0 ly]); daspect([1 1 1])

pause(0.2);
end