#############################################################################
# This script is based on the bandstructure example files
# https://kb.lumerical.com/en/index.html?diffractive_optics_pc_bandstructure.html
#############################################################################

runsweep;

# get fs data from the sweeps
sweepname="ky_sweep";
f=pinch(getsweepdata(sweepname,"f"),2,1); # get array of frequencies
fs= pinch(getsweepdata(sweepname,"fs"));
ky=getsweepdata(sweepname,"ky");

# simple imaging of fs vs k
image(ky,f/1e12,transpose(fs),"ky","f (THz)","bandstructure, logscale","logplot");

# plot bandstructure
num_band = 10;
tolerance = 1.5e-4;
bandstructure=matrix(num_band,length(ky)); # initialize matrix in which to store band frequency information

# extract resonance frequency for each k value
for (i=1:length(ky)){
  #use findpeaks to find num_band number of peaks
  temp = findpeaks(fs(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(temp(1),i)*tolerance;
  f_band=matrix(num_band);
  for(bandcount = 1:num_band) {
    if( fs(temp(bandcount),i) > minvalue) { 
      f_band(bandcount) = f(temp(bandcount)); 
    }
  }

#f_band_norm = f_band*a/c;  # normalize the frequency vector
bandstructure(1:num_band,i)=f_band;
}

bandstructure=transpose(bandstructure);  
plot(ky,bandstructure/1e12,"ky","f (THz)","bandstructure","plot points");