Emission from a point source

Sorry for the late reply. The emission feature is not really developed, I did not pursue it further. The mathematical details and further explanation can be found in supplementary chapter 5 of my paper: https://opg.optica.org/oe/fulltext.cfm?uri=oe-29-22-36201&id=461916#articleSupplMat

It is a similar story with the gyrating media in the other issue. I provided the foundation, but did not look more into polishing or extensive testing.

The three functions in emission.jl are meant to be used sequentially. The first one converts the 3-vector of currents to the fourier-space description (Dirac distribution around lattice point 0/0).

For the second one, you need the combined scattering matrices of all the layers before and after the point source, respectively. This then calculates how much light is emitted in the form of the amplitude vectors a0 and b0 which describe the light propagating forward and backward from the source. You can then use the scattering matrices to find out where the radiation goes.

The third one computes the power that was radiated from current and electric field, as it is dependent on the dielectric surroundings -> Purcell effect.

Note that what is simulated here is not a single emitter, but a coherent array with each unit cell containing one emitter.

Below you can find an old example file I used to test this. But it has been a while and I currently dont have much time to look into it. It computes radiated powers from emitters in a Si metasurface on SOI.

If you are interested and motivated to implement this, I would be happy if you could add the feature to this repo, or collaborate on it.

using RigorousCoupledWaveAnalysis,LinearAlgebra,ProgressMeter
radmat=InterpolPerm(RigorousCoupledWaveAnalysis.si_schinke)
geo=Circle(.5)
ox=ModelPerm(sio2_malitson)
si=InterpolPerm(si_schinke)
zno=SimpleLayer(50,ModelPerm(RigorousCoupledWaveAnalysis.sio2_malitson))
tsi=PatternedLayer(50,[ox,si],[geo])
act1=PatternedLayer(75,[ox,radmat],[geo])
act2=PatternedLayer(75,[ox,radmat],[geo])
bsi=SimpleLayer(80,InterpolPerm(RigorousCoupledWaveAnalysis.si_schinke))
box=SimpleLayer(1000,ModelPerm(RigorousCoupledWaveAnalysis.sio2_malitson))
mdl=rcwamodel([zno,tsi,act1,act2,bsi,box],ConstantPerm(1.0),InterpolPerm(RigorousCoupledWaveAnalysis.si_schinke))
wls=1000:20:1600
pitch=479*3
N=5
rad=zeros(length(wls))
up=zeros(length(wls))
dn=zeros(length(wls))
eps=zeros(length(wls))*1im
@showprogress for i=1:length(wls)
    λ=wls[i]
    grd=rcwagrid(N,N,pitch,ppitch1E-5,90,λ,ConstantPerm(1.0))
    s=scatMatrices(mdl,grd,λ)
    sup=halfspace(grd.Kx,grd.Ky,mdl.εsup,λ)
    sub=halfspace(grd.Kx,grd.Ky,mdl.εsub,λ)
    Shi=concatenate(s[1:4])
    Slo=concatenate(s[5:end])
    px,py,pz=RigorousCoupledWaveAnalysis.pointDipole(grd,0,1,0)
    a0,b0=RigorousCoupledWaveAnalysis.innerSource(grd,px,py,pz,Shi,Slo)
    rad[i]=real.(RigorousCoupledWaveAnalysis.dipoleRad(a0,b0,Slo,grd,sub,px,py,pz))/real.(sqrt(RigorousCoupledWaveAnalysis.get_permittivity(radmat,λ)))
    up[i]=.5RigorousCoupledWaveAnalysis.a2p(Shi.S12*b0,I,grd.Kx,grd.Ky,sup.Kz,grd.kin[3])/real.(sqrt(RigorousCoupledWaveAnalysis.get_permittivity(radmat,λ)))
    dn[i]=.5RigorousCoupledWaveAnalysis.a2p(Slo.S21*a0,I,grd.Kx,grd.Ky,sub.Kz,grd.kin[3])/real.(sqrt(RigorousCoupledWaveAnalysis.get_permittivity(radmat,λ)))
    #ste,stm=etmSource(grd.kin,N,N)
    #R,T=etm_reftra(ste,mdl,grd,λ)
    #result[i]=λ/1240*(1 .-R-T)
    eps[i]=RigorousCoupledWaveAnalysis.get_permittivity(radmat,λ)
end
plot(wls,rad,label="radiated from source")
plot!(wls,up,label="to superstrate")
plot!(wls,dn,label="to substrate")
plot!(wls,rad-up-dn,label="absorbed")

jonschlipf / RigorousCoupledWaveAnalysis.jl

Emission from a point source #3