Integral of of a vector field component field in a dielectric

[MKAbdElrahman]

I want to calculate the mode confinement for a slot mode, so I need to integrate the z-component of the poynitng field in the slot region and divide it by the same quantity in the whole domain.

I hope if there is a function that works like compute_energy_in_dielectric named integrate_field_in_dielectric

Here is my MPB code:

import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt

def slot_simulation(wavelength , n_core, n_clad , w,h,s,resolution = 100,sc_x = 2,sc_y = 2,num_modes = 2):
    omega      = 1/wavelength

    Si   = mp.Medium(index= n_core)
    SiO2 = mp.Medium(index= n_clad)
    
    geometry_lattice = mp.Lattice(size=mp.Vector3(sc_x, sc_y))
    geometry = [mp.Block(size=mp.Vector3(w, h,mp.inf),
                     center=mp.Vector3(+0.5*(s+w)), material=Si),
            mp.Block(size=mp.Vector3(w, h,mp.inf),
                     center=mp.Vector3(-0.5*(w+s)), material=Si)]
    ms = mpb.ModeSolver(
    geometry_lattice = geometry_lattice,
    geometry         = geometry,
    resolution       = resolution,
    num_bands        = num_modes,
    default_material = SiO2
    )
    ms.run()
    E = []
    energy=[]
    def store_fields(ms_temp, mode):
        E.append(ms_temp.get_efield(which_band=mode,bloch_phase=False))
        ms_temp.get_hfield(which_band=mode,bloch_phase=False)
        ms_temp.compute_field_energy()
        energy.append(ms_temp.compute_energy_in_objects([mp.Block(size=mp.Vector3(s, h),
                                                                  center=mp.Vector3(), material=SiO2)]))
        
    k = ms.find_k(
    mp.NO_PARITY, # p
    omega, # omega
    1, # band_min
    num_modes, # band_max
    mp.Vector3(0,0,1), # korig_and_kdir
    1e-4, # tol
    omega * 2.0, # kmag_guess
    omega * 0.1, # kmag_min
    omega * 4, # kmag_max
    store_fields, # band_funcs,
    )
    return ms, np.array(k)/omega, E
    

w = 0.180    # waveguide width
h = 0.300    # waveguide height 
s = 0.05     # slot width

n_core = 3.48 # core refractive index
n_clad = 1.44 # clad refractive index

wavelength = 1.55 # free space wavelength
num_modes = 2   # number of modes at the free space wavelength

ms,neff,E = slot_simulation(wavelength,n_core, n_clad,w,h, s,resolution = 100,sc_x = 2,sc_y = 2,num_modes = 2)

[stevengj]

If you only have two dielectric materials in the integration region, e.g. ε = 12 and ε = 1, just integrate the field multiplied by (ε–1), and then divide the result of the integral by (12–1).

[MKAbdElrahman]

I don't understand why this would work. I want to calculate the field confinement factor defined as the ratio of the time averaged power flow in the domain of interest to the time averaged power flow inside the total domain. My domain of interest is the slot between the two wave-guides.

I think a function that takes the current loaded field and a direction(X,Y,Z) and a domain defined by a set of geometrical shapes is important to add for MPB.
I did the integration by defining a mask function that returns true for all points in the integration domain

def mask(ms,obj):
    cell_size =  ms.geometry_lattice.size
    lx = cell_size[0]
    ly = cell_size[1]
    Nx , Ny  = ms.get_epsilon().shape
    x = np.linspace(-lx/2, lx/2, Nx)
    y = np.linspace(-ly/2, ly/2, Ny)
    xv, yv = np.meshgrid(x, y)
    r = [mp.Vector3(xv[i,j],yv[i,j]) for j in range(0,Ny) for i in range(0,Nx)]
    M  = [mp.is_point_in_object(ri,obj)  for ri in r ]
    return np.reshape(M,(Nx,Ny))

then, a direct sum

def integrate(ms,field,obj):
    M = mask(ms,obj)
    return (np.sum(field[M]) / np.sum(field))

[stevengj]

The point is that if you pass Dielectric to integrate_field_function (similar to this example then it will give you the ε at each point, so then you put ε(x)-1 inside your integrand, which is 0 in the air. Or, conversely, if you want the air region you can use ε(x)-εₕᵢ.

Of course, if your "slot" is just a rectangular block, then you can pass that volume directly to integrate_field_function as the where argument.

[stevengj]

(And you could also call mp.is_point_in_object(ri,obj) inside integrate_field_function.)

[MKAbdElrahman]

integrate_field_function is a method for Simulation type objects. I have only a ModeSolver object. The MPB documentation says there is function ModeSolver.compute_field_integral(f) where

• f is a function f(F, eps, r) that returns a number, possibly complex,
• F is the complex field vector at the given point.
  So, I can use it like this

slot_region =mp.Block(size = mp.Vector3(s,h)) 
all_domain =mp.Block(size = mp.Vector3(2,2))  
ms.get_poynting(2)
def f_all(F, eps, r):
    if mp.is_point_in_object(r,all_domain):
        return np.real(F.z)
    else:
        return 0   
def f_slot(F, eps, r):
    if mp.is_point_in_object(r,slot_region):
        return np.real(F.z)
    else:
        return 0     
ms.compute_field_integral(f_slot)/ms.compute_field_integral(f_all)

although, the result doesn't agree with my custom integrate function.

[stevengj]

Oh sorry, I missed that you were talking about MPB.