# define desired region of x and y
x = linspace(-5e-6,5e-6,200);
y = linspace(1e-6,50e-6,500);

# do far field projection
E2 = farfieldexact2d('T',x,y); 
E2 = sum(abs(E2)^2,3); # calculate E2 from Ex, Ey, Ez
# plot E field intensity profile as a function of x,y
image(x*1e6,y*1e6,E2,'x (microns)','y (microns)','E field intensity');

# calculate Hz
delta = 1e-9; # used to calculate the numerical derivative;
f = getdata("T","f");
Ey1 = pinch(farfieldexact2d("T",x-delta,y),3,2);
Ey2 = pinch(farfieldexact2d("T",x+delta,y),3,2); 
dEy_dx = (Ey2-Ey1)/(2*delta);

Ex1 = pinch(farfieldexact2d("T",x,y-delta),3,1);
Ex2 = pinch(farfieldexact2d("T",x,y+delta),3,1); 
dEx_dy = (Ex2-Ex1)/(2*delta);

# plot H field intensity profile as a function of x,y
Hz = (dEy_dx - dEx_dy) / (1i * 2* pi * f * mu0);
image(x*1e6,y*1e6,abs(Hz)^2,"x (um)","y (um)","Hz");