In[]:=
epsd=5.9;​​c=3*10^17;​​h=6.582119*10^-16;​​γp=0.04;​​Δz=2;​​ωp=9;​​ϵ0=8.85*10^(-12)*(10^-9)*(1.6*10^-19);​​ϵ∞=4.4;​​epsm[ω_]:=ϵ∞-ωp^2/(ω*(ω+I*γp));​​R=30;(*AssumedvalueforR*)​​r=R+Δz;​​dz=0;​​d21=2.0016*10^-19;​​​​η[q_,ω_]:=
2
ω
2
(h*c)
*epsd-
2
q
;​​ξ[q_,ω_]:=
2
ω
2
(h*c)
*epsm[ω]-
2
q
;​​​​​​xi[n_,x_]:=HankelH1[n,x];​​dxi[n_,x_]:=
1
2
(HankelH1[-1+n,x]-HankelH1[1+n,x]);​​psi[n_,x_]:=BesselJ[n,x];​​dpsi[n_,x_]:=
1
2
(BesselJ[-1+n,x]-BesselJ[1+n,x]);​​
In[]:=
In[]:=
green[k_,l_,n_,q_,ω_]:=
1
epsd
2
r
2
ω
2
η[q,ω]
dz(k-l)q

epsd
2
n
2
ω
psi[n,rη[q,ω]]xi[n,rη[q,ω]]+epsd
2
n
2
ω
ae[n,q,ω]
2
xi[n,rη[q,ω]]
+chqrdxi[n,rη[q,ω]]η[q,ω]-2
epsd
nωbo[n,q,ω]xi[n,rη[q,ω]]+chqr(dpsi[n,rη[q,ω]]+ce[n,q,ω]dxi[n,rη[q,ω]])η[q,ω],0,
1
epsdr
2
ω
c
dz(k-l)q

hxi[n,rη[q,ω]]
epsd
nωbo[n,q,ω]xi[n,rη[q,ω]]+chqr(dpsi[n,rη[q,ω]]+ce[n,q,ω]dxi[n,rη[q,ω]])η[q,ω],0,
1
epsd
2
r
2
ω
2
η[q,ω]
dz(k-l)q


2
c
2
h
2
n
2
q
xi[n,rη[q,ω]](psi[n,rη[q,ω]]+ce[n,q,ω]xi[n,rη[q,ω]])-2c
epsd
hnqrωbo[n,q,ω]dxi[n,rη[q,ω]]xi[n,rη[q,ω]]η[q,ω]+epsd
2
r
2
ω
dxi[n,rη[q,ω]](dpsi[n,rη[q,ω]]+ae[n,q,ω]dxi[n,rη[q,ω]])
2
η[q,ω]
,0,-
1
epsdr
2
ω
c
dz(k-l)q

h-
epsd
nωbo[n,q,ω]
2
xi[n,rη[q,ω]]
+chqrdxi[n,rη[q,ω]](psi[n,rη[q,ω]]+ce[n,q,ω]xi[n,rη[q,ω]])η[q,ω],0,
2
c
dz(k-l)q

2
h
xi[n,rη[q,ω]](psi[n,rη[q,ω]]+ce[n,q,ω]xi[n,rη[q,ω]])
2
η[q,ω]
epsd
2
ω
;
In[]:=
v={1/Sqrt[2]{I,0,1},1/Sqrt[2]{-I,0,1}};​​​​dens[k_,l_,ω_]:=dens[k,l,ω]=(((d21)^2*ω^2)/(8π^2*ϵ0*(h*c)^2))*ParallelSum[NIntegrate[Re[(2-KroneckerDelta[0,n])*(Conjugate[v[[k]]].green[k,l,n,q,ω].v[[l]])],{q,-10,10},AccuracyGoal->6,PrecisionGoal->6,MaxRecursion->20,Method->"LocalAdaptive"],{n,0,30}];
In[]:=
denn=Table[Interpolation[Table[{i,dens[k,l,i]},{i,0.01,5,0.01}]],{k,1,2},{l,1,2}];​​denb[ω_]:=Table[denn[[i]][[j]][ω],{i,1Length[denn]},{j,1Length[denn]}];
In[]:=
Plot[denb[ω],{ω,0.2,4},PlotStyle->Green,PlotRange->All]​​Table[Plot[denb[ω][[i]][[j]],{ω,0.2,4},PlotRange->All],{i,1,2},{j,1,2}]​​
Out[]=
Out[]=

,
,
,

Out[]=