The field obtained by sim.field_xy is not symmetric when the structure is symmetric

[Maxwell-github]

Thanks for providing this useful code.
Recently, I found a problem. The field obtained by sim.field_xy is not symmetric when the structure is symmetric. It has a one-pixel shift. I use a simple symmetric square to do the test. The coordinates are symmtric about the center but the field is not. I think it may have some problem. I also have tested odd sampling points N = 11 but the problem remains. Could you fix this problem?

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

from PIL import Image
import sys
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')

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

# material
PC_eps = 1.560484**2
N15_G1 = 1.506417**2

# DOE input
M = 5 # M * M unit
P = 500 # pixel size: nm
H = 1000 # DOE thickness: nm


# geometry
L = [P*M, P*M]            # nm / nm
N = 10  # number of sampling of each pixel
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 = N*M
torcwa.rcwa_geo.ny = N*M
torcwa.rcwa_geo.grid() 
torcwa.rcwa_geo.edge_sharpness = 1000.

x_axis = torcwa.rcwa_geo.x.cpu()
y_axis = torcwa.rcwa_geo.y.cpu()
lamb0 = torch.tensor(850.,dtype=geo_dtype,device=device)    # nm
z = torch.linspace(lamb0/2,lamb0,1,device=device) 
z_axis = z.cpu()

# Generate and perform simulation
order_N = 10
order = [order_N,order_N]

layer0_geometry = torcwa.rcwa_geo.square(W=P,Cx=L[0]/2,Cy=L[1]/2)

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.show()
layer0_thickness = H
sim = torcwa.rcwa(freq=1/lamb0,order=order,L=L,dtype=sim_dtype,device=device)
sim.add_input_layer(eps=1) # default: vacuum
sim.add_output_layer(eps=PC_eps)
sim.set_incident_angle(inc_ang=inc_ang,azi_ang=azi_ang)
layer0_eps = layer0_geometry*N15_G1 + (1.-layer0_geometry)
sim.add_layer(thickness=layer0_thickness,eps=layer0_eps)
sim.solve_global_smatrix()
sim.source_planewave(amplitude=[1.,0.],direction='forward')
[Ex, Ey, Ez], [Hx, Hy, Hz] = sim.field_xy(layer_num=1,x_axis=torcwa.rcwa_geo.x,y_axis=torcwa.rcwa_geo.y,z_prop=z)
Ex = torch.transpose(Ex,-2,-1).cpu()
Ey = torch.transpose(Ey,-2,-1).cpu()
Ez = torch.transpose(Ez,-2,-1).cpu()
Hx = torch.transpose(Hx,-2,-1).cpu()
Hy = torch.transpose(Hy,-2,-1).cpu()
Hz = torch.transpose(Hz,-2,-1).cpu()

plt.imshow(abs(Ex)-torch.flip(abs(Ex),[0]),cmap='jet',origin='lower',extent=[x_axis[0],x_axis[-1],y_axis[0],y_axis[-1]])
plt.colorbar()
plt.show()