#############################################################################
# This script performs the Gamma-X sweep, plots the bandstructure and dispersion
# for waveguide3D.fsp
#
# Properties:     
# a:            period used to normalize the frequency (fs_all_norm=fs_all*a/c)
#               in this case the hexagonal lattice constant
# fs_all:       Frequencies of bands in units of Hz
# fs_all_norm:  Frequencies of bands in units of Hz * a / c
#############################################################################

runsweep; # run Gamma-X sweep

# get fs data from the sweep
sweepname="Gamma-X";
kx=getsweepdata(sweepname,"kx");
resonance=getsweepresult(sweepname,"fs");
fs_all=resonance.fs;

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

# calculate the lightline for n=1.45
period_x=a*sqrt(3);
f_lightline = c*kx*(2*pi/period_x)/(2*pi)/1.45*a/c;
f_lightline_upper = c*(1-kx)*(2*pi/period_x)/(2*pi)/1.45*a/c;

fs_all_norm = fs_all*a/c;

# plot bandstructure 
plot(kx,transpose(fs_all_norm),"k (Gamma-X)","f (Hz*ax/c)","Bandstructure","plot points, color=blue");
plot(kx,transpose(fs_all_norm),f_lightline,f_lightline_upper,
     "kx (normalized)","f (normalized)",
     "Bandstructure & Lightline","plot points, color=blue");

# if there is only one mode, do group velocity calculation
sz=size(fs_all);
if (sz(2)==1) { 
  ?"Calculating dispersion.";

  # get data from sweep 
  slope=pinch(getsweepdata(sweepname,"log_slope"));

  n_bands=sz(2);   # number of bands
  p=sz(1);         # number of k points

  vg = matrix(p);
  for(i=1:p) {
    if(i == p) {
      f2 = fs_all(i,1);
      k2 = kx(i)*2*pi/period_x;
    } else {
      f2 = fs_all(i+1,1);
      k2 = kx(i+1)*2*pi/period_x;
    }
    if(i == 1) {
      f1 = fs_all(i,1);
      k1 = kx(i)*2*pi/period_x;	    
    } else {
      f1 = fs_all(i-1,1);
      k1 = kx(i-1)*2*pi/period_x;
    }
    vg(i) = abs(2*pi*(f2-f1)/(k2-k1));
  }
  loss = abs(10*log10(exp(1))*slope/vg);

  plot(kx,vg/c,"kx (units of 2*pi/a)","group velocity (vg/c)","group velocity vs kx");
  plot(kx,loss*1e-3,"kx (units of 2*pi/a)","loss dB/mm","loss vs kx");	
}