###############################################################################
# OLED_total_QE_analysis_3D_square.lsf
#
# This script file cacalates the emission enhancement, spectrum, angular distribution
# by projecting the data to the far field for analysis.
# Users will need to run the parameter sweep before analysing the results
###############################################################################

clear; # clear existing script variables to prevent confusion

# run the parameter sweep if there are no sweep results
if (havesweepresult!=1){
 runsweep; save;
}

#########################################################
# User inputs
datafile = "OLED_farfield_results.ldf"; # file to store far field results
project_in_air = 1; # 1 for air, 0 for glass
plot_all_wavelengths = 0; # 0 by default, 1 for generating plots for all wavelengths
res = 201; # set this value to match the resolution set in the analysis variables
           # of the far_field_change_index analysis group, reduce this no. to speed up

# color filter parameters
lambda0 = [460; 510; 620 ]*1e-9;
FWHM = 15e-9;
#########################################################

# load sweep parameters and results
px = getsweepdata('sweep','px');
py = getsweepdata('sweep','py');
weight = getsweepdata('sweep','weight');
dipole_orientation = getsweepdata('sweep','dipole_orientation');
simpoints = length(px); # total number of points in the sweep
sym_45 = sum(dipole_orientation==1)>sum(dipole_orientation==2); # check to see if rotation symmetry is used
E2_far = getsweepresult('sweep','E2_far');
E2_far_air = E2_far.E2_far2;
E2_far_glass = E2_far.E2_far1;
ux = E2_far.ux;
uy = E2_far.uy;
f = E2_far.f;

# initialize result matrices
result_x = matrix(res,res,length(f));
result_y = matrix(res,res,length(f));
result_z = matrix(res,res,length(f));

# possibility of reducing the number of frequency points to study
p1 = 1;         # index of first frequency point in range of interest and
p2 = length(f); # index of last frequency point in range of interest

# loop over all simulation locations for the patterned case
# this corresponds to the first 7 points in the sweep
for(i=1:simpoints-2) {
  ?"Analyzing point " + num2str(i) + " of " + num2str(simpoints);
  temp = matrix(res,res);
  for(fpoint=p1:p2) { # loop over each frequency point
    if(project_in_air) { # project in air
      if(dipole_orientation(i)==1){
         # x oriented dipoles
         temp = pinch(pinch(E2_far_air,4,i),3,fpoint);
         result_x(1:res,1:res,fpoint) = pinch(result_x,3,fpoint) + weight(i)*temp;
         if ((sym_45) and (px(i)==py(i))){ 
	 result_y(1:res,1:res,fpoint) = pinch(result_y,3,fpoint) + weight(i)*transpose(temp);	  
         }
       } else {
          if(dipole_orientation(i)==2){
          # y oriented dipoles
          result_y(1:res,1:res,fpoint) = pinch(result_y,3,fpoint) + weight(i)*pinch(pinch(E2_far_air,4,i),3,fpoint);
          } else {
	 # z oriented dipoles
	 result_z(1:res,1:res,fpoint) = pinch(result_z,3,fpoint) + weight(i)*pinch(pinch(E2_far_air,4,i),3,fpoint);
          }
       }
    } else { # project in glass
       if(dipole_orientation(i)==1){
         # x oriented dipoles
         temp = pinch(pinch(E2_far_glass,4,i),3,fpoint);
         result_x(1:res,1:res,fpoint) = pinch(result_x,3,fpoint) + weight(i)*temp;
         if ((sym_45) and (px(i)==py(i))){
	result_y(1:res,1:res,fpoint) = pinch(result_y,3,fpoint) + weight(i)*transpose(pinch(pinch(E2_far_glass,4,i),3,fpoint));	  
         }
       } else {
         if (dipole_orientation(i)==2){
           # y oriented dipoles
	result_y(1:res,1:res,fpoint) = pinch(result_y,3,fpoint) + weight(i)*pinch(pinch(E2_far_glass,4,i),3,fpoint);
         } else {
	 # z oriented dipoles
	 result_z(1:res,1:res,fpoint) = pinch(result_z,3,fpoint) + weight(i)*pinch(pinch(E2_far_glass,4,i),3,fpoint);
         }
       }
    }
  } # end loop over frequency points  
} # end loop over sweep points in patterned case

if (sym_45){
# symmetry was used, unfold to 1/4 of unit cell
  for(fpoint=p1:p2) { # loop over frequency points
  result_x(1:res,1:res,fpoint) = pinch(result_x,3,fpoint) + transpose(pinch(result_x,3,fpoint));
  result_y(1:res,1:res,fpoint) = pinch(result_y,3,fpoint) + transpose(pinch(result_y,3,fpoint));
  result_z(1:res,1:res,fpoint) = pinch(result_z,3,fpoint) + transpose(pinch(result_z,3,fpoint));
  }
}

# unfold to full cell
result_x = result_x + 
         result_x(res:-1:1,1:res,1:length(f)) +
         result_x(1:res,res:-1:1,1:length(f)) +
         result_x(res:-1:1,res:-1:1,1:length(f));
result_y = result_y + 
         result_y(res:-1:1,1:res,1:length(f)) +
         result_y(1:res,res:-1:1,1:length(f)) +
         result_y(res:-1:1,res:-1:1,1:length(f));
result_z = result_z + 
         result_z(res:-1:1,1:res,1:length(f)) +
         result_z(1:res,res:-1:1,1:length(f)) +
         result_z(res:-1:1,res:-1:1,1:length(f));
result = (1/3)*(result_x+result_y+result_z);

# calculate results when no PC is present
?"Calculating results without PC";
result0_x = matrix(res,res,length(f));
result0_y = matrix(res,res,length(f));
result0_z = matrix(res,res,length(f));


for(i=simpoints-1:simpoints) {
  ?"Analyzing point " + num2str(i) + " of " + num2str(simpoints);
  for(fpoint=p1:p2) { # loop over frequency points
    if(project_in_air) { # project in air
      temp = pinch(pinch(E2_far_air,4,i),3,fpoint);
      if(dipole_orientation(i) == 1){
        result0_x(1:res,1:res,fpoint) = pinch(result0_x,3,fpoint) + temp; 
        result0_y(1:res,1:res,fpoint) = pinch(result0_y,3,fpoint) + transpose(temp);
      } else{
      result0_z(1:res,1:res,fpoint) = pinch(result0_z,3,fpoint) + temp;
      }
    } else { # project in glass
      temp = pinch(pinch(E2_far_glass,4,i),3,fpoint);
      if(dipole_orientation(i) == 1){
        result0_x(1:res,1:res,fpoint) = pinch(result0_x,3,fpoint) + temp; 
        result0_y(1:res,1:res,fpoint) = pinch(result0_y,3,fpoint) + transpose(temp); 
      } else{
        result0_z(1:res,1:res,fpoint) = pinch(result0_z,3,fpoint) + temp;
      }
    }
  } # end loop over frequency points
} # end loop over sweep points for non-patterened case

result0 = (1/3)*(result0_x+result0_y+result0_z);

#########################################################
#plot the far field profile at all frequencies
if(plot_all_wavelengths) {
  for(fpoint=1:length(f)) {
    temp = num2str(round(c/f(fpoint)*1e9))+"nm";
    image(ux,uy,pinch(result,3,fpoint),"","",temp,"polar");
    pause(0.1); # for better viewing
    exportfigure("OLED_"+temp);
  }
}

# calculate extraction efficiency at each wavelength
power_radiated = 0.5*sqrt(eps0/mu0)*farfield3dintegrate(result,ux,uy);
power_radiated0 = 0.5*sqrt(eps0/mu0)*farfield3dintegrate(result0,ux,uy);

# calculate and plot the enhancement 
enhancement = power_radiated/power_radiated0;
plot(c/f*1e9,enhancement,"wavelength (nm)","Enhancement");

# apply color filters
filter = matrix(length(f),length(lambda0));
for(i=1:length(lambda0)) {
  temp = exp( -(c/f-lambda0(i))^2/(FWHM/2/sqrt(log(2)))^2);
  temp = temp/integrate(temp,1,f);
  filter(1:length(f),i) = temp;
}
plot(c/f*1e9,filter,"wavelength (nm)","filters");

filtered_result = matrix(length(ux),length(uy),length(lambda0));
df = f(2)-f(1);
for(i=1:length(lambda0)) {
  for(fpoint=1:length(f)) {
    filtered_result(1:length(ux),1:length(uy),i) = pinch(filtered_result,3,i) 
                           + filter(fpoint,i)*df*pinch(result,3,fpoint);
  }
  temp = num2str(round(lambda0(i)*1e9))+"nm";
  image(ux,uy,pinch(filtered_result,3,i),"","","lambda0 = "+temp,"polar");
}

# save the far field results to file
savedata(datafile,result,result_x,result_y,result_z,result0,result0_x,result0_y,result0_z,ux,uy,f);