###########################################################################
# Scriptfile: vector_pupil_example.lsf
#
# Description:
# This example script file shows how to create beams with radial/azimuthal polarization
# using a vectorial pupil function.
# 
###########################################################################

clear;
closeall;
switchtolayout;

# Settings:
source_name = "source";
T_monitor = "T";
source_theta = 0; # azimuthal angle
source_phi = 0; # polar angle
NA = 0.2; # maximum NA of the beam
pol = 2; # 1 for radial, 2 for azimuthal polarization
n = 501; # number of points in each direction of k-space

#Define pupil function:
ux = linspace(-NA, NA, n);
uy = linspace(-NA, NA, n);
Ux = meshgridx(ux, uy);
Uy = meshgridy(ux, uy);
Rho = sqrt(Ux^2+Uy^2);
p = (Rho <= NA);
theta = atan2(Uy, Ux);

if (pol == 1) {
    ## Radial polarization
    E1 = cos(theta)*p;
    E2 = sin(theta)*p;
} else {
    if (pol == 2){
        ## Azimuthal polarization
        E1 = -sin(theta)*p;
        E2 = cos(theta)*p;
    } else{ ?"Error: polarization option not defined"; break; }
}


pupil = matrixdataset("pupil");
pupil.addparameter("ux",ux);
pupil.addparameter("uy",uy);
pupil.addattribute("E1", E1);
pupil.addattribute("E2", E2);
visualize(pupil);
matlabsave("vector_pupil.mat",pupil); # create matlab file that can be imported in source using GUI

# Apply pupil function
select("source");
set("angle theta", source_theta);
set("angle phi", source_phi);
set("use custom pupil function",true);
importdataset(pupil); # import directly in source

#Run simulation
run;

# Plot fields recorded by a monitor in front of the source
E_forward = getresult(T_monitor,"E");
visualize(E_forward);
vectorplot(E_forward); # use vector plot to identify polarization

# Plot radial and azimuthal components of E field recorded by monitor:
Ex = pinch(getdata(T_monitor, "Ex"));
Ey = pinch(getdata(T_monitor, "Ey"));
Ez = pinch(getdata(T_monitor, "Ez"));
x = getdata(T_monitor, "x");
y = getdata(T_monitor, "y");

X = meshgridx(x,y);
Y = meshgridy(x,y);
PHI = atan2(Y,X);
Eradial = Ex*cos(PHI) + Ey*sin(PHI);
Ephi = -Ex*sin(PHI) + Ey*cos(PHI);
image(x*1e6,y*1e6,abs(Eradial)^2,"x (microns)","y (microns)","Intensity Eradial");
image(x*1e6,y*1e6,abs(Ephi)^2,"x (microns)","y (microns)","Intensity Ephi");

# Plot farfield:
ff = farfield3d(T_monitor);
ux_ff = farfieldux(T_monitor);
uy_ff = farfielduy(T_monitor);
polarimage(ux_ff,uy_ff,ff,"ux","uy","Far field result"); # Far field