clear; closeall;

f = c/500e-9;
n = [1.5; 2.13; 1.87; 1.94; 1.75; 0.644+5.28i];
d = [0; 0.09; 0.177; 0.048; 0.029; 0]*1e-6;
n_dipoles = 301; # number of dipoles
# dipoles evenly placed through active layer
z = linspace(sum(d(1:3))+5e-10,sum(d(1:4))-5e-10,n_dipoles); 
orientation = "random";

result = stackpurcell(n,d,f,z,orientation);

# plot results
density = result.density;
theta = density.getparameter("theta");
downward_emission = pinch(density.getattribute("downward"));
image(z*1e6,theta,transpose(downward_emission),"z (microns)","theta","far field power density");

# separate quantum efficiency and extraction efficiency
sin_theta = matrix(1000,n_dipoles);
for (i=1:n_dipoles){sin_theta(1:1000,i) = sin(theta*pi/180);}

F = pinch(result.power.getattribute('purcell_factor'));
Frad = integrate(downward_emission*sin_theta,1,theta*pi/180)*2*pi;
plot(z*1e6,Frad/F,F,Frad,'z (microns)','','','loglog');
legend('Extraction Efficiency','Quantum Efficiency',"Far Field Angular Density");