clear;
switchtolayout;
select("Gaussian-Source");
#Disable the multifrequency calculation for the first run
set("multifrequency beam calculation",0);
run;

# get monitor data
x=getdata("Fields-Monitor","x");
y=getdata("Fields-Monitor","y");
f=getdata("Fields-Monitor","f");
Ez=pinch(getdata("Fields-Monitor","Ez"));

#plot field profile at min and max frequency
fi=1;           # first frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency Off");

fi=ceil(length(f)/2);   # center frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency Off");

fi=length(f);   # last frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency Off");





clear;
switchtolayout;
select("Gaussian-Source");
#Enable the multifrequency calculation for the second run
set("multifrequency beam calculation",1);
run;

# get monitor data
x=getdata("Fields-Monitor","x");
y=getdata("Fields-Monitor","y");
f=getdata("Fields-Monitor","f");
Ez=pinch(getdata("Fields-Monitor","Ez"));

#plot field profile at min and max frequency
fi=1;           # first frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency On");

fi=ceil(length(f)/2);   # center frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency On");

fi=length(f);   # last frequency
Ez_real=real(pinch(Ez,3,fi));
image(x*1e6,y*1e6,Ez_real,"x (um)","y (um)","Ez at "+num2str(f(fi)/1e12)+"THz Multi-Frequency On");