###############################################################################
# OLED_internal_QE_analysis_3D_square.lsf
#
# This script file cacalates the radiated power based on the parameter sweep results.
# It calculates the radiated power recorded by the trans_box and dipolepower command for
# with and without PC patterning. The results can be unwrapped if lattice symmetry techniques
# are applied. Users will need to run the parameter sweep before analysing the resutls
###############################################################################

clear; # clear existing script variables to prevent confusion
save("OLED_3D_square.fsp");
datafile = "OLED_dipolepower_results.ldf"; # file to store dipole power results

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

#########################################################
# load sweep parameters and results
px = getsweepdata('sweep','px');
py = getsweepdata('sweep','py');
weight = getsweepdata('sweep','weight');
dipole_orientation = getsweepdata('sweep','dipole_orientation');
dipolePower = getsweepresult('sweep','dipolePower');
source_power = dipolePower.source_power;
dipole_power = dipolePower.dipole_power;
dipole_box_power = dipolePower.dipole_power_box;
simpoints = length(px);
sym_45 = sum(dipole_orientation==1)>sum(dipole_orientation==2); # check to see if rotation symmetry is used

select("source1");
nf = length(source_power)/simpoints;
f = linspace(get("frequency start"),get("frequency stop"),nf);
sp = matrix(length(f));   # sourcepower
sp0 = matrix(length(f));  # sourcepower, no PC
dp1 = matrix(length(f));  # dipolepower command
dp2 = matrix(length(f));  # trans_box
dp01 = matrix(length(f)); # dipolepower command, no PC
dp02 = matrix(length(f)); # trans_box, no PC

# loop over all simulation locations
for(i=1:simpoints-2) { # last 2 sims are no PC patterning
  sp = sp + weight(i)*source_power(1:nf,i);
  dp1 = dp1 + weight(i)*dipole_power(1:nf,i);
  dp2 = dp2 + weight(i)*dipole_box_power(1:nf,i);
  # if symmetry is used and dipole is along diagonal
  if ((px(i)==py(i))and (sym_45)){?[px(i),py(i)];
    # need to duplicate results for y orientation if sym is used 	
    if (dipole_orientation(i)==1){
      sp = sp + weight(i)*source_power(1:nf,i);
      dp1 = dp1 + weight(i)*dipole_power(1:nf,i);
      dp2 = dp2 + weight(i)*dipole_box_power(1:nf,i);
    }
  }
}
if (sym_45){
# if 45 degree rotational symmetry was used, unfold to 1/4 of unit cell
# a factor of 2 for 45 deg rot sym, a factor of 4 for a quarter of unit cell, 
# a factor of 1/3 for three dipole orientations
sp = sp*8/3; dp1 = dp1*8/3; dp2 = dp2*8/3;
} else {
# do not need to unfold to 1/4 cell
sp = sp*4/3; dp1 = dp1*4/3; dp2 = dp2*4/3;
}

# calculate results when no PC is present (last two simpoints)
sp0 = sp0 + 2*source_power(1:nf,simpoints-1)+source_power(1:nf,simpoints);
dp01 = dp01 + 2*dipole_power(1:nf,simpoints-1)+dipole_power(1:nf,simpoints);
dp02 = dp02 + 2*dipole_box_power(1:nf,simpoints-1)+dipole_box_power(1:nf,simpoints);
sp0 = (1/3)*sp0; dp01 = (1/3)*dp01; dp02 = (1/3)*dp02;

# calculate the power radiated by the dipoles with and without
# patterning, compared to the homogeneous case
plot(c/f*1e9,dp1,dp2,dp01,dp02,"Radiated power");
legend("PC dipole","PC box","no PC dipole","no PC box");

# save the dipole power results to file
savedata(datafile,dp1,dp2,dp01,dp02,f);