ZeroDivisionError: float division by zero.

[niubiworker]

I was trying to use Meep to calculate the transmission loss of photonic crystal waveguides, but after running it, I encountered the following error: ZeroDivisionError: float division by zero. Why does this problem occur?

import meep as mp
import matplotlib.pyplot as plt
from meep.materials import LiNbO3
resolution = 10 # pixels/a

a = 0.775         # units of um
r = 0.315*a         # units of um
h = 0.3       # units of um

period_num = 50
wave_length = 1.55
fcen = 1/wave_length
df = 0.1
pml_height = 0.5*wave_length
y_num = 6

geometry_up_boundary = y_num*np.sqrt(3)*a+2*pml_height
geometry_left_boundary = period_num*a+2*pml_height
size = mp.Vector3(geometry_left_boundary,h+10*pml_height,geometry_up_boundary*2)
w = np.sqrt(3)*a

geometry = [mp.Block(center=mp.Vector3(0,0,0),size=mp.Vector3(geometry_left_boundary*2,h,geometry_up_boundary), material=LiNbO3)]

pml_layers = [mp.PML(pml_height)]

sources = [mp.Source(mp.GaussianSource(frequency=fcen,fwidth=df),
                     component=mp.Ez,
                     center=mp.Vector3(-geometry_left_boundary/2+pml_height,0,0),
                     size=mp.Vector3(0,h,w))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=size,
                    geometry=geometry,
                    boundary_layers=pml_layers,
                    sources=sources,
                    eps_averaging=True,
                    symmetries=[mp.Mirror(mp.Y,phase=1)],
                    Courant=0.02)

nfreq = 10
ref1_local = -period_num/2*a+pml_height+0.5*a

ref1_fr = mp.FluxRegion(center=mp.Vector3(ref1_local,0,0),
                        size=mp.Vector3(0,2*h,2*w),
                        direction=mp.X)
ref1 = sim.add_flux(fcen,df,nfreq,ref1_fr)

tran_fr = mp.FluxRegion(center=mp.Vector3(ref1_local+period_num*a,0,0),
                        size=mp.Vector3(0,2*h,2*w),
                        direction=mp.X)
tran = sim.add_flux(fcen,df,nfreq,tran_fr)

# run
pt = mp.Vector3(ref1_local+period_num*a-0.5*a,0,0)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3))
print(sim.get_flux_data(tran))
straight_ref1_data = sim.get_flux_data(ref1)
straight_tran_flux = mp.get_fluxes(tran)

sim.reset_meep()

geometry = [mp.Block(center=mp.Vector3(0,0,0),size=mp.Vector3(geometry_left_boundary*2,h,geometry_up_boundary), material=LiNbO3)]
for i in np.linspace(-geometry_left_boundary,geometry_left_boundary,period_num*2+1):
    for j in np.linspace(-geometry_up_boundary,geometry_up_boundary,2*y_num+1):
        if j != 0:
            geometry += [mp.Cylinder(center=mp.Vector3(i, 0, j, ), height=mp.inf, material=mp.air, radius=r,axis=mp.Vector3(y=1))]
            geometry += [mp.Cylinder(center=mp.Vector3(i + 0.5 * a, 0, np.sqrt(3) / 2 * a + j, ), height=mp.inf, material=mp.air,radius=r,axis=mp.Vector3(y=1))]
            geometry += [mp.Cylinder(center=mp.Vector3(i + 0.5 * a, 0, -np.sqrt(3) / 2 * a + j, ), height=mp.inf, material=mp.air,radius=r,axis=mp.Vector3(y=1))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=size,
                    geometry=geometry,
                    boundary_layers=pml_layers,
                    sources=sources,
                    symmetries=[mp.Mirror(mp.Y,phase=1)],
                    Courant=0.02)

ref1 = sim.add_flux(fcen,df,nfreq,ref1_fr)
tran = sim.add_flux(fcen,df,nfreq,tran_fr)
sim.load_minus_flux_data(ref1,straight_ref1_data)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3))

pc_ref1_flux = mp.get_fluxes(ref1)
pc_tran_flux = mp.get_fluxes(tran)

flux_freqs = mp.get_flux_freqs(ref1)

wl = []
Rs = []
Ts = []
for i in range(nfreq):
    wl = np.append(wl, 1 / flux_freqs[i])
    Rs = np.append(Rs, -pc_ref1_flux[i] / straight_tran_flux[i])
    Ts = np.append(Ts, pc_tran_flux[i] / straight_tran_flux[i])

plt.plot(wl, Rs, "bo-", label="reflectance")
plt.plot(wl, Ts, "ro-", label="transmittance")
plt.plot(wl, 1 - Rs - Ts, "go-", label="loss")
plt.xlim(min(wl),max(wl))
plt.xlabel("wavelength (μm)")
plt.legend(loc="upper right")
plt.show()```

[leSemaleon]

I had this problem too.
Rs = np.append(Rs, -pc_ref1_flux[i] / straight_tran_flux[i])
Ts = np.append(Ts, pc_tran_flux[i] / straight_tran_flux[i])
Here you have zero in denominator