#####################################################
# surface_plasmon.lsf
#
# This script calculates the effective index and 
# propagation loss for a surface plasmon mode
# as well as the analytic result and calculates
# the error as a function of grid accuracy. It also
# compares the results when a mesh override region is
# used at the interface between the two materials.
#
# Copyright 2016 Lumerical Solutions Inc.
#####################################################

# load simulation file
load("surface_plasmon.lms");

# load initial TM mode data to track during sweep
loaddata("surface_plasmon.ldf");

#############################################################
# calculate the analytic effective index and propagation loss
#############################################################
lam = 632.8e-9;
ko = 2*pi/lam;
ei = 1;
nm = 0.23783 + 1i*3.384755;
em = nm^2;

ksp2 = ko^2*ei*em/(ei+em);
ksp = sqrt(ksp2);
Nsp = ksp/ko;
N = real(Nsp); # analytic effective index
L = 40*pi/lam/log(10)*imag(Nsp)/1000; # analytic propagation loss

##################################################################################
# sweep through number of mesh cells along the y direction from 100 cells doubling
# for each point up to 12800 mesh cells
##################################################################################
mesh_cell = matrix(1,8); # initialize vector of number of mesh cells to use in each simulation
mesh_cell(1) = 100; # first value
for(i=2:length(mesh_cell)){ # double each subsequent number of mesh cells
   mesh_cell(i)=mesh_cell(i-1)*2;
}

# initialize variables to hold the simulated effective index and propagation loss results
N1 = matrix(1,length(mesh_cell));
L1 = matrix(1,length(mesh_cell));

for(ii = 1:length(mesh_cell)){
    switchtolayout;
    setnamed("MODE","mesh cells y",mesh_cell(ii));
    findmodes;
    N1(ii) = real(getdata(bestoverlap("TM1"),"neff"));
    L1(ii) = getdata(bestoverlap("TM1"),"loss")/1000;
}
switchtolayout;

########################################################
# add a mesh override region and perform the sweep again
########################################################
addmesh;
set("override y mesh",1);
set("set mesh multiplier",1);
set("y mesh multiplier",20);
set("x",0);
set("x span",2e-6);
set("y min",0);
set("y max",0.03e-6);

# since a mesh override region is also used, the actual number of mesh cells
# is larger than the base number of mesh cells we set in the FDE solver object
# initialize a variable to store the actual number of mesh cells used in each
# simulation with the mesh override 
mesh_cell2 = matrix(1,8);

# initialize variables to store the effective index and propagation loss
N2 = matrix(1,length(mesh_cell));
L2 = matrix(1,length(mesh_cell));

for(ii = 1:length(mesh_cell)){
    switchtolayout;
    setnamed("MODE","mesh cells y",mesh_cell(ii));
    findmodes;
    N2(ii) = real(getdata(bestoverlap("TM1"),"neff"));
    L2(ii) = getdata(bestoverlap("TM1"),"loss")/1000;
    mesh_cell2(ii)=getnamed("MODE","actual mesh cells y");
}

#############################
# calculate error amplitude %
#############################
# error for uniform mesh (without mesh override)
p1 = 100*abs((N1-N)/N);
p2 = 100*abs((L1-L)/L);

# error for graded mesh (when mesh override is added)
p3 = 100*abs((N2-N)/N);
p4 = 100*abs((L2-L)/L);

########################################
# plot results against analytic solution
########################################
# create vectors containing analytic values at all points to plot against simulated data
N=matrix(1,length(mesh_cell2))+N;
L=matrix(1,length(mesh_cell2))+L;

plot(mesh_cell2,N2,N,"grid points","effective index","graded mesh","");
legend("calculated","theory");

plot(mesh_cell2,L2,L,"grid points","propagation loss(dB/mm)","graded mesh","");
legend("calculated","theory");

plot(mesh_cell2,p3,p4,"grid points","error amplitude(%)","graded mesh","logplot");
legend("effective index","propagation loss");

# compare uniform mesh and graded mesh error amplitudes for effective index
plotxy(mesh_cell,p1,mesh_cell2,p3,"grid points","error amplitude(%)","effective index error amplitude","logplot");
legend("uniform mesh","graded mesh");