#############################################################################
# This script gets the results from the Gamma-X sweep
# and plots the frequency spectrum fs from the bandstructure object
# results over k, and the bandstructure.
# The bandstructure information is extracted using tolerance and num_band
# specified by the user at the beginning of the script.
#
# Properties:     
# a:            period used to normalize the frequency (f_norm=f*a/c)
# f_band:       Frequencies of bands in units of Hz
# f_band_norm:  Frequencies of bands in units of Hz * a / c
#############################################################################
clear;
# User Defined properties:
tolerance = 1e-3; #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
#
#############################################################################
nsweep=10;
runsweep; # run Gamma-X sweep

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

# get fs data from the sweeps
sweepname="Gamma-X";
spectrum=getsweepresult(sweepname,"spectrum");
f=c/spectrum.lambda;
fs=spectrum.fs;
nfs=size(fs);
nn=nfs(1);
# simple imaging of fs vs k
plot(f*1e-12,transpose(fs(1:nn,nsweep)),"f (THz)","Amplitude (au)");
plot(f*1e-12,transpose(fs(1:nn,nsweep)),"f (THz)","Amplitude (au)","","logplot");
image(1:nsweep,f,transpose(fs),"k (Gamma-X)","f (Hz)","bandstructure, logscale","logplot");
image(1:nsweep,f,transpose(fs),"k (Gamma-X)","f (Hz)","bandstructure, linearscale");
setplot("colorbar min",0);
setplot("colorbar max",max(fs)*1e-4);

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

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

bandstructure=transpose(bandstructure);  
plot(1:nsweep,bandstructure,"k (Gamma-X)","f (Hz*a/c)","bandstructure","plot points");
bandstructure_f=bandstructure*c/a;
plot(1:nsweep,bandstructure_f*1e-12,"k (Gamma-X)","f (THz)","bandstructure","plot points");