#########################################################
# substrate_integrated_waveguide.lsf
# 
# This script performs the Gamma-X sweep in FDTD and plots 
# the bandstructure and attenuation loss of the SIW's
# fundamental TE10 mode. If multiple modes
# Uses FDTD to calculates the return loss, directivity 
# pattern, max directivity, gain, and radiation efficiency
# of a rectangular open-ended waveguide (OEWG) probe
# antenna mounted on a finite and infinite ground plane.
#
# Copyright 2016 Lumerical Solutions Inc
##########################################################
#############################################################################
# 
#
# 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 = 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 =1;    #number of bands to search for in the bandstructure
#
#############################################################################
f1=8e9;f2=24.5e9;
switchtolayout;
select("::model");set("f1",f1);set("f2",f2);
setglobalsource("frequency start", f1);setglobalsource("frequency stop",f2);
runsweep; # run Gamma-X sweepx

# get fs data from the sweep
sweepname="Gamma - X";
ky=getsweepdata(sweepname,"ky");
spectrum=getsweepresult(sweepname,"spectrum");
slope=pinch(getsweepdata(sweepname,"log_slope"));
fields=getsweepdata(sweepname,"fields");
f=c/spectrum.lambda;
fs_all=spectrum.fs;

# get a from model
S = getnamed("::model::SIW","S");
period_y=S;

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

# loop over sweep results
for (i=1:length(ky)){
  #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*S/c;  # normalize the frequency vector
norm_bandstructure(1:num_band,i)=f_band;#_norm;
}
norm_bandstructure=transpose(norm_bandstructure);  

sz=size(bandstructure);
p=sz(2);         # number of k points

vg = matrix(p);
  for(i=1:p) {
    if(i == p) {
      f2 = bandstructure(1,i);
      k2 = ky(i)*2*pi/period_y;
    } else {
      f2 = bandstructure(1,i+1);
      k2 = ky(i+1)*2*pi/period_y;
    }
    if(i == 1) {
      f1 = bandstructure(1,i);
      k1 = ky(i)*2*pi/period_y;	    
    } else {
      f1 = bandstructure(1,i-1);
      k1 = ky(i-1)*2*pi/period_y;
    }
    vg(i) = abs(2*pi*(f2-f1)/(k2-k1));
  }
  loss = 10*log10(exp(1))*slope/vg;
  
  plot(ky*2*pi/S,norm_bandstructure/1e9,"Wavenumber (1/m)","Frequency (GHz)","Dispersion Diagram","plot points");legend("TE10");
  plot(norm_bandstructure/1e9,loss/1e3,"Frequency (GHz)","loss dB/mm","loss vs Frequency");legend("TE10");