#############################################################################
# This script performs the Gamma-X, X-M, and M-Gamma sweeps
# and plots the bandstructure and dispersion for the all-dielectric
# zero-index metamaterial in the FDTD file dielectric_zim_band.fsp.
# At the end of the sweeps, the script finds the field profiles of
# the electric and magnetic dipole Mie resonances within the Silicon rods.
#
# Properties:     
# a:            period used to normalize the frequency (f_norm=f*a/c)
#               in this case the hexagonal lattice constant
# f_band:       Frequencies of bands in units of Hz
# f_band_norm:  Frequencies of bands in units of Hz * a / c
# Copyright 2017 Lumerical Solutions Inc
#############################################################################
# User Defined properties:
tolerance = 8e-5; #tolerance for finding peaks and accepting bands
#                 setting this too low will result in noisy data where sidelobes of
#                 peaks are interpreted as new bands
#                 setting it too high will mean that some bands are not found 
num_band =10;     #number of bands to search for in the bandstructure
#
#############################################################################

#Turn on dipole cloud.
switchtolayout;
select("dipole_cloud");set("enabled",1);
select("single_dipole");set("enabled",0);
runsweep; # run all sweeps

# get fs data from the sweeps
sweepname="Gamma - X";
spectrum_g_x=getsweepresult(sweepname,"spectrum");
sweepname="X - M";
spectrum_x_m=getsweepresult(sweepname,"spectrum");
sweepname="M-Gamma";
spectrum_m_g=getsweepresult(sweepname,"spectrum");

#get size for data
s1=size(spectrum_g_x.fs);s2=size(spectrum_x_m.fs);s3=size(spectrum_m_g.fs);
st=s1(2)+s2(2)+s3(2);
a = getnamed("::model","ax");
f=c/spectrum_g_x.lambda;
#construct bandstructure
fs_all=matrix(length(f),st); # initialize matrix to store fs data in
fs_all(1:length(f),1:s1(2))=spectrum_g_x.fs;
fs_all(1:length(f),(s1(2)+1):(s1(2)+s2(2)))=spectrum_x_m.fs;
fs_all(1:length(f),(s1(2)+s2(2)+1):st)=spectrum_m_g.fs;

# construct bandstructure matrix
bandstructure=matrix(num_band,st); # initialize matrix in which to store band frequency information

# loop over sweep results
for (i=1:st){
  #use findpeaks to find num_band number of peaks
  temp = findpeaks(fs_all(1:length(f),i),num_band);

  #collect data for any peaks that are more than 'tolerance' of the maximum peak (to avoid minor peaks like sidelobes)
  minvalue = fs_all(temp(1),i)*tolerance;
  f_band=matrix(num_band);
  for(bandcount = 1:num_band) {
    if( fs_all(temp(bandcount),i) > minvalue) { 
      f_band(bandcount) = f(temp(bandcount)); 
    }
  }
f_band_norm = f_band;  # normalize the frequency vector
bandstructure(1:num_band,i)=f_band_norm;
}
bandstructure=transpose(bandstructure);  
#plot band structure
plot(1:st,bandstructure/1e12,"k (Gamma-X-M-Gamma)","f (THz)","bandstructure","plot points, color=blue");
setplot("show legend",0);

#################################################################################
#Turn off dipole cloud and turn on single dipole source
switchtolayout;
select("dipole_cloud");set("enabled",0);
select("single_dipole");set("enabled",1);

#Electric dipole
set("dipole type", "electric");set("phi",0);set("theta",0);run;
E=getresult("XY Monitor","E");
f=E.f;ind_f=find(f/1e12,211);
Ez=pinch(E.Ez,4,ind_f);
image(E.x/1e-6,E.y/1e-6,-real(Ez));setplot("title","Electric Dipole Ez");

#Magnetic Dipole
switchtolayout;
select("single_dipole");set("dipole type", "magnetic");set("phi",90);set("theta",90);run;
H=getresult("XY Monitor","H");
f=H.f;ind_f=find(f/1e12,211);
Hy=pinch(H.Hy,4,ind_f);
image(H.x/1e-6,H.y/1e-6,real(Hy));setplot("title","Magnetic Dipole Hy");
setplot("colorbar min", -7);setplot("colorbar max", 8);