# convergence with number of cells in tapered region
if(1){
  switchtolayout;

  # reset values of other settings in case convergence testing was
  # performed previously for other settings
  setnamed("EME","number of modes for all cell groups",10);
  setnamed("mesh","dy",0.01e-6);
  setnamed("mesh","dz",0.01e-6);
  
  npts = 10; # number of points to sweep
  ncells = linspace(2,20,npts); # values to sweep

  S52 = matrix(npts); # initialize result matrix

  for(i = 1:npts){ # sweep
    switchtolayout;
    setnamed("EME","cells",[1;ncells(i);1]);
    run;
    emepropagate;
    S = getresult("EME","user s matrix");
    S52(i) = S(2,1);
  }

  plot(ncells,abs(S52)^2,"Number of cells in tapered region",
	"|S52|^2");

  # delta_T
  delta_T = (abs(S52(2:npts))^2-abs(S52(1:npts-1))^2)^2/abs(S52(2:npts))^4;

  # delta_T_N
  delta_T_N = (abs(S52)^2-abs(S52(npts))^2)^2/abs(S52(npts))^4;

  # plot delta_T and delta_T_N
  plotxy(ncells(2:npts),delta_T, ncells,delta_T_N,
         "Number of cells in tapered region","delta_T");
  legend("delta_T","delta_T_N");
}

# the convergence with number of eigenmodes
if(1){
  switchtolayout;
  
  # reset values of other settings in case convergence testing was
  # performed previously for other settings

  setnamed("mesh","dy",0.01e-6);
  setnamed("mesh","dz",0.01e-6);
  
  start_mode = 4; #set smaller number of modes for convergence test
  mode_interval = 1; #set mode interval for convergence test

  setnamed("EME","number of modes for all cell groups",25); #should be the largest number of modes desired for test    

  run;
  setemeanalysis("mode convergence sweep", 1);
  setemeanalysis("start mode", start_mode); 
  setemeanalysis("mode interval", mode_interval); 
  emesweep("mode convergence sweep");
  
  S = getemesweep("S_mode_convergence_sweep");
  nmodes = S.modes;
  S52 = S.s52;
  npts = length(S52);
  
  plot(nmodes,abs(S52)^2,"Number of eigenmodes","|S52|^2");
  
  # delta_T
  delta_T = (abs(S52(2:npts))^2-abs(S52(1:npts-1))^2)^2/abs(S52(2:npts))^4;

  # delta_T_N
  delta_T_N = (abs(S52)^2-abs(S52(npts))^2)^2/abs(S52(npts))^4;

  # plot delta_T and delta_T_N
  plotxy(nmodes(2:npts),delta_T, nmodes,delta_T_N,
         "Number of eigenmodes","delta_T");
  legend("delta_T","delta_T_N");

}

# convergence with transverse mesh in mesh override region
if(1){
  switchtolayout;
  
  # reset values of other settings in case convergence testing was
  # performed previously for other settings
  setnamed("EME","number of modes for all cell groups",10);
  setnamed("EME","cells",[1;19;1]);

  npts = 10; # number of points to sweep
  step_size = linspace(0.05e-6,0.01e-6,npts); # values to sweep

  S52 = matrix(npts); # initialize result matrix

  for(i = 1:npts){ # sweep
    switchtolayout;
    setnamed("mesh","dy",step_size(i));
    setnamed("mesh","dz",step_size(i));
    run;
    emepropagate;
    S = getresult("EME","user s matrix");
    S52(i) = S(5,2);
  }

  plot(step_size,abs(S52)^2,"mesh step size","|S52|^2");

  # delta_T
  delta_T = (abs(S52(2:npts))^2-abs(S52(1:npts-1))^2)^2/abs(S52(2:npts))^4;

  # delta_T_N
  delta_T_N = (abs(S52)^2-abs(S52(npts))^2)^2/abs(S52(npts))^4;

  # plot delta_T and delta_T_N
  plotxy(step_size(2:npts),delta_T, step_size,delta_T_N,
         "Mesh step size in mesh override region","delta_T");
  legend("delta_T","delta_T_N");
}