#############################################################################
# This script combines the results from the sweep
# and plots the frequency spectrum. It then obtains effective index and beta.
# The bandstructure information is extracted using tolerance and num_band
# specified by the user at the beginning of the script.
#
# Properties:     
# f_band:       Frequencies of bands in units of Hz
#############################################################################
# User Defined properties:
tolerance = 1e-2; #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 = 2;    #number of bands to search for in the bandstructure

#############################################################################

runsweep; # run all sweeps


# get fd data from the sweeps
sweepname="sweep kx";
spectrum = getsweepresult(sweepname,"spectrum");
fd = pinch(spectrum.fd);
f = spectrum.f;
beta = spectrum.kx;

# simple imaging of fd vs k
image(c/f*1e6,beta,fd,"lambda (microns)","beta (1/m)","beta","logplot");
# or plot
plot(f*1e-12,fd,"f THz","spectrum");

# extract the bandstructure
bandstructure=matrix(num_band,length(beta)); # initialize matrix in which to store band frequency information

# loop over sweep results
for (i=1:length(beta)){
  #use findpeaks to find num_band number of peaks
  temp = findpeaks(fd(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 = fd(temp(1),i)*tolerance;
  f_band=matrix(num_band);
  for(bandcount = 1:num_band) {
    if( fd(temp(bandcount),i) > minvalue) { 
      f_band(bandcount) = f(temp(bandcount)); 
    }
  }

bandstructure(1:num_band,i)=f_band;
}

bandstructure=transpose(bandstructure);  

# exctract neff vs f from the bandstructure
f1 = pinch(bandstructure,2,1);
f2 = pinch(bandstructure,2,2);
# resort beta (assume no mode crossings)

for(i=1:length(f1)) {
  	if(f1(i) > f2(i)) {
    	temp = f2(i);
    	f2(i) = f1(i);
    	f1(i) = temp;
  	}

}

neff2 = beta / ( 2*pi*f2/c);
if(max(abs(f1)) > 0) {
  neff1 = beta / ( 2*pi*f1/c);
} else {
  neff1 = 0*neff2;
}

# set neff2 to neff1 for values where 2 peaks were not found
neff1 = (neff1==0)*neff2 + (neff1!=0)*neff1;


plotxy(c/f1*1e6,beta,c/f2*1e6,beta,"wavelength (microns)","beta","");
legend("beta1","beta2");
plotxy(c/f1*1e6,neff1,c/f2*1e6,neff2,"wavelength (microns)","neff","");
legend("neff1","neff2");