##############################################################
# scriptfile: build_S_matrix.lsf
#
# Description:  This script collects the data form the simulations
#               to build the S-matrix and store it in a ldf file
#
# Copyright 2016 Lumerical Solutions Inc.
#########################################################

# clear existing workspace
clear;

#########################################################
# 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;
S_filename = pwd + "/" +run_filename + "_S_matrix";

# turn off auto-redrawing
redrawoff;

# monitor name
mname = "analysis";

# prevent redraws
redrawoff;

# callocate storage ES and EP
# arguments are ux_in, uy_in, polarization_in (S,P), ux_out, uy_out, polarization_out (S,P)
S = matrix(length(ux0),length(uy0),2,length(ux0),length(uy0),2);

# loop over all input angles
for(i=1:length(ux0)) {
for(j=1:length(uy0)) {

  # output progress
  ?num2str(i)+", " + num2str(j);

  # loop over input polarizations
    for(p=0:90:90) {

      # find filename for this simulation
      filename = base_filename+"_"+num2str(i)+"_"+num2str(j)+"_"+num2str(p);
      uz0 = sqrt( 1-ux0(i)^2-uy0(j)^2 );
      uxy0 = sqrt( ux0(i)^2 + uy0(j)^2 );
      #determine angles in degrees
      theta = acos(uz0)*180/pi;
      phi = atan2(-uy0(j),-ux0(i))*180/pi;

      # only include angles within the NA of the source
      if(real(uxy0) <= NA_max0*1.001) {
        load(filename);
        if(!havedata(mname)) {
          runanalysis(mname);
          save(filename);
        }
        
        # frequency and k vectors of the input beams
        f = getdata(mname,"f");
        k0 = 2*pi*f/c;
        kx0 = -k0*ux0(i);
        ky0 = -k0*uy0(j);
        kz0 = -k0*uz0; 
        R = 1+transmission("analysis::fields"); # monitor in front of the source

        # get the Ep and Es of the light
        ux = getdata(mname,"ux");
        uy = getdata(mname,"uy");
        zs = getdata("source1","z");

        if(p==0) { # p-polarized
          pol_in = 2;
        } else {
          pol_in = 1;
        }

        # note: extra - sign from source, only for positive source propagation
        # 1i from diffraction, 1m of phase to 0, correct for source position. Now have phase at global 0
        # this is a plane wave decomposition of the resulting fields, normalized to 
        # the global position (0,0,0)
        S(i,j,pol_in,1:length(ux),1:length(uy),1) = 
              1i*exp(-1i*k0-1i*kz0*(0-zs))*pinch(getdata(mname,"ES"))*sqrt(R);
        S(i,j,pol_in,1:length(ux),1:length(uy),2) = 
              1i*exp(-1i*k0-1i*kz0*(0-zs))*pinch(getdata(mname,"EP"))*sqrt(R);
      }

    }
}
}

# reload original file or switch back to original directory
if(original_filename != "") {
  load(original_filename);
} else {
  cd(current_dir);
}

# save S parameters
savedata(S_filename,S,ux0,uy0,f);