# scriptfile: dipole_gf.lsf
############################################################################################
### Green's function  ######################################################################
############################################################################################
runanalysis;
f=getdata("trans_box","f");

component = "Ez0"; # should match the dipole orientation!
eps_r=1;

# get the field recorded automatically by the source
Ez = pinch(getdata("source1",component));
t = getdata("source1","time0");
Ez_w = czt(Ez,t,2*pi*f)*(t(2)-t(1));

# get the field from the source for Greens function Calculation
moment = getnamed("source1","total amplitude");
m = moment*getdata("source1","time_signal");
ts = getdata("source1","time");
m_w = czt(m,ts,2*pi*f)*(ts(2)-ts(1));
w = 2*pi*f;
G=(Ez_w)*c^2*eps0/(w^2*m_w);

#plot(f,real(G), "f", "real(G)", "real(G)");  This is not meaningful since in theory it approaches infinity
plot(f/1e12,imag(G),"f (THz)", "imag(G)", "imag(G)");

############################################################################################
### Radiated power  ########################################################################
############################################################################################

# calculate the power from the local Green's function, normalized
# to the power in a homogeneous medium

P=imag(G)*w^3*abs(m_w)^2/2/c^2/eps0/abs(sourcenorm(f))^2;
p_G = P / sourcepower(f);

# calculate the power using the built-in dipolepower script function,
# normalized to the power in a homogeneous region
p_D = dipolepower(f) / sourcepower(f);

# calculate power radiated as measured by the box of monitors
p_B = getdata("trans_box","T");

plot(c/f*1e9,p_B,p_D,p_G,"wavelength (nm)","normalized radiated power","radiated power compared to homogenous medium");
legend("box of monitors","dipolepower function", "from local Green's function");

############################################################################################
### Partial local density of states ########################################################
############################################################################################

# calculate rho from Green's function
rho_z= [6*w/pi/c^2]*imag(G); 

# calculate rho directly from dipolepower
rho_z_d= 12*eps0*dipolepower(f)*abs(sourcenorm(f))^2/pi/w^2/abs(m_w)^2;

plot(f/1e12,rho_z,rho_z_d, "f (THz)", "rho","partial local density of states");
legend("from Green's function","from dipole power");

############################################################################################
### spontenous decay rate ##################################################################
############################################################################################

gamma= w*pi/3/hbar/eps0*abs(m_w)^2*rho_z;
plot(f/1e12,gamma, "f (THz)", "gamma", "spontenous decay rate");

############################################################################################
### Total local density of states ########################################################
############################################################################################

# Three different orientations of the dipole can be examined and the results can be incoherently added:

#rho_total=rho_x+rho_y+rho_z;