dipole source spectrum problem

[chenyudong92]

Hi,

I just want to get transmission/reflection spectrum for a dipole source with structure. I don't know why the sum of transmission spectrum over 1. below is my script.

#%%
###
#module import
###
from __future__ import division
import meep as mp
import math
import numpy as np
import matplotlib.pyplot as plt
import cmath        #complex number math
###

####
#material design
####
resolution=40
n_SiN=1.984      #Si3N4 material index
n_material=n_SiN
n_SiO2=1.452    #SiO2 material index
n_air=1         #air material index
n_Si=3.646        #Si material index
ng=n_air        #the index of the source
dpml = 2      #the length of the PML boundary layer
dpad = 2      #the empty region
dpadx = 0     #the empty region
h_wg=0.34        #total height
h_grating=0.2      #the grating height
gp=0.566      #length of the grating period
gdc=0.5     #grating duty cycle (percentage)
h_sub=1.9     #height of substrate
h_sub1=3

mon_ps=0    #monitor position

wvl_min = 0.8          # min wavelength
wvl_max = 0.9           # max wavelength

theta=0     #the source angle (CCW: clock wise)

L_tooth=gp*gdc      #length of the tooth region
L_groove=gp*(1-gdc) #length of the groove region

L_total=5     #total length of the grating

N=int(L_total/gp)    #number of grating
PL_groove=[]
grating_offset_position=0 #move grating start position
###
#source control parameter
###
sx=(L_total+dpml+dpadx)*2
D_SMF=125   #single mode fiber diameter (um)
sigma = 5.2       # beam radius (beam width/2)
source_width=sx
# structure_offset=-D_SMF*np.sin(theta/180*np.pi)/2
# offset_y=D_SMF*np.sin(theta/180*np.pi)+h_wg+structure_offset
#structure_offset=-11/2
#offset_y=11+structure_offset
structure_offset=0
offset_y=0
#offset_y=2
####
#boundary setting
####
pml_layers = [mp.PML(thickness=dpml)]
sx=(L_total+dpml+dpadx)*2
sy=(dpml+dpad)*2+h_sub+h_wg+h_sub1*2+offset_y
cell_size=mp.Vector3(sx,sy,0)
###
#source range setting
###
nfreq = 100      #the number of the interpolation points between fmax and fmin
                #related with the monitor setting
#wvl_min = 1.45          # min wavelength
#wvl_max = 1.65           # max wavelength
fmin = 1/wvl_max        # min frequency
fmax = 1/wvl_min        # max frequency
fcen = 0.5*(fmin+fmax)  # center frequency
df = fmax-fmin          # frequency width

###
#define the function of the rotation source
###
tilt_angle = math.radians(theta) # angle of tilted beam
                              #counter clock wise(CCW) is the postive angel

k = mp.Vector3(y=1).rotate(mp.Vector3(z=1),tilt_angle).scale(fcen*ng)
#mp.Vector3().rotate(mp.Vector3(),tilt_angle), consisder as below.
#r(x,y,z) rotate [r1(x,y,z),angle]
#whic means,vector r rotate angle with vector r1
#.scale(fcen)
#the rotation vector is unit vector, we need get the k value of the center frequency
###
#Gaussian source setting
###
#sigma = 6         # beam width
#source_width=sx
#offset_y=D_SMF*np.sin(theta/180*np.pi)
#the source center, if we want to change angle, the position should near the PML to absourb another side
src_pt=mp.Vector3(mon_ps-0.5,h_wg+offset_y+2,0)
src_pt1=mp.Vector3(mon_ps+0.5,h_wg+offset_y+2,0)     #the source center, if we want to change angle, the position should near the PML to absourb another side
#%%

'''
def gaussian_beam(sigma, k, x0):
    def _gaussian_beam(x):
        return cmath.exp(1j*2*math.pi*k.dot(x-x0)-(x-x0).dot(x-x0)/(2*sigma**2)) #a.dot(b)=np.dot(a,b), inner produce
    return _gaussian_beam

sources=[mp.Source(#mp.GaussianSource(fcen, fwidth=df),    #continous source and Gassian source is the same
                    mp.ContinuousSource(fcen),
                     component=mp.Ez,   #for TE mode
                     center=src_pt,
                     amplitude=1,
                     size=mp.Vector3(sx,0,0),     #source width
                     amp_func=gaussian_beam(sigma,k,src_pt1)
                     )]
'''

#beam_x0 = mp.Vector3(0,3.0)    # beam focus (relative to source center)

rot_angle = theta  # CCW rotation angle about z axis (0: +y axis)

beam_kdir = mp.Vector3(0,-1,0).rotate(mp.Vector3(0,0,1),math.radians(rot_angle))  # beam propagation direction

beam_w0 = sigma   # beam waist radius

beam_E0 = mp.Vector3(0,0,1) 
'''
sources = [mp.GaussianBeamSource(#src=mp.ContinuousSource(fcen),
                                src=mp.GaussianSource(fcen, fwidth=df),

                                 center=src_pt,

                                 size=mp.Vector3(sx,0,0),   #source open window size

                                 beam_x0=mp.Vector3(0,0,0),  #beam focus: distance from source port to beam wasit

                                #component=mp.Ez

                                 beam_kdir=beam_kdir,

                                 beam_w0=beam_w0,

                                 beam_E0=beam_E0)]
'''
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),  #source type is Gassusian, the frequecy center is fcen, the band width is df.
                        #mp.ContinuousSource(fcen),
                                    amplitude=-1,

                                    component=mp.Ez, #the generation at the Hz component 

                                    size=mp.Vector3(0,0,0),

                                    center=src_pt),
            mp.Source(mp.GaussianSource(fcen, fwidth=df),  #source type is Gassusian, the frequecy center is fcen, the band width is df.
                        #mp.ContinuousSource(fcen),
                                    amplitude=1,

                                    component=mp.Ez, #the generation at the Hz component 

                                    size=mp.Vector3(0,0,0),

                                    center=src_pt1),]
#simulation seting
sim = mp.Simulation(cell_size=cell_size, #the whole running cell

                    geometry=[], #the structure detail set, it is a 1 by N matrix (list), which be desgined above.

                    sources=sources,# source set, there is blank

                    k_point=k,

                    boundary_layers=pml_layers,# boundary setting

                    default_material=mp.Medium(index=ng),

                    resolution=resolution) #resolution setting

###
#input monitor setting
###
mon_pt_in = mp.Vector3(0,-(sy/2-dpml),0)   #monitor position
flux_mon_in = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_in, size=mp.Vector3(x=sx-2*dpml)))    #set the monitor position and size

mon_pt_in_1 = mp.Vector3(0,(sy/2-dpml),0)   #monitor position
flux_mon_in_1 = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_in_1, size=mp.Vector3(x=sx-2*dpml)))    #set the monitor position and size

mon_pt_in_2 = mp.Vector3(-(sx/2-dpml),0,0)   #monitor position
flux_mon_in_2 = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_in_2, size=mp.Vector3(y=sy-2*dpml)))    #set the monitor position and size

mon_pt_in_3 = mp.Vector3((sx/2-dpml),0,0)   #monitor position
flux_mon_in_3 = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_in_3, size=mp.Vector3(y=sy-2*dpml)))    #set the monitor position and size

f = plt.figure(dpi=120)
sim.plot2D(ax=f.gca())
plt.show()
#%%
#%%

sim.run(until=30)
# plot the intenisty
#%%
non_pml_vol = mp.Volume(center=mp.Vector3(), size=cell_size)

f = lambda x: np.abs(x) ** 2

field_parameters = {'post_process':f,'cmap':'RdBu'}

sim.plot2D(fields=mp.Ez,field_parameters=field_parameters)

#plt.savefig('gaussian_beam.png')   #save picture
plt.show()

ez_data = sim.get_array(vol=non_pml_vol, component=mp.Ez)

plt.figure()
plt.imshow(np.flipud(np.transpose(np.real(ez_data))), interpolation='spline36', cmap='RdBu')
plt.axis('off')
plt.show()
#%%
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt_in, 1e-3))   #running the simulation unitl the fields decay to 1e-9 at position mon_pt;
# %%
sim.display_fluxes(flux_mon_in)
# %%
# %%
sim.display_fluxes(flux_mon_in)
# %%
flx_in=sim.display_fluxes(flux_mon_in)    #showing the flux
flx_in_1=sim.display_fluxes(flux_mon_in_1)
flx_in_2=sim.display_fluxes(flux_mon_in_2)
flx_in_3=sim.display_fluxes(flux_mon_in_3)
input_flux = mp.get_fluxes(flux_mon_in)    #save the flux field
input_flux_1 = mp.get_fluxes(flux_mon_in_1)    #save the flux field
input_flux_2 = mp.get_fluxes(flux_mon_in_2)
input_flux_3 = mp.get_fluxes(flux_mon_in_3)
#%%

######
#structure region
######
####
#psoition setting
####
for i in range(N):
    PL_groove=np.append(PL_groove,L_groove/2+gp*i)

print(PL_groove)

#%%
###
blk_wg=mp.Block(    material=mp.Medium(index=n_material),      #the material setting
                    size=mp.Vector3(mp.inf,h_wg,mp.inf),       #the block size, x=dmpl+dsub, y=inf, z=inf
                    center=mp.Vector3(0,h_wg/2+structure_offset,0)
                )

blk_sub=mp.Block(    material=mp.Medium(index=n_SiO2),      #the material setting
                    size=mp.Vector3(mp.inf,h_sub,mp.inf),       #the block size, x=dmpl+dsub, y=inf, z=inf
                    center=mp.Vector3(0,-h_sub/2+structure_offset,0)
                )

blk_sub1=mp.Block(  material=mp.Medium(index=n_Si),      #the material setting
                    size=mp.Vector3(mp.inf,h_sub1,mp.inf),       #the block size, x=dmpl+dsub, y=inf, z=inf
                    center=mp.Vector3(0,-h_sub1/2-h_sub+structure_offset,0)
                )
###
blk_grating=[]
for i in range (N):
    blk_grating=np.append(blk_grating,
                mp.Block(material=mp.Medium(index=n_air),      #the material setting
                size=mp.Vector3(L_groove,h_grating,mp.inf),       #the block size, x=dmpl+dsub, y=inf, z=inf
                center=mp.Vector3(PL_groove[i]+grating_offset_position,(h_grating/2+(h_wg-h_grating))+structure_offset,0)
                ))

#geometry_without_sub=[blk_wg]       #structure wihtout substrate
#geometry_without_sub.extend(blk_grating)
geometry_with_sub=[blk_wg,blk_sub,blk_sub1]  #structure with substrate
#geometry_with_sub.extend(blk_grating)
sim.reset_meep()
sim = mp.Simulation(cell_size=cell_size, #the whole running cell

                    geometry=geometry_with_sub, #the structure detail set, it is a 1 by N matrix (list), which be desgined above.

                    sources=sources,# source set, there is blank

                    k_point=k,

                    boundary_layers=pml_layers,# boundary setting

                    resolution=resolution) #resolution setting

###
#monitor setting    left trans
###
mon_pt = mp.Vector3(-(sx/2-dpml),0)   #monitor position
#nfreq = 100      #the number of the interpolation points between fmax and fmin
flux_mon = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy-2*dpml)))    #set the monitor position and size
###

###
#monitor setting    right trans
###
mon_pt_right = mp.Vector3((sx/2-dpml),0)   #monitor position
#nfreq = 100      #the number of the interpolation points between fmax and fmin
flux_mon_right = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_right, size=mp.Vector3(y=sy-2*dpml)))    #set the monitor position and size
###

#monitor setting bottom trans
###
mon_pt_bottom = mp.Vector3(0,-(sy/2-dpml),0)   #monitor position
#nfreq = 100      #the number of the interpolation points between fmax and fmin
flux_mon_bottom = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_bottom, size=mp.Vector3(x=sx-2*dpml)))    #set the monitor position and size

#monitor setting top trans
###
mon_pt_top = mp.Vector3(0,sy/2-dpml,0)   #monitor position
#nfreq = 100      #the number of the interpolation points between fmax and fmin
flux_mon_top = sim.add_flux(fcen, df, nfreq,                                        #Fourier Transform frequncy and function nunber setting
                        mp.FluxRegion(center=mon_pt_top, size=mp.Vector3(x=sx-2*dpml)))    #set the monitor position and size

f = plt.figure(dpi=120)
sim.plot2D(ax=f.gca())
plt.show()
#%%

###
#filed display
###
non_pml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(sx,sy/2,0))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-3))   #running the simulation unitl the fields decay to 1e-9 at position mon_pt;
ez_data = sim.get_array(vol=non_pml_vol, component=mp.Ez)

###
#left monitor calculation
###
sim.display_fluxes(flux_mon)    #showing the flux
trans_flux = mp.get_fluxes(flux_mon)    #save the flux field
flux_freqs = mp.get_flux_freqs(flux_mon)
###

###
#right monitor calculation
###
sim.display_fluxes(flux_mon_right)    #showing the flux
trans_flux_right = mp.get_fluxes(flux_mon_right)    #save the flux field
flux_freqs_right = mp.get_flux_freqs(flux_mon_right)
###

###
#bottom monitor calculation
###
sim.display_fluxes(flux_mon_bottom)    #showing the flux
trans_flux_bottom = mp.get_fluxes(flux_mon_bottom)    #save the flux field
flux_freqs_bottom = mp.get_flux_freqs(flux_mon_bottom)

###
#bottom monitor calculation
###
sim.display_fluxes(flux_mon_top)    #showing the flux
trans_flux_top = mp.get_fluxes(flux_mon_top)    #save the flux field
flux_freqs_top = mp.get_flux_freqs(flux_mon_top)

###
plt.figure(dpi=120)
plt.imshow(np.flipud(np.transpose(np.real(ez_data))), interpolation='spline36', cmap='RdBu')
#np.transpose:
#np.flipud:
plt.axis('on')
plt.show()

#%%

tran=[]
tran_bottom=[]
tran_top=[]
tran_right=[]
wl=[]
#transmission
input_flux_fix=[]
for i in range(len(flux_freqs)):
    input_flux_fix=np.append(input_flux_fix,input_flux[i]-input_flux_1[i]+input_flux_2[i]-input_flux_3[i])

for i in range(len(flux_freqs)):
    tran=np.append(tran,trans_flux[i]/(input_flux_fix[i]))
    tran_bottom=np.append(tran_bottom,trans_flux_bottom[i]/input_flux_fix[i])
    tran_top=np.append(tran_top,-trans_flux_top[i]/input_flux_fix[i])
    #tran_right=np.append(tran_right,1-(-trans_flux_top[i]+trans_flux_bottom[i]+trans_flux[i])/input_flux_fix[i])
    tran_right=np.append(tran_right,-trans_flux_right[i]/input_flux_fix[i])
    wl=np.append(wl,1/flux_freqs[i])

if mp.am_master():

    plt.figure();

    #plt.plot(wl,Rs,'bo-',label='reflectance')#plot refelectance with blue

    plt.plot(wl,tran, color='orange',label='Left')#plot transmittane with red
    plt.plot(wl,tran_bottom,'b-',label='Tran')#plot transmittane with red
    plt.plot(wl,tran_top,'r-',label='Refl')#plot transmittane with red
    plt.plot(wl,tran_right,'g-',label='Right')#plot transmittane with red
    plt.plot(wl,tran_bottom+tran_top+tran_right+tran,'g--',label='total')#plot loss with green

    #plt.axis([5.0, 10.0, 0, 1])#set the axis range

    plt.xlabel("wavelength (μm)")#set the x-axis unit label

    plt.ylabel("Transimissivity")

    plt.legend(loc="upper right")#set the lengend position

    plt.grid()

    plt.show()#plot the set before
#%% 

[stevengj]

You are computing the power P radiated by a dipole in one structure, and comparing it to the power P₀ radiated by the same dipole in vacuum, and are surprised that P/P₀ > 1?

It's perfectly physically possible to have P > P₀, as explained in the Meep FAQ. (In the case of a resonant structure, this is known as "Purcell enhancement.")

[chenyudong92]

Hi, Johnson

Thanks for your answer. But, When I try to use dipole source (two Gaussian source with opposite amplitude) to get the transmission/reflection spectrum with a single layer Si without other structures. The spectrum will change when I change the simulation region. I am not sure how to get power P₀ radiated by the same dipole in vacuum to give a stable spectrum. [cid:3daf9a9e-0ec1-4dcc-846a-7ec6b59f67b6] [cid:cfb1dd78-e48b-4e38-9f35-b2346485404f] [cid:35047ec6-6323-4451-ad8d-9e2177bd589d]

Sincerely, Yudong Chen M.S. of Electrical Engineering at SMU Ph.D. Student of Electrical Engineering at UTA Office at 2nd floor at Nano Fab in UTA

---

发件人: Steven G. Johnson notifications@github.com 发送时间: 2020年12月10日 14:53 收件人: NanoComp/meep meep@noreply.github.com 抄送: Chen, Yudong yudong.chen@mavs.uta.edu; Author author@noreply.github.com 主题: Re: [NanoComp/meep] dipole source spectrum problem (#1448)

You are computing the power P radiated by a dipole in one structure, and comparing it to the power P₀ radiated by the same dipole in vacuum, and are surprised that P/P₀ > 1?

It's perfectly physically possible to have P > P₀, as explained in the Meep FAQhttps://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmeep.readthedocs.io%2Fen%2Flatest%2FFAQ%2F%23how-does-the-current-amplitude-relate-to-the-resulting-field-amplitude&data=04%7C01%7Cyudong.chen%40mavs.uta.edu%7C7bd2deb2fa014cd9201708d89d45381b%7C5cdc5b43d7be4caa8173729e3b0a62d9%7C1%7C0%7C637432267903980777%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=95DXaKBwS2nb%2BUkwS6fAMAQr8s8SUFqVWUsSSQnNKMs%3D&reserved=0. (In the case of a resonant structure, this is known as "Purcell enhancement.")

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FNanoComp%2Fmeep%2Fissues%2F1448%23issuecomment-742758039&data=04%7C01%7Cyudong.chen%40mavs.uta.edu%7C7bd2deb2fa014cd9201708d89d45381b%7C5cdc5b43d7be4caa8173729e3b0a62d9%7C1%7C0%7C637432267903990769%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=G9ApCIXLyqjWcMjzmVM4ZBR0635Pz7%2FgUuDSu%2B4BAFM%3D&reserved=0, or unsubscribehttps://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FANXYTVYDYYN6EAULTJRF32DSUERKLANCNFSM4UVPXLZQ&data=04%7C01%7Cyudong.chen%40mavs.uta.edu%7C7bd2deb2fa014cd9201708d89d45381b%7C5cdc5b43d7be4caa8173729e3b0a62d9%7C1%7C0%7C637432267904000766%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=KwYclzu365P51Cm8XvKJn9aVgmvCdxrdSD%2F3O2ELrro%3D&reserved=0.

[stevengj]

You probably need to re-think the physics of what you are doing. The concept of a transmission or reflection spectrum is usually for a planewave source, and it's not clear what it means for a dipole source (because the presence of the structure changes the power emitted by a dipole).