clear; closeall;

f = [c/500e-9; c/700e-9];
n = [1.5; 2.13; 1.87; 1.94; 1.75; 0.644+5.28i];
d = [0; 0.09; 0.177; 0.048; 0.029; 0]*1e-6;
theta = 0;
field = stackfield(n,d,f,theta);
Es = pinch(field.Es,6,2); #extracting the y-component of Es. This value might be different for other polarizations
plot(field.z*1e6,abs(pinch(pinch(Es,4,1))),abs(pinch(pinch(Es,4,2))),"z (um)","E","Electric field within stack");
holdon;

switchtolayout;
run;

E = sqrt(getelectric("profile"));
y = getresult("profile","y");
plot(y*1e6,pinch(pinch(E,4,1)),pinch(pinch(E,4,2)));
legend("500 nm-analytical","700 nm-analytical","700 nm-FDTD","500 nm-FDTD");
holdoff;