Pythonのtorcwaを使ってメタサーフェイス用RCWAの計算をする

import torch

print(torch.__version__) # 2.7.1+cu128
print(torch.cuda.is_available()) # True

!pip install torcwa

!pip install torcwa

# Import
import numpy as np
import torch
from matplotlib import pyplot as plt
import scipy.io
import tqdm

import torcwa

# Hardware
# If GPU support TF32 tensor core, the matmul operation is faster than FP32 but with less precision.
# If you need accurate operation, you have to disable the flag below.
torch.backends.cuda.matmul.allow_tf32 = False
sim_dtype = torch.complex64
geo_dtype = torch.float32
device = torch.device('cuda') #CPUの場合はcpuを選択

# Simulation environment
# light
inc_ang = 0.*(np.pi/180)    # radian
azi_ang = 0.*(np.pi/180)    # radian

# material
substrate_eps = 1.45**2
air = 1.0**2
a_Si = 3.55**2
lam = 1064. # nm

# geometry
L = [500., 500.]            # nm / nm
torcwa.rcwa_geo.dtype = geo_dtype
torcwa.rcwa_geo.device = device
torcwa.rcwa_geo.Lx = L[0]
torcwa.rcwa_geo.Ly = L[1]
torcwa.rcwa_geo.nx = 500
torcwa.rcwa_geo.ny = 500
torcwa.rcwa_geo.grid()
torcwa.rcwa_geo.edge_sharpness = 1000.
z = torch.linspace(-500,1500,101,device=device)

x_axis = torcwa.rcwa_geo.x.cpu()
y_axis = torcwa.rcwa_geo.y.cpu()
z_axis = z.cpu()

# layers
layer0_geometry = torcwa.rcwa_geo.circle(R=90.,Cx=L[0]/2.,Cy=L[1]/2.)
layer0_thickness = 800.

# View layers
plt.imshow(torch.transpose(layer0_geometry,-2,-1).cpu(),origin='lower',extent=[x_axis[0],x_axis[-1],y_axis[0],y_axis[-1]])
plt.title('Layer 0')
plt.xlim([0,L[0]])
plt.xlabel('x (nm)')
plt.ylim([0,L[1]])
plt.ylabel('y (nm)')
plt.colorbar()

# Generate and perform simulation
order_N = 5
order = [order_N,order_N]
pillar_list = torch.linspace(50.,147,50,dtype=geo_dtype,device=device)

txx = []
for pillar_radius in tqdm.tqdm(range(len(pillar_list))):
    sim = torcwa.rcwa(freq = 1/lam, order = order, L = L , dtype = sim_dtype, device = device)
    sim.add_input_layer(eps = substrate_eps)
    sim.set_incident_angle(inc_ang = inc_ang, azi_ang = azi_ang)
    layer0_geometry = torcwa.rcwa_geo.circle(R = pillar_list[pillar_radius], Cx=L[0]/2., Cy=L[1]/2.)
    layer0_eps = layer0_geometry*a_Si + (1.-layer0_geometry)
    sim.add_layer(thickness = layer0_thickness, eps = layer0_eps)
    sim.solve_global_smatrix()
    txx.append(sim.S_parameters(orders=[0,0],direction='forward',port='transmission',polarization='xx',ref_order=[0,0]))
txx = torch.cat(txx)

# View Transmittance and phase
fig,ax = plt.subplots()
ax.plot(pillar_list.cpu(),torch.abs(txx).cpu()**2,label='Transmittance (a.u.)',color='blue')
ax2 = ax.twinx()
ax2.plot(pillar_list.cpu(),np.unwrap(torch.angle(txx).cpu())*180/np.pi,color = 'orange',linestyle='--',label='Phase (rad)')
ax.set_title('Spectrum (order: '+str(order_N)+')')
ax.set_xlabel('Pillar radius (nm)')
ax.set_ylabel('Transmittance (a.u.)',color='blue')
ax2.set_ylabel('Phase (rad)',color='orange')
plt.grid()