Simulation of a layer of a substance absorbing at a certain frequency

[leSemaleon]

Hello,

I have material that is not among the standard ones. I wrote an approximation for it using functions. But the transmission spectrum does not correspond to reality. The transmission spectrum is as if there is a transparent material in the path of the radiation, that is, it is equal to 100 percent everywhere in the tested frequency range. But this shouldn't happen.
The pictures show the imaginary and real parts of epsilon for the material and the transmission spectrum


The real and imaginary parts of epsilon are obtained using the method Medium.epsilon using my material function
Layer thickness is 100 microns. The transmission spectrum does not depend on thickness for some reason. I tried with resolution 1, 10, and 100, the result is the same.

[smartalecH]

How did you add your materials and geometry to your meep simulation? Some definitions, ie those involving material functions, require you to additionally pass the materials into the extra_materials field of the simulation constructor.

[leSemaleon]

I created a class for my material, or rather susceptibility, which inherits from the class in the geom. It also contains all the methods found in susceptibilities like Drude. Using this susceptibility, I set the material using mp.Meidum(epsilon=1, E_Susceptibility = [mySusceptibility])

[smartalecH]

Right but how did you assign it to your geometry? It's usually easier if you just include your code.

[leSemaleon]

Ok. There is a code. 0 Iteration without geometry nd 1 iteration with

def main(iteration=0, refer_refl_data=[]):

    # meep unit = 1 μm
    dpml = 1000
    dAirIn = 200
    dAirOut = 200
    dFilm = 500

    accuracityOfDecay = 1e-9

    sx = 0
    sy = 0 
    sz = dpml + dAirIn + dFilm + dAirOut + dpml

    field_component = mp.Ex

    cell_size = mp.Vector3(sx, sy, sz)
    
    gridstep = 0.1
    resolution = 1/gridstep

    lmin = 100         # source min wavelength
    lmax = 1000         # source max wavelength
    fmin = 1/lmax       # source min frequency
    fmax = 1/lmin       # source max frequency
    fcen = 0.5*(fmin+fmax)
    df = fmax-fmin

    nfreq = NFREQ  # number of frequencies at which to compute flux    

    sources = [
        mp.Source(
            mp.GaussianSource(fcen, fwidth=df),
            component=field_component,
            center=mp.Vector3(0, 0, -0.5*sz+dpml+0.2*dAirIn),
            size=mp.Vector3(sx,sy,0),
        )
    ]
    
    boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z,side=mp.Low),
                      mp.PML(thickness=dpml,direction=mp.Z,side=mp.High)]

    geometry = []
    if iteration:
        film = mp.Block(
                    size = mp.Vector3(sx, sy, dFilm),
                    center = mp.Vector3(0, 0, -0.5*sz+dpml+dAirIn+0.5*dFilm),
                    material = usrmp.RDX
                )
        geometry.append(film)

    sim = mp.Simulation(
        cell_size = cell_size,
        boundary_layers = boundary_layers,
        sources = sources,
        geometry = geometry,
        resolution = resolution,
        dimensions = 1,
    )
          
    refl_fr = mp.FluxRegion(
        center=mp.Vector3(0,0,-sz/2+dpml+0.5*dAirIn), 
        size=mp.Vector3(sx,sy,0)
    )

    refl = sim.add_flux(fcen, df, nfreq, refl_fr)

    tran_fr = mp.FluxRegion(
        center=mp.Vector3(0,0,sz/2-dpml-0.5*dAirOut), 
        size=mp.Vector3(sx,sy,0)
    )

    tran = sim.add_flux(fcen, df, nfreq, tran_fr)
    
    if iteration:
        sim.load_minus_flux_data(refl, refer_refl_data)
    
    sim.run(until_after_sources=mp.stop_when_fields_decayed(30, field_component, 
                                                            mp.Vector3(0,0,sz/2-dpml-0.5*dAirOut), 
                                                            accuracityOfDecay))
    
    refl_data = sim.get_flux_data(refl)
    refl_flux = mp.get_fluxes(refl)
    tran_flux = mp.get_fluxes(tran)
    flux_freqs = mp.get_flux_freqs(refl)

    sim.reset_meep()
    data = [resolution, dFilm]
    return refl_data, refl_flux, tran_flux, flux_freqs, data

[leSemaleon]

I do everything as in the example, except to shorten it by denoting the simulation as a function and calling it twice. I dont know why it doesnt work/

[leSemaleon]

any ideas on how to fix this?

[WeiKuoLi]

Your dipole fields are not updated accordingly. Meep thus think the whole medium is vacuum and gives you 100% transmission. Without knowing your custom code it would not be possible knowing what went wrong.

[leSemaleon]

The previous code is missing a description of RDX.
I'm adding it now

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import meep as mp
import numpy as np
from scipy import special
import math

import sys
sys.path.append('/anaconda3/envs/meep/lib/python3.11/site-packages/meep')
import geom

class RDXSusceptibility(geom.Susceptibility):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        
    def eval_susceptibility(self, freq):
        sigma = np.expand_dims(
            geom.Matrix(diag=self.sigma_diag, offdiag=self.sigma_offdiag), axis=0
        )
        
        f = [0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1, 
              1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, 
              1.75, 1.80, 1.85, 1.90, 1.95, 2, 2.05, 2.10, 2.15, 2.20, 2.25, 2.30]
        n = [1.5+0.048j, 1.5+0.056j, 1.54+0.065j, 1.62+0.077j, 1.7+0.099j, 1.77+0.127j, 
             1.83+0.165j, 1.85+0.231j, 1.76+0.314j, 1.5+0.332j, 1.3+0.275j, 1.24+0.188j, 
             1.29+0.133j, 1.36+0.118j, 1.39+0.113j, 1.4+0.095j, 1.45+0.080j, 1.5+0.080j, 
             1.55+0.092j, 1.57+0.109j, 1.55+0.12j, 1.54+0.121j, 1.53+0.12j, 1.52+0.117j, 
             1.52+0.113j, 1.53+0.108j, 1.54+0.105j, 1.57+0.108j, 1.59+0.117j, 1.6+0.131j, 
             1.6+0.146j, 1.57+0.156j, 1.54+0.157j, 1.52+0.151j, 1.51+0.147j, 1.5+0.14j, 
             1.5+0.137j, 1.5+0.135j, 1.49+0.132j, 1.49+0.124j]
        f0 = freq*3e8/1e-6*1e-12
        
        # n0 = np.interp(f0, f, n)
        n1 = np.interp(f0,f,np.real(n))
        k1 = np.interp(f0,f,np.imag(n))
        
        eps1 = n1**2-k1**2
        eps2 = 2*n1*k1
        
        return (eps1 + 1j*eps2) * sigma

RDX = mp.Medium(epsilon=1, E_susceptibilities=[RDXSusceptibility(sigma=1)])   

The function from the previous code is used using the code below

start = datetime.now()

table = np.zeros((NFREQ, 2))

refer_refl_data, refer_refl_flux, refer_tran_flux, *t = main(0, [])
refl_data, result_refl_flux, result_tran_flux, flux_freqs, data = main(1, refer_refl_data)
fl = []
Ts = []
for j in range(NFREQ):
    fl = np.append(fl, flux_freqs[j]*3e8/1e-6*1e-12)
    Ts = np.append(Ts, result_tran_flux[j]/refer_tran_flux[j] if refer_tran_flux[j]!=0 else 2)
table[:,0] = np.transpose(fl)
table[:,1] = np.transpose(Ts)
    
end2 = datetime.now()