Clear["Global`*"]
(* visualize rectangular waveguide solutions:
propagation direction is z, we consider TE modes. 0<x<a, 0<y<b, b<a *)
(* © Marko Horbatsch, York University, Toronto, Canada *)

c=3 10^8;

cm=10^-2; 
(* SI we could define the inch instead of cm as our length unit for the guide parameters a and b *)

a=1 cm;
b=2/3 cm;

pi=N[Pi];
tpi=2pi

ωmn[m_,n_]:=pi c Sqrt[(m/a)^2+(n/b)^2]

{ωmn[1,0],ωmn[0,1]}
(* the interval for single-mode operation of the guide *)

f10=ωmn[1,0]/tpi

f01=ωmn[0,1]/tpi

f20=ωmn[2,0]/tpi
(* this checks whether the single-mode interval is indeed 15-22.5 GHz as stated above *)

B0=1;
m=1;
n=0;
(* we choose an arbitrary unit for the stationary field B_z[x,y], and specify the mode as TM_{1,0} *)

f0=1.25 f10
(* we dial a frequency on the generator within the single-mode TM_{1,0} operation interval - ignoring also TE modes! *)

ω0=tpi f0;
k0=Sqrt[ω0^2-ωmn[m,n]^2]/c
(* the wavenumber in inverse meters *)

λ=tpi/k0
(* the wavelength is in the cm range: *)

vGroup=c Sqrt[1-(f10/f0)^2]
(* the group velocity states how fast energy travels along z *)

Clear[BLong];
BLong[m_,n_]:=B0 Cos[m pi x/a]Cos[n pi y/b]

Bz=BLong[1,0]
(* in the TE_{1,0} mode the y-dependence in Bz is constant *)

Ez=0;
(* We now use Griffiths (9.180), which are the Maxwell curl equations to derive the transverse electric field components *)

Ex=I/((ω0/c)^2-k0^2)(k0 D[Ez,x]+ω0 D[Bz,y])

Ey=I/((ω0/c)^2-k0^2)(k0 D[Ez,y]-ω0 D[Bz,x])//Simplify

Bx=I/((ω0/c)^2-k0^2)(k0 D[Bz,x]-ω0/c^2 D[Ez,y])//Simplify

By=I/((ω0/c)^2-k0^2)(k0 D[Bz,y]+ω0/c^2 D[Ez,x])//Simplify

(* we have three field components that do not vanish: (i) Bz[x] only, i.e., no y-dependence in this mode; (ii) Ey[x] was derived; and (iii) Bx[x] was derived. For other modes we should expect dependence on y. *)
(* Now go for the real fields which depend on time as well: *)
EyR=Re[Ey Exp[I(k0 z-ω0 t)]]//ComplexExpand //Simplify

T=tpi/ω0

Clear[res];
res[z1_,t1_]:=(EyR/.x->a/2)/.{t->t1,z->z1}
(* Looking at x=a/2 where the field is strongest *)

Animate[Plot[res[z,t],{z,0,2λ}],{t,0,3T},AnimationRunning->False]

Clear[res1];
res1[x1_,t1_]:=(EyR/.z->0)/.{t->t1,x->x1}
(* Now fix z and animate the x dependence of Ey over time *)

Peak=MaxValue[res1[x,t],{x,t}]

Manipulate[Plot[res1[x,t],{x,0,a},PlotRange->{{0,a},{-Peak,Peak}}],{t,0,3T}(*,AnimationRunning->False*)]

(* Compute the real magnetic field components for visualization: *)

BxR=Re[Bx Exp[I(k0 z-ω0 t)]]//ComplexExpand //Simplify

BzR=Re[Bz Exp[I(k0 z-ω0 t)]]//ComplexExpand //Simplify

res2[z1_,x1_]:=BzR/.{t->T/2,z->z1,x->x1}
res3[z1_,x1_]:=BxR/.{t->T/2,z->z1,x->x1}

PL1=VectorPlot[{res2[z,x],res3[z,x]},{z,0,λ},{x,0,a},VectorStyle->Red,FrameLabel->{z,x}]

Plot[EyR/.{x->a/2,t->T/2},{z,0,λ}]
(* comparing this graph with the previous shows how the B field curls around the places where the electric field vanishes *)

PL2=ContourPlot[EyR/.t->T/2,{z,0,λ},{x,0,a},ContourLabels->All];
Show[PL2,PL1,FrameLabel->{z,x}]

(* Now assemble an animation by explicit plot construction: *)
nFrames=32;
tMax=2T;
PLa={};
PLb={};
Clear[res2z,res2x];
res2z[z1_,x1_,j_]:=BzR/.{t->j tMax/nFrames,z->z1,x->x1};
res2x[z1_,x1_,j_]:=BxR/.{t->j tMax/nFrames,z->z1,x->x1};
Do[
PLa=AppendTo[PLa,VectorPlot[{res2z[z,x,j],res2x[z,x,j]},{z,0,λ},{x,0,a},VectorStyle->Red]];
PLb=AppendTo[PLb,ContourPlot[EyR/.t->j tMax/nFrames,{z,0,λ},{x,0,a}]];
,{j,1,nFrames}];

ListAnimate[Table[Show[PLb[[j]],PLa[[j]],FrameLabel->{z,x}],{j,1,nFrames}],AnimationRunning->False]
(* watch for:
(i) the B-field pattern moves along with the equicontours of the electric field
(ii) look at the time evolution of the Bfield arrows for fixed location z.
 *)

(* Note that we have a magnetic field without y-dependence. This implies that the magnetic field is a parallel field to the surfaces y=0 and y=b, and we see from the above plots that it is also parallel to x=0 and x=a. This implies that there are surface currents in the waveguide material which are responsible for a non-zero magnetic field inside the guide, and vanishing magnetic field outside. The surface current vector K is calculated as the cross product of the area vector with the auxiliary field H. More accurately: apply Griffiths,Eq.(7.63): Hinside-Houtside = (K cross n-hat) at y=0 proceeding from Houtside=0 towards y>0. This implies that K = y-hat cross H. 
Looking at the surface current in the panel at y=0; presumably the reverse occurs at y=b;
the area vector direction is y-hat.
Since we are not concerned with units here, we drop the μ0.
*)
Cross[{0,1,0},{Bvx,0,Bvz}]
Cross[{1,0,0},{Bvx,0,Bvz}]

PL3=StreamDensityPlot[{-res3[z,x],res2[z,x]},{z,0,λ},{x,0,a}(*,VectorStyle->Blue*),FrameLabel->{z,x}];

Show[PL3,PL1]
(* to reconcile this plot with the explanation above, note that the y-axis is out of the screen for the chosen x and z axes. This implies we are stepping from behind the screen into the waveguide towards y>0. By the right-hand rule the surface current 'turns on' the H field. At the y=b interface the opposite surface current direction is turning off the H field as one steps out of the waveguide towards y>b.

Note: darker (brown) shading means less current, brighter (white) = stronger current.
 *)

(* To confirm the conclusion that at y=b the surface current has the opposite direction, which follows from Eq.(7.63), since now H_above=H_outside=0, and Hbelow=H_inside, and thus, a minus sign when applying the cross product with n-hat=y-hat, we shoudl look at the surface current on the panels y=0, y=b. The x-z plot above shows that there are strong fields with alternating direction as one moves along z which exist throughout those panels. The cross product calculated above, n_hat=x_hat cross H *)

PL4=VectorPlot[{res2[z,0],0},{z,0,λ},{y,0,b},VectorStyle->Red,FrameLabel->{z,y}];
PL5=StreamDensityPlot[{0,-res2[z,0]},{z,0,λ},{y,0,b}(*,VectorStyle->Blue*)];
Show[PL5,PL4,FrameLabel->{z,y}]

(* For x=a we have the reverse situation: *)
PL6=VectorPlot[{res2[z,a],0},{z,0,λ},{y,0,b},VectorStyle->Red,FrameLabel->{z,y}];
PL7=StreamDensityPlot[{0,-res2[z,a]},{z,0,λ},{y,0,b}(*,VectorStyle->Blue*)];
Show[PL7,PL6,FrameLabel->{z,y}]