#############################################################################
# This script performs the Gamma-kz sweep, plots the bandstructure and dispersion
# for tri_2D_kz.fsp
#
# 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
#############################################################################
# User Defined properties:
tolerance = 1e-4; #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 = 15;    #number of bands to search for in the bandstructure
#
#############################################################################

runsweep; # run sweep

# get fs data from the sweep
sweepname="Gamma-kz";
kz=getsweepdata(sweepname,"kz");
spectrum=getsweepresult(sweepname,"spectrum");
f=c/spectrum.lambda;  # get array of frequencies
fs_all=spectrum.fs;

# get a from model
a = getnamed("::model::hex_pc","a");

# simple imaging of fs vs k
image(kz,f*a/c,transpose(fs_all),"k (Gamma-kz)","f (Hz*a/c)","bandstructure, logscale","logplot");
image(kz,f*a/c,transpose(fs_all),"k (Gamma-kz)","f (Hz*a/c)","bandstructure, linearscale");
setplot("colorbar min",0);
setplot("colorbar max",max(fs_all)*1e-4);

# extract bandstructure
bandstructure=matrix(num_band,length(kz)); # initialize matrices in which to store band frequency information
norm_bandstructure=matrix(num_band,length(kz)); 

# loop over sweep results
for (i=1:length(kz)){
  #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)); 
    }
  }

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

# plot bandstructure
norm_bandstructure=transpose(norm_bandstructure);  
plot(kz,norm_bandstructure,"k (Gamma-kz)","f (Hz*a/c)","bandstructure","plot points");