#######################################################
# File: usr_convergence_analysis.lsf
# 
# Description: This file can be used to analyze the 
#  results of parameter sweeps in usr_convergence.fsp
#
# Copyright 2011, Lumerical Solutions, Inc.
#######################################################

# the convergence with PML distance
if(true) {
  # get the data
  span = getsweepdata("sweep_PML_distance","x_span");
  f = getsweepdata("sweep_PML_distance","f");
  s_scatter = getsweepdata("sweep_PML_distance","scattered");
  s_abs = -getsweepdata("sweep_PML_distance","absorbed");
  s_ext = pinch(s_scatter+s_abs);

  # plot the spectra
  plot(c/f*1e9,s_ext*1e9,"wavelength (nm)","extinction cross section (nm)");

  # delta_s
  delta_s_ext = s_ext(1:length(f),2:length(span))-s_ext(1:length(f),1:length(span)-1);
  delta_s_ext = sqrt( integrate( delta_s_ext^2, 1, c/f) / integrate( s_ext(1:length(f),2:length(span))^2, 1, c/f) );
  
  # delta_s_N
  s_ext_N = meshgridx(pinch(s_ext,2,length(span)),span);
  delta_s_ext_N = s_ext_N-s_ext;
  delta_s_ext_N = sqrt( integrate( delta_s_ext_N^2, 1, c/f) / integrate( s_ext_N^2, 1, c/f) );

  # plot delta_s and delta_s_N
  plotxy(span(2:length(span))*1e6,delta_s_ext,
         span*1e6,delta_s_ext_N,
         "simulation span (microns)","delta_s");
  legend("delta_s","delta_s_N");
}

# the convergence with PML layers
if(true) {
  # get the data
  layers = getsweepdata("sweep_PML_layers","Layers");
  f = getsweepdata("sweep_PML_layers","f");
  s_scatter = getsweepdata("sweep_PML_layers","scattered");
  s_abs = -getsweepdata("sweep_PML_layers","absorbed");
  s_ext = pinch(s_scatter+s_abs);

  # plot the spectra
  plot(c/f*1e9,s_ext*1e9,"wavelength (nm)","extinction cross section (nm)");

  # delta_s
  delta_s_ext = s_ext(1:length(f),2:length(layers))-s_ext(1:length(f),1:length(layers)-1);
  delta_s_ext = sqrt( integrate( delta_s_ext^2, 1, c/f) / integrate( s_ext(1:length(f),2:length(layers))^2, 1, c/f) );

  # delta_s_N
  s_ext_N = meshgridx(pinch(s_ext,2,length(layers)),layers);
  delta_s_ext_N = s_ext_N-s_ext;
  delta_s_ext_N = sqrt( integrate( delta_s_ext_N^2, 1, c/f) / integrate( s_ext_N^2, 1, c/f) );

  # plot delta_s and delta_s_N
  plotxy(layers(2:length(layers)),delta_s_ext,
         layers,delta_s_ext_N,
         "PML layers","delta_s");
  legend("delta_s","delta_s_N");
}

# the convergence with mesh accuracy
if(true) {
  # get the data
  mesh = getsweepdata("sweep_mesh_accuracy","mesh_accuracy");
  f = getsweepdata("sweep_mesh_accuracy","f");
  s_scatter = getsweepdata("sweep_mesh_accuracy","scattered");
  s_abs = -getsweepdata("sweep_mesh_accuracy","absorbed");
  s_ext = pinch(s_scatter+s_abs);

  # plot the spectra
  plot(c/f*1e9,s_ext*1e9,"wavelength (nm)","extinction cross section (nm)");

  # delta_s
  delta_s_ext = s_ext(1:length(f),2:length(mesh))-s_ext(1:length(f),1:length(mesh)-1);
  delta_s_ext = sqrt( integrate( delta_s_ext^2, 1, c/f) / integrate( s_ext(1:length(f),2:length(mesh))^2, 1, c/f) );

  # delta_s_N
  s_ext_N = meshgridx(pinch(s_ext,2,length(mesh)),mesh);
  delta_s_ext_N = s_ext_N-s_ext;
  delta_s_ext_N = sqrt( integrate( delta_s_ext_N^2, 1, c/f) / integrate( s_ext_N^2, 1, c/f) );

  # plot delta_s and delta_s_N
  plotxy(mesh(2:length(mesh)),delta_s_ext,
         mesh,delta_s_ext_N,
         "mesh accuracy setting","delta_s");
  legend("delta_s","delta_s_N");

}

# the convergence with inner dx, including comparison with Mie theory
if(true) {
  # get the data
  dx = getsweepdata("sweep_inner_mesh","dx");
  mesh_step = getsweepdata("sweep_inner_mesh","mesh_step");
  max_step = getsweepdata("sweep_inner_mesh","max_steps");
  f = getsweepdata("sweep_inner_mesh","f");
  s_scatter = getsweepdata("sweep_inner_mesh","scattered");
  s_abs = -getsweepdata("sweep_inner_mesh","absorbed");
  s_ext = pinch(s_scatter+s_abs);

  # plot the spectra
  plot(c/f*1e9,s_ext*1e9,"wavelength (nm)","extinction cross section (nm)");

  # delta_s
  delta_s_ext = s_ext(1:length(f),2:length(mesh_step))-s_ext(1:length(f),1:length(mesh_step)-1);
  delta_s_ext = sqrt( integrate( delta_s_ext^2, 1, c/f) / integrate( s_ext(1:length(f),2:length(mesh_step))^2, 1, c/f) );

  # delta_s_N
  s_ext_N = meshgridx(pinch(s_ext,2,length(mesh_step)),mesh_step);
  delta_s_ext_N = s_ext_N-s_ext;
  delta_s_ext_N = sqrt( integrate( delta_s_ext_N^2, 1, c/f) / integrate( s_ext_N^2, 1, c/f) );

  # read in the data from the file anc calculate delta_s_theory
  M = readdata("nanowire_au_jc_theory_from_mcm_fit.txt");
  s_theory = pinch(M,2,2)*1e-9;
  s_theory = meshgridx(s_theory,mesh_step);
  delta_s_theory= s_theory-s_ext;
  delta_s_theory = sqrt( integrate( delta_s_theory^2, 1, c/f) / integrate( s_theory^2, 1, c/f) );

  # plot delta_s, delta_s_N and delta_s_theory
  plotxy(dx(2:length(mesh_step))*1e9,delta_s_ext,
        dx*1e9,delta_s_ext_N,
        dx*1e9,delta_s_theory,
        "dx (nm)","delta_s");
  legend("delta_s","delta_s_N","delta_s_theory");

  # plot the best result at smallest dx vs theory
  plot(c/f*1e9,pinch(s_ext,2,length(dx))*1e9,pinch(s_theory,2,length(dx))*1e9,
       "wavelength (nm)","extinction cross section (nm)");
  legend("FDTD","theory");

}