# load("nanoarray.ldev");
monitorR = "::model::DGTD::R";
monitorT = "::model::DGTD::T";

fluxR = getresult(monitorR,"flux");
fluxT = getresult(monitorT,"flux");

R = 0.5*(fluxR.T_front + fluxR.T_back);
T = 0.5*(fluxT.T_front + fluxT.T_back);
lambda = fluxR.lambda;
assert("Monitors R, T should have matching lambda", (lambda(1) == fluxT.lambda(1))*(lambda(end) == fluxT.lambda(end)));
if(lambda(1) > lambda(end)) {
    lambda = flip(lambda,1);
    R = flip(R,1);
    T = flip(T,1);
}

##
## Plot: T and R
##

plot(lambda*1e9, R, T, R+T);
legend("R","T","R+T");
setplot("x label","wavelength(nm)");
setplot("y label","T and R");
setplot("x min",1e9*min(lambda));
setplot("x max",1e9*max(lambda));
setplot("y min",0);
setplot("y max",1);

##
## Plot: T and R (matching FDTD)
##

plot(lambda*1e9, R, T, R+T);
legend("R","T","R+T");
setplot("x label","wavelength(nm)");
setplot("y label","T and R");
setplot("x min",400);
setplot("x max",800); # to match FDTD example
setplot("y min",0);
setplot("y max",1);

##
## Plot: T normalized (matching FDTD)
##
scale = pi*(getnamed("::model::geometry::hole","radius"))^2/(getnamed("::model::simulation region","x span")*getnamed("::model::simulation region","y span"));
plot(lambda*1e9, T/scale);
legend("");
setplot("title","T normalized to pi*r^2/unit_cell_area");
setplot("x label","wavelength(nm)");
setplot("y label","T (normalized)");
setplot("x min",400);
setplot("x max",800); # to match FDTD example
setplot("y min",0);
setplot("y max",2.5); # to match FDTD example

##
## Near Fields
##
monitorName = "::model::DGTD::profile_at_700nm";
profileFields = getresult(monitorName,"fields");
deltaZ = 2.5e-9;

##
## Plot |E|^2, z=-0.06um
##
xC = linspace(-0.2e-6,0.2e-6,203);
yC = linspace(-0.2e-6,0.2e-6,201);
zC = -0.06e-6 + linspace(-deltaZ,deltaZ,13);
fallbackValue = 0.;

data3D = interptet(profileFields.elements, [profileFields.x, profileFields.y, profileFields.z], profileFields.E2, xC, yC, zC, fallbackValue);
data2D = sum(data3D,3)/size(data3D,3);
image(xC*1e6,yC*1e6,data2D,"x(um)","y(um)","|E|^2 at reflected surface");
setplot("colorbar min",5);
setplot("colorbar max",260);

##
## Plot |E|^2, z=+0.06um
##
xC = linspace(-0.2e-6,0.2e-6,203);
yC = linspace(-0.2e-6,0.2e-6,201);
zC = 0.06e-6 + linspace(-deltaZ,deltaZ,13);
fallbackValue = 0.;

data3D = interptet(profileFields.elements, [profileFields.x, profileFields.y, profileFields.z], profileFields.E2, xC, yC, zC, fallbackValue);
data2D = sum(data3D,3)/size(data3D,3);
image(xC*1e6,yC*1e6,data2D,"x(um)","y(um)","|E|^2 at transmitted surface");
setplot("colorbar min",0);
setplot("colorbar max",27);

##
## Plot |E|^2, y=0um
##
xC = linspace(-0.2e-6,0.2e-6,303);
yC = [0.];
zC = linspace(-0.3e-6,0.3e-6,451);
fallbackValue = 0.;

data2D = interptet(profileFields.elements, [profileFields.x, profileFields.y, profileFields.z], profileFields.E2, xC, yC, zC, fallbackValue);
image(xC*1e6,zC*1e6,data2D,"x(um)","z(um)","|E|^2 at y=0");
setplot("colorbar min",0);
setplot("colorbar max",10);