###########################################################
# sers_pt_ag.lsf
#
# This script calculates and plots the SERS enhancement
# factor for XZ and YZ cross sections of a Pt nanosphere
# above a silver substrate
#
# Copyright 2016 Lumerical Solutions Inc.
###########################################################

clear;
closeall;

#get the electric intensity
exz=pinch(getelectric("XZ"));
eyz=pinch(getelectric("YZ"));

#find the maximum
mxz=max((exz)^2);
myz=max((eyz)^2);
# get some raw data
x1=getdata("XZ","x");
z1=getdata("XZ","z");
y2=getdata("YZ","y");
z2=getdata("YZ","z");
f=getdata("YZ","f");
#
EFxz=exz^2;#enhancement factor
index = find(EFxz,mxz);
# convert index to row, col indices
matrix_size = size(EFxz);
indices = matrix(length(matrix_size));
# do for each dimension
for (i = 1:length(matrix_size)) {
 mod_dividend = index;
 mod_divisor  = matrix_size(i);
 mod_remainder = mod(mod_dividend,mod_divisor);
 if (mod_remainder == 0) { mod_remainder = matrix_size(i); }
 indices(i) = mod_remainder;
 # remove this dimension from further calculations
 index = (index+(matrix_size(i)-mod_remainder))/matrix_size(i);
}
?"max at x1 ="+num2str(x1(indices(1))*1e6)+" um";
?"multi indice access: EF factor xz ("+num2str(indices(1))+","+
                          num2str(indices(2))+","+
                          num2str(indices(3))+")="+
                          num2str(EFxz(indices(1),indices(2),indices(3)));
# users may need to modify those values in order to have proper view of the resulting images
zmin=-5e-9;
zmax= 20e-9;
xym = 20e-9;
#
nx1=find(x1,-xym)-1;
nx2=find(x1, xym)+1;
nz1=find(z1,zmin)-1;
nz2=find(z1,zmax)+1;
image(x1(nx1:nx2)*1e9,z1(nz1:nz2)*1e9,EFxz(nx1:nx2,nz1:nz2,indices(3)),"x nm","z nm","EF xz");
image(x1(nx1:nx2)*1e9,z1(nz1:nz2)*1e9,exz(nx1:nx2,nz1:nz2,indices(3)),"z nm","z nm","Intensity xz");

# in another cross section
EFyz=eyz^2;
index = find(EFyz,myz);
# convert index to row, col indices
matrix_size = size(EFyz);
indices = matrix(length(matrix_size));
# do for each dimension
for (i = 1:length(matrix_size)) {
 mod_dividend = index;
 mod_divisor  = matrix_size(i);
 mod_remainder = mod(mod_dividend,mod_divisor);
 if (mod_remainder == 0) { mod_remainder = matrix_size(i); }
 indices(i) = mod_remainder;
 # remove this dimension from further calculations
 index = (index+(matrix_size(i)-mod_remainder))/matrix_size(i);
}
?"multi indice access: EF factor yz ("+num2str(indices(1))+","+
                          num2str(indices(2))+","+
                          num2str(indices(3))+")="+
                          num2str(EFyz(indices(1),indices(2),indices(3)));
?"max at y1 ="+num2str(y2(indices(1))*1e6)+" um";
ny1=find(y2,-xym)-1;
ny2=find(y2, xym)+1;
nz1=find(z2,zmin)-1;
nz2=find(z2,zmax)+1;
image(y2(ny1:ny2)*1e9,z2(nz1:nz2)*1e9,EFyz(ny1:ny2,nz1:nz2,indices(3)),   "y nm","z nm","EF yz");
image(y2(ny1:ny2)*1e9,z2(nz1:nz2)*1e9,eyz(ny1:ny2,nz1:nz2,indices(3)),"y nm","z nm","Intensity yz");