###########################################################################
# Scriptfile: negative_index_chen_analysis.lsf
#
# Description: This file plots the normalized transmission and reflection
#	as a function of frequency for the simulation negative_index_chen.fsp
#	which models the structure described in
#	H.-T. Chen et al., "Active terahertz metamaterial devices", 
#       Nature 444, 597-600 (2006). It also plots |E| .
#
#
# Copyright 2007, Lumerical Solutions, Inc.
###########################################################################

# get the frequency, R and T 
f = getdata("T","f");
T = transmission("T");
R = -transmission("R");

# plot the results for R and T
plot(f*1e-12,R,T,"frequency (THz)","Transmission");
legend("R","T");

# see if we can calculate complex r and t coefficients
Pz_R = getdata("R","Pz");
Pz_T = getdata("T","Pz");
?"max(|imag(Pz)/real(Pz)|) in transmission: " + num2str(max(abs(imag(Pz_T)/real(Pz_T))));
?"max(|imag(Pz)/real(Pz)|) in reflection: " + num2str(max(abs(imag(Pz_R)/real(Pz_R))));

# calculate r and t based on E fields at x,y = (0,0);
t = getdata("T","Ey");
x = getdata("T","x");
y = getdata("T","y");
t = t(find(x,0),find(y,0),1,1:length(f));
r = getdata("R","Ey");
x = getdata("R","x");
y = getdata("R","y");
r = r(find(x,0),find(y,0),1,1:length(f));

# compare |r|^2 and n2/n1*|t|^2 to R and T

# The next 2 lines are commented out since there are multiple grating orders in the 
# high index slab. Due to the multiple grating orders, the method of calculting r based 
# on E(0,0) does not work.

# plot(f*1e-12,R,abs(r)^2,"frequency (THz)","Transmission");
# legend("R","|r|^2");

n2 = 1;
n1 = 3.5;
plot(f*1e-12,T,n2/n1*abs(t)^2,"frequency (THz)","Transmission");
legend("T","n2/n1*|t|^2");

# plot the phase of r and t
plot(f*1e-12,unwrap(angle(r)),unwrap(angle(t)),"frequency (THz)","phase of r and t (radians)");
legend("angle(r)","angle(t)");

# find the positions of the resonances (lowest point in transmission)
res_pos = findpeaks(1/T,2);
?"Resonances are found to be " + num2str(f(res_pos(1))*1e-12) + " THz, and "+ num2str(f(res_pos(2))*1e-12) + " THz";

# plot |E|^2 at lower resonance frequency
res_pos_low = min(res_pos);
mname="profile_surface";
x = getdata(mname,"x");
y = getdata(mname,"y");
E2 = getelectric(mname);
image(x*1e6,y*1e6,sqrt(pinch(E2,4,res_pos_low)),"x (microns)","y (microns)","|E| at " + num2str(f(res_pos_low)*1e-12) + " THz");

# plot surface current at resonance frequency
fi = res_pos_low;

# first create a spatial filter to include only data over the surface of the conductor
nx = pinch(getdata("material_properties","index_x")); # get index over space from index monitor
nx = real(nx);  # get real refractive index
filter = nx > 3.5;

# Get field data from monitors 1 cell above and below the conductive sheet
m = "profile_above";
Hx_above = pinch(pinch(getdata(m,"Hx"),4,fi));
Hy_above = pinch(pinch(getdata(m,"Hy"),4,fi));
m = "profile_below";
Hx_below = pinch(pinch(getdata(m,"Hx"),4,fi));
Hy_below = pinch(pinch(getdata(m,"Hy"),4,fi));
x = getdata(m,"x");
y = getdata(m,"y");
z=0;

# Calculate surface current and multiply by filter to get the result only on the conductor surface
Ky = real(Hx_above-Hx_below)*filter;
Kx = -real(Hy_above-Hy_below)*filter;
Kz=0*Kx;

# plot H field profiles at resonance
image(x*1e6,y*1e6,Ky,"x (um)","y (um)","Ky = Hx (A/m)");
image(x*1e6,y*1e6,Kx,"x (um)","y (um)","Kx = -Hy (A/m)");

# plot vector plot of surface current
K = rectilineardataset("K",x,y,z);
K.addattribute("K",Kx,Ky,Kz);
vectorplot(K);
setplot("down sample",1);
setplot("Nx",50);
setplot("Ny",50);
setplot("enable scale for vector",1);
setplot("scale factor",10);
setplot("invert background color",1);
setplot("title",m);