Problem with the python example contained in the installation directory

[sgrigno2]

Hi all,

First of all, thanks all for all your work that allows me to have a running S4 on my machine. I am working on the example PythonTest.py contained in the S4 folder. I am experiencing inconsistencies between my results and the one published in the "Highly sensitive biosensors based on all-dielectric nanoresonators" by Bontempi and coauthors.

By running the PythonTest.py without any modifications, I got the same results posted by sajidmc (victorliu#61). The problem is that those results are not the ones shown in the publications mentioned above. Trying to obtain a more fair comparison, I slightly modified the PythonTest.py script. In particular, I have changed the material of the bottom layer from Air to Silicon, increased the number of harmonics from 20 up to 120, and the frequencies resolutions (to guarantee convergencies). I am attaching the modified script and the results as well.
As you can see, S4 cannot reproduce the published results in a lower wavelength range. I am wondering why. I am wondering if you can share your thoughts. Any suggestions are more than welcome!

Below the code

CODE:
import S4
import numpy as np
import matplotlib.pyplot as plt
import time

t = time.time()

p = 0.96
d = 0.73
h = 0.22
S = S4.New(Lattice=((p, 0), (0, p)), NumBasis=100)

Define material

S.SetMaterial(Name='SiO2', Epsilon=(1.45 + 0.0j)**2)
S.SetMaterial(Name='Si', Epsilon=(3.487+0.0j)**2)
S.SetMaterial(Name='Vacuum', Epsilon=(1 + 0j)**2)

Define structure

Top semi-infinite latyer

S.AddLayer(Name='AirAbove', Thickness=0, Material='Vacuum')
#S.AddLayer(Name = 'AirDisp',Thickness = 1, Material = 'Vacuum')

Si-Nnaodisck: Air block plus Si circle

S.AddLayer(Name='Si_disks', Thickness=h, Material='Vacuum')
S.SetRegionCircle(Layer='Si_disks', Material='Si', Center=(0, 0), Radius=d/2)

SiO2 substrate

S.AddLayer(Name='Glass_Below', Thickness=1.99, Material='SiO2')

Semiinfite layer: Semi infite Silicon

S.AddLayer(Name='SiBelow', Thickness=0, Material='Si')

Excitation: plane wave

S.SetExcitationPlanewave(
IncidenceAngles=(
0,# polar angle in [0,180)
0 # azimuthal angle in [0,360)
),
sAmplitude=0,
pAmplitude=1,
Order=0
)

Set options

S.SetOptions(
PolarizationDecomposition=True,
PolarizationBasis='Normal'
)

#wavelength_um = 1;
#freq = 1 / float(wavelength_um);
#S.SetFrequency(freq) #unit 2pic_const / a
#forward,backward = S.GetPowerFlux(Layer = 'AirAbove', zOffset = 0) # reflected power
#forward = S.GetPowerFlux(Layer = 'SiBelow',zOffset = 0)
#print('f='+str(1/freq)+', Power_forward= '+str(forward[0].real)+', Power_backward='+str(-backward.real))

frequency sweep

wavelength_space = np.linspace(1.3, 1.7, 1500)

memory allocation

R = np.zeros((len(wavelength_space)))
T = np.zeros((len(wavelength_space)))
for ii in range(len(wavelength_space)):
lam = wavelength_space[ii]
f = 1 / float(lam)
S.SetFrequency(f)
(forward1, backward1) = S.GetPowerFlux(Layer='AirAbove', zOffset=0)
(forward2, backward2) = S.GetPowerFlux(Layer='SiBelow', zOffset=0)
reflection = -backward1
transmission = forward2
R[ii] = np.real(reflection)
T[ii] = np.real(transmission)

plt.plot(wavelength_space, R, label="R")
plt.plot(wavelength_space, T, label="T")
plt.plot(wavelength_space, T+R, label="T+R")
plt.xlabel('Wavelength (um)')
plt.ylabel('Transmission')
plt.legend(loc="upper left")
plt.ylim(0, 1.1)
plt.show()