# Garcia_vidal.lsf
##################################################
# Plots E field intensity as function of x and y. 
# Finds the focal length (distance from surface to the maximum E field intensity.
# Finds the focal width (FWHM at maximum E field intensity along the y direction)
####################################################

# choose area to plot and spatial resolution
nx=500;#to dupilcate online figure, use 1000
ny=250;#to dupilcate online figure, use 500
x=linspace(0.2e-6,150e-6,nx);
y=linspace(-20e-6,20e-6,ny);

# do far field projection
E2 = farfieldexact2d('monitor1',x,y);
E2 = sum(abs(E2)^2,3);

# find the focal length
y0=find(y,0);
xp = findpeaks(pinch(E2,2,y0));

?"Focal length: " + num2str(x(xp)*1e6) + " microns";

# find the focal width
m = E2(xp,y0);
cs_y = pinch(E2,1,xp)/m;

for(i=y0;cs_y(i)>.5;i=i+1){1;}
focal_width = y(i-1);
for(i=y0;cs_y(i)>.5;i=i-1){1;}
focal_width = (focal_width - y(i+1))*1e6;

?"Focal width: " + num2str(focal_width) + " microns";

# plot E field intensity profile as a function of x,y
image(x*1e6,y*1e6,E2,'x (microns)','y (microns)','E field intensity');
setplot("colorbar max",m);

# Cut along the line y = focal_length
plot(x*1e6,pinch(E2,2,y0)/m,'x (microns)','E field intensity (au)','');
setplot('y max',1);

# plot cut along the line x = 0
plot(y*1e6,cs_y,'y (microns)','E field intensity (au)','x='+num2str(x(xp)*1e6)+" microns");