# usr_bloch_angle_sweep.lsf
# Do a parameter sweep over wavelength and angle of incidence
# using bloch boundary conditions.  The simulation results
# must be interpolated onto a rectangular grid. 

clear; closeall;


########################################################
# User setup parameters

# set range of angles
ntheta = 90;                   
theta = linspace(10,45,ntheta);

# set frequency range
f_max = c/0.5e-6;
f_min = c/1e-6;

# set number of simulations to run.
nk = 30;  # number of kInPlane points to sweep

# set additional configuration settings
index_source = getresult("source","index");  # refractive index at source location
n = index_source.index; # refractive index at source location

fname = "sweep_";		       # file name 			
run_sims = 1;                            # 1-> run sims, 0-> analyze pre-run sims



########################################################
# Calculate the kInPlane values to be used in the simulation.  
# First, determine largest and smallest kInPlane values (these will occur at the limits of Fmin/max and THETAmin/max
# Next, create a linearly spaced vector from K min to K max
# Note that due to the sign conventions used in FDTD Solutions, postive angles (theta) correspond to 
# negative kInPlane values.

kTemp = -2*pi*n/c * [f_min*sin(min(theta)*pi/180), f_min*sin(max(theta)*pi/180), f_max*sin(min(theta)*pi/180), f_max*sin(max(theta)*pi/180)];
kInPlane = linspace(min(kTemp),max(kTemp),nk);


########################################################
# Next calculate the appropriate min/max frequency for the source.  
#
# The calculation of fStop (the higher frequency) is simple: just use the max frequency of interest for all simulations.
#
# The calculation of fStart (the lower frequency) is more complicated.  
# For smaller abs(kInPlane) values, where the effective source angle at f min is smaller than the maximum theta value 
# of interest, we use fStart=fmin
# For larger abs(kInPlane) values, where the effective source angle at f min is larger than the maximum theta value 
# of interest, we calculate a different fStart value to ensure the source angle is not larger than 90 degrees.  
# The filter is a 1 or 0 vector to help prevent fStart result a source angle larger than 90 degree.
#
# Also, IF statements are required to deal with cases where theta min/max are positive or negative

fStop  = linspace(f_max,f_max,nk);
 
if (max(theta) > 0) { 
  if (min(theta) >= 0) {  # case where theta min & max are both positive
    filter = kInPlane >  (-2*pi*n/c*f_min*sin(max(theta)*pi/180));
    fStart = linspace(f_min,f_min,nk)                   *  filter +
             -kInPlane*c/2/pi/n/sin(max(theta)*pi/180)  * !filter ;
  } else {                # case where theta min is negative and theta max is positive
    filter1 = kInPlane >  (-2*pi*n/c*f_min*sin(max(theta)*pi/180));
    filter2 = kInPlane <  (-2*pi*n/c*f_min*sin(min(theta)*pi/180));
    fStart = linspace(f_min,f_min,nk)                   * (filter1 & filter2) +
             -kInPlane*c/2/pi/n/sin(max(theta)*pi/180)  * !filter1 +
             -kInPlane*c/2/pi/n/sin(min(theta)*pi/180)  * !filter2 ;
  }
} else {                  # case where theta min & max are both negative
    filter = kInPlane <  (-2*pi*n/c*f_min*sin(min(theta)*pi/180));
    fStart = linspace(f_min,f_min,nk)                   *  filter + 
             -kInPlane*c/2/pi/n/sin(min(theta)*pi/180)  * !filter ;
}



########################################################
# Once the fStart, fStop freequency is calculated for each simulation in the sweep, we calculate
# the corresponding source angle for each simulation.
fCenter = (fStart+fStop)/2;
thetaSIM = asin(-kInPlane*c/2/pi/n/fCenter)*180/pi;



########################################################
# Run the simulations
if (run_sims) {

  # override the simulation bandwidth, so it is not inhereted from
  # the source.  This ensures all other simulation parameters
  # (such as the mesh size, material fits, monitor frequency vector, etc)
  # don't change between individual simulations.
  # We add 10% to the range to avoid potential interpolation problems
  switchtolayout;
  setnamed("FDTD","set simulation bandwidth",1);
  setnamed("FDTD","simulation frequency min",f_min);  setnamed("FDTD","simulation frequency max",f_max);

  # set source frequency range and angle
  for (i=1:nk) {
    switchtolayout;
    setnamed("source","frequency start", fStart(i));
    setnamed("source","frequency stop",  fStop(i));
    setnamed("source","angle",           thetaSIM(i));

    save(fname+num2str(i));
    addjob(fname+num2str(i));
  }
  runjobs;  # run all simulations
} else {
  ?"Analyzing pre-run simulations";
}


########################################################
# Get original data from simulations
load(fname+"1"); 
mname = "T"; 

# initialize matrices
f  = getdata(mname,"f");      # monitor frequency vector
nf = length(f);
T_raw =       matrix(nk,nf);        # raw transmission data from each simulation     
sourceAngle = matrix(nk,nf);        # corresponding source angle data
T_final =     matrix(ntheta,nf);    # final (interpolated) transmission data

# get data from each simulation
for (i=1:nk) {
  load(fname+num2str(i)); 
  T_raw(i,1:nf) = transmission(mname);  
  sourceAngle(i,1:nf) = getsourceangle("source",f); 
}



########################################################
# interpolate to standard source angle vector
for (i=1:nf) {
  T_final(1:ntheta,i) = interp(pinch(T_raw,2,i),pinch(sourceAngle,2,i),theta);
}



########################################################
# plot some results
image(f/1e12,theta,transpose(T_final),"frequency (THz)","theta (deg)","Transmission");   ## LUM CHANGE
setplot("colorbar min",0);  setplot("colorbar max",1);      ## LUM CHANGE

plot(f/1e12,transpose(real(sourceAngle)),
     "f (THz)","theta (deg)","Source angle data - raw");
setplot("y min",-90); setplot("y max",90);

plot([f_min, f_min, f_max, f_max,f_min]/1e12, 
     [max(theta),min(theta),min(theta),max(theta),max(theta)],
     "f (THz)","theta (deg)","Source angle data - desired");
setplot("y min",-90); setplot("y max",90);

plot(f/1e12,meshgridy(f,theta),
     "f (THz)","theta (deg)","Source angle data - final");
setplot("y min",-90); setplot("y max",90);

plot(f/1e12,meshgridy(f,kInPlane),
     "f (THz)","kInPlane","kInPlane - raw");
plotxy(fStart/1e12,kInPlane,fStop/1e12,kInPlane,
       "f (THz)","kInPlane","kInPlane - desired");
legend("start frequency","stop frequency");