#########################################
# Problem setup
#

# geometry: material n1 film of specified thickness
n1 = 1.5; # lower halfspace material
n2 = 1.0; # upper halfspace material (emission medium)

# source: monochrome 500nm dipole
# position: in material n2 delta nm above interface
wavelength = 500e-9;
delta = 80e-9;

currentDensityInAmpPerM2 = 1.0;
excitonFraction = 1.0;
singletRatio = 0.25;
relativeDecayRate = 1.0;

# angular_res: resolution for emission angle (farfield angle)
angular_res = 173;
theta = linspace(0,90,num=angular_res+1);
theta = theta(1:angular_res);
sin_theta = sin(pi/180*theta);
cos_theta = cos(pi/180*theta);

#########################################
# Nominal Power: (rate of photon production)*(energy per photon)
#
electronsPerNm2PerSec = currentDensityInAmpPerM2 * 6.2415e18;   # eg: 1 Coulomb/s
#planckTimesCInJouleNm = 198.6446582 * 1e-18; 
probabilityChargeCarrierResultsInPhoton = excitonFraction * singletRatio;
normalizedPowerInWattPerM2 = electronsPerNm2PerSec * probabilityChargeCarrierResultsInPhoton * h * c / wavelength;

#########################################
# Stackdipole
#
n = [n1; n2];
d = [0, 2*delta];
f = [c/wavelength];
z = [delta];
spectrum = [1.0];

data_Vert = stackdipole(n,d,f,z,spectrum,"vert",angular_res);
data_Horz = stackdipole(n,d,f,z,spectrum,"horz",angular_res);

P_Horz_sdipl = data_Horz.radiance;
P_Vert_sdipl = data_Vert.radiance;

data_uniform = stackdipole(n,d,f,z,spectrum,"rand",angular_res);
P_uniform_sdipl = data_uniform.radiance;

#########################################
# FDTD
#

# filename prefix to use
name_prefix = "fdtd_dipole_halfspace";

# choose to rerun the simulation or not
rerun_simulations = 0; # set to 0 to not rerun simulations

load(name_prefix); # load the initial file
if(rerun_simulations) {

    # dipole orientation 1 (along x-axis)
    switchtolayout;
    setnamed("source1","theta",90);
    setnamed("source1","phi",0);
    save(name_prefix+"_0");
    run;

    # dipole orientation 2 (along y-axis)
    switchtolayout;
    setnamed("source1","theta",90);
    setnamed("source1","phi",90);
    save(name_prefix+"_90");
    run;

    # dipole orientation 3 (along z-axis)
    switchtolayout;
    setnamed("source1","theta",0);
    setnamed("source1","phi",0);
    save(name_prefix+"_z");
    run;
}

# turn on smoothing of the far field projection
farfieldfilter(1);

# collect data from the horizontal-oriented simulations
load(name_prefix+"_0");
monname = "above";
ux = farfieldux(monname);
uy = farfielduy(monname);
phi = 0;

E2_0 = farfield3d(monname);
T_0 = transmission(monname);
angDistrib_0 = E2_0 / (farfield3dintegrate(E2_0,ux,uy)); 
angDistrib_0 = farfieldspherical(angDistrib_0,ux,uy,theta,phi);
fdtd_pHorz = T_0 * angDistrib_0;
fdtd_pHorz_dipolepower = dipolepower(c/wavelength);
fdtd_pHorz_sourcepower = sourcepower(c/wavelength);

load(name_prefix+"_90");
ux = farfieldux(monname);
uy = farfielduy(monname);
phi = 0;

E2_90 = farfield3d(monname);
T_90 = transmission(monname);
angDistrib_90 = E2_90 / (farfield3dintegrate(E2_90,ux,uy)); 
angDistrib_90 = farfieldspherical(angDistrib_90,ux,uy,theta,phi);
fdtd_sHorz = T_90 * angDistrib_90;
fdtd_sHorz_dipolepower = dipolepower(c/wavelength);
fdtd_sHorz_sourcepower = sourcepower(c/wavelength);

# Since T_0, T_90 include equal portions of power "transmitted from each polarization"
fdtd_pHorz = fdtd_pHorz * 0.5; 
fdtd_sHorz = fdtd_sHorz * 0.5;

# collect data from the vertical-oriented simulation
load(name_prefix+"_z");
ux = farfieldux(monname);
uy = farfielduy(monname);
phi = 0;

T_z = transmission(monname);
E2_z = farfield3d(monname);
angDistrib_z = E2_z / (farfield3dintegrate(E2_z,ux,uy)); 
angDistrib_z = farfieldspherical(angDistrib_z,ux,uy,theta,phi);
fdtd_pVert = T_z * angDistrib_z;
fdtd_pVert_dipolepower = dipolepower(c/wavelength);
fdtd_pVert_sourcepower = sourcepower(c/wavelength);


# convert from simulation units to specified nominal power
P_pHorz_fdtd = fdtd_pHorz*relativeDecayRate/( relativeDecayRate*fdtd_pHorz_dipolepower/fdtd_pHorz_sourcepower + (1-relativeDecayRate) )*normalizedPowerInWattPerM2;
P_sHorz_fdtd = fdtd_sHorz*relativeDecayRate/( relativeDecayRate*fdtd_sHorz_dipolepower/fdtd_sHorz_sourcepower + (1-relativeDecayRate) )*normalizedPowerInWattPerM2;
P_Vert_fdtd = fdtd_pVert*relativeDecayRate/( relativeDecayRate*fdtd_pVert_dipolepower/fdtd_pVert_sourcepower + (1-relativeDecayRate) )*normalizedPowerInWattPerM2;
P_uniform_fdtd = (2.0*(fdtd_pHorz+fdtd_sHorz)+fdtd_pVert)/3*relativeDecayRate/( relativeDecayRate*(fdtd_pHorz_dipolepower+fdtd_sHorz_dipolepower+fdtd_pVert_dipolepower)/(fdtd_pHorz_sourcepower+fdtd_sHorz_sourcepower+fdtd_pVert_sourcepower) + (1-relativeDecayRate) )*normalizedPowerInWattPerM2;

# reload the template simulation
load(name_prefix);

#########################################
# Plot Comparisons
#
closeall;

res_dipole = P_Horz_sdipl;
res_fdtd = (P_pHorz_fdtd+P_sHorz_fdtd);
res_title = "Horizontal";

plot(theta,res_dipole,res_fdtd,"emission angle","power/steradian",res_title);
legend("stackdipole","fdtd");
?res_title;
?max(abs(res_dipole-res_fdtd));
?max(abs(res_dipole-res_fdtd))/max(abs(res_dipole));

res_dipole = P_Vert_sdipl;
res_fdtd = P_Vert_fdtd;
res_title = "Vertical";

plot(theta,res_dipole,res_fdtd,"emission angle","power/steradian",res_title);
legend("stackdipole","fdtd");
?res_title;
?max(abs(res_dipole-res_fdtd));
?max(abs(res_dipole-res_fdtd))/max(abs(res_dipole));

res_dipole = P_uniform_sdipl;
res_fdtd = P_uniform_fdtd;
res_title = "Uniform";

plot(theta,res_dipole,res_fdtd,"emission angle","power/steradian",res_title);
legend("stackdipole","fdtd");
?res_title;
?max(abs(res_dipole-res_fdtd));
?max(abs(res_dipole-res_fdtd))/max(abs(res_dipole));