redrawoff;
slash = "/";
if (operatingsystem == "windows") {slash = "\\";}
#########################################
# Set constants
mkm = 1e-6;
groupscope("::model");
select("FDTD");
min_mesh_step=get("min mesh step");
select("sub-Au");
sub_Au_z_max = get("z max");
#########################################
# Create results dir if needed
curdir = pwd;
dirname_base = filebasename(currentfilename)+"."+fileextension(currentfilename)+".results";
dirname = curdir+slash+dirname_base;

files = splitstring(dir,endl);
isCreated = false;
for(i=1:length(files)) {
    if (files{i} == dirname_base) { isCreated = true;}
}
if (isCreated){ ?"output dir already exist";
} else { system("mkdir "+dirname);}
#########################################
# Get FDTD material data
groupscope("::model::mon");
min_f = getglobalmonitor("minimum frequency");
max_f = getglobalmonitor("maximum frequency");
fpoints = getglobalmonitor("frequency points");
fvec=linspace(min_f,max_f,fpoints);
#Au_name = "Au (Gold) - Johnson and Christy Copy 2";
Au_name = "Au-McPeak";
n_fdtd=getfdtdindex(Au_name,fvec,min_f,max_f);
n_fdtd;
#########################################
?"Save data";
cd(dirname);
?dirname;
format long;
z_arr = [5, 20];
for(x=1:10) {
    for (media=1:2){
        if (media == 1){shift_sign = 1;} else {shift_sign = -1;}
        for(i=1:length(z_arr)){
            zshift_all = z_arr(i);
            monitor_name = "mon_x"+num2str(x)+"mkm_media"+num2str(media)+"_zshift"+num2str(zshift_all)+"nm";
            select(monitor_name);
            if (haveresult(monitor_name) == 1){
                E = getresult(monitor_name, "E");
                H = getresult(monitor_name, "H");
                len = length(E.lambda);
                #Save to file
                fname = dirname+slash+monitor_name+".txt";
                if (fileexists(fname)){rm(fname);}
                write(fname,"lambda, dip.power, Ex, Ey, Ez, Hx, Hy, Hz, n_Au");
                for (i=1:len){
                    str = num2str(E.lambda(i))
                    +", "+num2str(dipolepower(E.f(i)))
                    +", "+num2str(E.Ex(i))
                    +", "+num2str(E.Ey(i))
                    +", "+num2str(E.Ez(i))
                    +", "+num2str(H.Hx(i))
                    +", "+num2str(H.Hy(i))
                    +", "+num2str(H.Hz(i))
                    +", "+num2str(n_fdtd(i))
                    ;
                    write(fname,str);
                }
            }
        }
    }
}
cd(curdir);
redrawon;