clear;
closeall;
run;
# to get the monitor data
Ez = pinch(getdata("monitor1","Ez"));
f=getdata("monitor1","f");
# prepare the data for fft
Ez1 = pinch(Ez,3,1);  #lambda=0.5 um
Ez2 = pinch(Ez,3,2);  #lambda=0.75 um
Ez3 = pinch(Ez,3,3);  #lambda=1 um
Ez4 = pinch(Ez,3,4);  #lambda=1.25 um
Ez5 = pinch(Ez,3,5);  #lambda=1.5 um
#
x = getdata("monitor1","x");
y = getdata("monitor1","y");
nx=length(x);
ny=length(y);
X=meshgridx(x,y);
Y=meshgridy(x,y);
filt=cos(pi*X/(max(x)-min(x)))*cos(pi*Y/(max(y)-min(y)));
# plot spatial fields
image(x*1e6,y*1e6,abs(Ez5)^2,"x (um)","y (um)","lambda=1.5um","logplot");


# period, ax is the lattice constant
a=0.325;lambda=1e-6;ax=a*lambda;
nczt=2*120;
kx = linspace(-3*pi/ax,3*pi/ax,nczt);
ky = kx;
Ezf1 = czt(Ez1,x,y,kx,ky);
image(kx*1e-6,ky*1e-6,abs(Ezf1)^2,"kx 1/um","ky 1/nm","lambda=0.5um EFC","logplot");
Ezf10 = czt(Ez1*filt,x,y,kx,ky);
image(kx*1e-6,ky*1e-6,abs(Ezf10)^2,"kx 1/um","ky 1/nm","lambda=0.5um EFC","logplot");
#using fft;fft will use 1028 data points
# use the central k region 200 points
kn=(512-100):512+100;
kxf = fftk(x,3);
kyf = fftk(y,3);
Efft1 = fft(Ez1*filt,3);Eff1=  Efft1/max(abs(Efft1));
Eff10=abs(Eff1)>1/3;
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff1(kn,kn))^2,"kx 1/um","ky 1/nm","Lambda=0.5um");

Efft2 = fft(Ez2*filt,3);Eff2= Efft2/max(abs(Efft2));
Eff20=abs(Eff2)>1/3;
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff2(kn,kn))^2,"kx 1/um","ky 1/nm","Lambda=0.75um");
# you can plot Eff30 or Eff3
Efft3 = fft(Ez3*filt,3);Eff3= Efft3/max(abs(Efft3));
Eff30=abs(Eff3)>1/3;
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff30(kn,kn))^2,"kx 1/um","ky 1/nm","Lambda=1.0um");

Efft4 = fft(Ez4*filt,3);Eff4= Efft4/max(abs(Efft4));
Eff40=abs(Eff4)>1/3;
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff40(kn,kn))^2,"kx 1/um","ky 1/nm","Lambda=1.25um");

Efft5 = fft(Ez5*filt,3);Eff5= Efft5/max(abs(Efft5));
Eff50=abs(Eff5)>1/3;
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff5(kn,kn))^2,"kx 1/um","ky 1/nm","Lambda=1.5um");
setplot("colorbar min",0.4);
setplot("colorbar max",0.45);

# f4 and f5 are in higher bands, the graph is not good since the coupling is not good
# getting the 3 contours into 1 graph
image(kxf(kn)*1e-6,kyf(kn)*1e-6,abs(Eff30(kn,kn))-abs(Eff40(kn,kn))+0.4*abs(Eff50(kn,kn)),"kx 1/um","ky 1/nm","EFCs ");
# to see higher Brillouin zone
kn1=(512-150):512+150;
image(kxf(kn1)*1e-6,kyf(kn1)*1e-6,(abs(Eff5(kn1,kn1)))^0.01,"kx 1/um","ky 1/nm","Brillouin Zone,Lambda=1.5um");
setplot("colorbar min",0.93);
setplot("colorbar max",1);