##############################################################
# scriptfile: run_simulations.lsf
#
# Description:  This script creates a sweep over source settings
#               for plane wave illumination within the illumination
#               NA
#
# Copyright 2016 Lumerical Solutions Inc.
#########################################################

#########################################################
# COMMON SETTINGS
# These must be the same in the run and build S-matrix files

NA_max0 = 0.9; # define max NA for input beam

run_filename = "metal_sim"; # name for temporary files, no extension
rundir = "rundir"; # name of subdirectory to store temporary files

#########################################################

# get the supercell periodicity from the FDTD region
select("FDTD");
Lx = get("x span");
Ly = get("y span");

# get the frequency from the source
select("source1");
f = get("center frequency");
k0 = 2*pi*f/c;

# choose the in plane k vector (typically 0,0)
kx0_center = 0;
ky0_center = 0;

# construct the corresponding grating orders
n = [ ceil( -(k0+kx0_center)*Lx/(2*pi) ) : floor( (k0+kx0_center)*Lx/(2*pi) )];
m = [ ceil( -(k0+ky0_center)*Ly/(2*pi) ) : floor( (k0+ky0_center)*Ly/(2*pi) )];

# construct the in plane component of the unit vectors corresponding to 
# the grating orders
ux0 = (n*2*pi/Lx + kx0_center)/k0;
uy0 = (m*2*pi/Ly + ky0_center)/k0;

# get the current filename and directory
original_filename = currentfilename;
current_dir = pwd;
base_filename = pwd + "/" + rundir +"/" +run_filename;

# turn off auto-redrawing
redrawoff;

# clear any jobs on queue
clearjobs;

# loop over all incoming angles and polarizations
simcount = 0;
for(i=1:length(ux0)) {
  for(j=1:length(uy0)) {
    for(p=0:90:90) {
      uz0 = sqrt( 1-ux0(i)^2-uy0(j)^2 );
      uxy0 = sqrt( ux0(i)^2 + uy0(j)^2 );
      theta = acos(uz0)*180/pi;
      phi = atan2(uy0(j),ux0(i))*180/pi;

      # ONLY include incoming angles that land within the desired incoming NA
      if(real(uxy0) <= NA_max0*1.001) {
        simcount = simcount + 1;
        switchtolayout;
        select("source1");
        set("angle theta",theta);
        set("angle phi",phi);
        set("polarization angle",p);
        redraw;
        filename = base_filename+"_"+num2str(i)+"_"+num2str(j)+"_"+num2str(p);
        save(filename);
 
        # add jobs to the queue
        addjob(filename);
      } 
    }
  }
}

# resave the original file
save(original_filename);

# indicate the total number of jobs to run
?"Queued " + num2str(simcount) + " simulations";

# run the jobs
runjobs;