Issue: non-zero derivatives using autograd after grating cutoff

[jlin0351]

• grcwa version: 0.1.2
• Python version: 3.11.4
• Operating System: MacOS

Hi!

Description

I am trying to calculate the derivative of reflection coefficients with respect to incident plane wave angle in a 1D grating structure. I would like to be able to calculate the derivative at a fixed angle (theta = 0.) but with varying wavelength. I am using autograd to calculate the theta derivative, but I notice that the derivative beyond the cutoff of the m = -1 diffraction order (when the wavelength just exceeds the grating period) is non-zero even though the reflection coefficients themselves, calculated with GRCWA, are identically zero.

Example

Here is an example code:

import grcwa
grcwa.set_backend('autograd') 
import autograd.numpy as np
from autograd import grad
import matplotlib.pyplot as plt


def rNeg1(inc_angle, wavelength, grating_grid):    
    # Grating
    dy = 1e-4 # Small slice in the y-direction to simulate 2D grating
    grating_period = 1 
    L1 = [grating_period,0]
    L2 = [0,dy]
    nG = 30

    # Incoming wave
    freq = 1/wavelength 
    Qabs = 1e5
    freqcmp = freq*(1+1j/2/Qabs)
    theta = inc_angle # radians
    phi = 0.

    # setup RCWA
    obj = grcwa.obj(nG,L1,L2,freqcmp,theta,phi,verbose=0)


    ## CREATE LAYERS ##
    # Grid resolution
    Nx = 30
    Ny = 1 

    # Layer depth
    eps_vacuum = 1
    vacuum_depth = grating_period
    grating_depth = grating_period

    # Construct layers
    obj.Add_LayerUniform(vacuum_depth,eps_vacuum)
    obj.Add_LayerGrid(grating_depth,Nx,Ny)
    obj.Add_LayerUniform(vacuum_depth,eps_vacuum)
    obj.Init_Setup()

    # Construct patterns
    obj.GridLayer_geteps(grating_grid)


    ## INCIDENT WAVE ##
    planewave={'p_amp':0,'s_amp':1,'p_phase':0,'s_phase':0}
    obj.MakeExcitationPlanewave(planewave['p_amp'],planewave['p_phase'],planewave['s_amp'],planewave['s_phase'],order = 0)

    # solve for R and T
    R_byorder = obj.RT_Solve(normalize=1, byorder=1)[0]
    Fourier_orders = obj.G

    order = -1
    return np.sum(R_byorder[Fourier_orders[:,0]==order])


# Use autograd to calculate theta derivative
PDrNeg1 = grad(rNeg1)


# Plot near cutoff
inc_angle = 0.
wavelength = 0.7
grating_grid = [1]*10 + [9]*10 + [1]*10

n = 100
vals = np.zeros(n, dtype=float)
wavelengths = np.linspace(0.999,1.001,n)
for idx,lam in enumerate(wavelengths):
    vals[idx] = PDrNeg1(inc_angle, lam, grating_grid)

fig,ax = plt.subplots(1)
ax.plot(wavelengths,vals, lw='2', color='r')
ax.set(xlabel="Wavelength/grating period", ylabel=r"$\partial r_{-1}(0)/\partial \theta'$")
ax.axvline(x=1, color='k', linewidth=0.5)
ax.axhline(y=0, color='k', linewidth=1)

and the plot output on my end:

I've tried inputting a grating with randomised permittivities in grating_grid, and increasing Nx and nG but the result is still a non-zero derivative after cutoff. I was wondering if you would have any insights into this behaviour?

Thank you for your time, and thanks for making this package!

[weiliangjinca]

    # Incoming wave
    freq = 1/wavelength 
    Qabs = 1e5
    freqcmp = freq*(1+1j/2/Qabs)

Note that you are solving Maxwell's equation at a complex frequency, so the reflection is likely nonzero beyond the cutoff. For this value of Qabs, it's normal to see this level of low reflection. If you further increase Qabs, I think you can see the derivative drops.

[jlin0351]

Ah I see, that makes sense. Thank you!