Using the Integrate module

In this notebook, we demonstrate how to use the pySCATMECH Integrate module. Its main purpose is to calculate integrated scatter from BRDF models or integrated cross section from local BRDF models over a complex geometry.

In the following, we use some symbols and definitions:

• \(f_\mathrm{r}\) : bidirectional reflectance distibution function (BRDF)

• \(\mathrm{d}\sigma/\mathrm{d}\Omega\) : differential scattering cross section (DSC)

• \(\theta\) : polar angle

• \(\phi\) : azimuth angle

• \(\Omega\) : solid angle

• \(\rho\) : reflectance

• \(\sigma\) : cross section

• \(\theta_\mathrm{i}\) : incident polar angle

When we perform numerical integration of the BRDF over a solid angle to obtain reflectance, we approximate \begin{equation} \rho = \int_\Omega f_\mathrm{r}(\theta_\mathrm{i},\theta,\phi)\;\cos\theta\sin\theta \; \mathrm{d}\theta \; \mathrm{d}\phi \end{equation} with \begin{equation} \rho \approx \sum_j f_\mathrm{r}(\theta_\mathrm{i},\theta_j,\phi_j)\; w_j \end{equation} where \(\theta_j\) and \(\phi_j\) are sampled polar and azimuth angles, and \(w_j\) are weighting factors.

For local models that define the differential scattering cross section (DSC), we approximate the integrated cross section \(\sigma\) as \begin{equation} \sigma = \int_\Omega \frac{\mathrm{d}\sigma(\theta_\mathrm{i},\theta,\phi)}{\mathrm{d}\Omega}\;\sin\theta \; \mathrm{d}\theta \; \mathrm{d}\phi \end{equation} with \begin{equation} \rho \approx \sum_j \frac{\mathrm{d}\sigma(\theta_\mathrm{i},\theta_j,\phi_j)}{\mathrm{d}\Omega}\; \frac{w_j}{\cos\theta_j} \end{equation}

In this notebook, we will be demonstrating how to integrate the BRDF to obtain the reflectance or the DSC to obtain cross section.

We always start by importing the libraries we need. We are going to do some graphing, so we need matplotlib.pyplot, and since we are demonstrating it, we will need pySCATMECH.integrate.

import matplotlib.pyplot as plt
from pySCATMECH.integrate import *

Defining geometries

We can define various solid angles. Let’s start with a hemisphere:

hemi = Hemisphere()

The Integrator class will be used to perform integration. Its constructor takes a point spacing and a detector geometry.

The Integrator.PlotSamplingPoints() function plots the integration points in a manner in which the symbol diameter is proportional to the weighting factor. There are some optional parameters that we will discuss later.

There are three different types of samplings available; type = 1 is the default. For type=, points are spaced evenly in angle space. For type=2, points are spaced evenly in solid angle. For type=3, points are evenly spaced in projected solid angle.

Here, we show the three different integration types for the hemispherical detector.

Notice that pySCATMECH provides deg, which is \(\pi/180\). Also, notice that the sampling weights become smaller where the density of points becomes greater.

integrator1 = Integrator(2*deg, hemi, type=1)
integrator1.PlotSamplingPoints()

integrator2 = Integrator(2*deg, hemi, type=2)
integrator2.PlotSamplingPoints()

integrator3 = Integrator(2*deg, hemi, type=3)
integrator3.PlotSamplingPoints()

Let’s build some more interesting detection geometries. We start by demonstrating some builtin shapes. One can collect radiation over a circular cone. Here, the cone is centered on \(\theta=45^\circ\) and \(\phi=90^\circ\) and has a half angle of \(20^\circ\).

cone = CircularCone(theta=45*deg, phi=90*deg, alpha=20*deg)

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

Here, we have an elliptical cone. The half angle in the x direction is 20°, the half angle in the y direction is 45°.

cone = EllipticalCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg))

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

We can rotate the previous elliptical cone by angle \(\gamma\), in this case 45°.

cone = EllipticalCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg), gamma=45*deg)

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

We can have a rectangular cone. In this case, it is centered on the surface normal (\(\theta=0,\phi=0\)) and has half angles of 20° and 45°.

cone = RectangularCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg))

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

And we can rotate that rectangular cone by \(\gamma=45^\circ\):

cone = RectangularCone(theta=0*deg, phi=0*deg, alpha=(20*deg, 45*deg), gamma=45*deg)

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

We can also define arbitrary paths in projected cosine space with a list of \((x,y)\) pairs. The pairs must represent consecutive vertices around the perimeter. Pretty complicated geometries can be generated this way.

cone = ProjectedPolygon([(0.5,0.5),
                         (0.5,-0.5),
                         (0,-0.2),
                         (-0.5,-0.5),
                         (-0.5,0.5),
                         (0,0.2)])

integrator = Integrator(1*deg, cone)
integrator.PlotSamplingPoints()

We can also combine shapes with logical operations & (and), | (or), and ~ (not). Thus, we can place holes to let radiation enter and exit the hemisphere.

hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)

detector =  Hemisphere() & ~(hole1 | hole2)

integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()

Here is a more complex geometry, created by combining different elements.

circle = CircularCone(theta=0, phi=0, alpha=70*deg)
hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)
hole3 = CircularCone(theta=0, phi=0, alpha=30*deg)
center = CircularCone(theta=0, phi=0, alpha=20*deg)

detector = circle & ~(hole1 | hole2 | hole3) | center

integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()

Yet another example:

slot = RectangularCone(theta=0, phi=0, alpha=(90*deg, 6*deg))

detector3 = ~slot

integrator = Integrator(1*deg, detector3)
integrator.PlotSamplingPoints()

Integrate reflectance

So, let’s integrate a BRDF model. We start by creating a model with a set of parameters. We are going to use Microroughness_BRDF_Model with a gaussian rough surface. Then, we integrate it over the hemisphere (with a port provided to let the incident radiation enter and specularly reflected radiation exit) using Reflectance.

# The model parameters as a dictionary...
parameters = {'lambda' : 0.532,
              'substrate' : 4.05+0.05j,
              'type' : 0,
              'psd' : {None : 'Gaussian_PSD_Function',
                       'sigma' : 0.05,
                       'length' : 1
                       }
             }

# Define the model and set its parameters...
model = BRDF_Model("Microroughness_BRDF_Model", **parameters)

# Define the incident angle...
thetaIncident = 6*deg

# The detector is the hemisphere with a hole cut out of it to let in the
# incident light and to let out the specular reflection...
detector = Hemisphere() & ~CircularCone(theta=0*deg, alpha=10*deg)

# Get the integration points and show them...
integrator = Integrator(1*deg, detector)
integrator.PlotSamplingPoints()

# Calculate and print the reflectance ...
ref = integrator.Reflectance(model, thetaIncident)
print("Integrated reflectance = %g" % ref)

Integrated reflectance = 0.239009

We can see what this looks like as a function of a parameter. In the following, we plot the integrated reflectance as a function of the root-mean-square (RMS) roughness. Notice that we space the points out logarithmically. Also, note how we can define the incident polarization.

# Get the RMS roughnesses in logarithmic spacing...
sigmas = np.exp(np.linspace(np.log(0.0001), np.log(.01), 20))

# Initialize a list of reflectances...
refs = []

# s-polarized incident radiation...
spol = StokesVector(1, 1, 0, 0)

for sigma in sigmas:
    # Set the RMS roughness...
    model.setParameter("psd.sigma", sigma)

    # Calculate the reflectance...
    ref = integrator.Reflectance(model, thetaIncident, incpol=spol)

    # Add it to the list..
    refs.append(ref)

# Plot it...
plt.figure()
plt.loglog(sigmas,refs)
plt.xlabel("$\sigma$ (µm)")
plt.ylabel("Reflectance")
plt.title("Reflectance versus RMS Roughness")
plt.show()

Integrate cross section

Models that inherit Local_BRDF_Model calculate the differential scattering cross section (DSC) for individual particles or localized defects. We can use IntegrateCrossSection to obtain the scattering cross section. For this example, we use the Bobbert_Vlieger_BRDF_Model, which is an accurate theory for the scattering by a sphere above a surface.

BVparameters = {'lambda' : 0.532,
              'substrate' : "silicon",
              'type' : 0,
              'sphere' : 1.59,
              'radius' : 0.05,
              'spherecoat' : "No_StackModel",
              'stack' : {None : 'SingleFilm_StackModel',
                        'material' : '(1.5,0)',
                        'thickness' : '0.0016'},
              'delta' : 0} # use defaults for operational parameters

BVmodel = Local_BRDF_Model("Bobbert_Vlieger_BRDF_Model",**BVparameters)

thetaIncident = 70*deg

circle = CircularCone(theta=0, phi=0, alpha=70*deg)
hole1 = CircularCone(theta=70*deg, phi=0, alpha=5*deg)
hole2 = CircularCone(theta=70*deg, phi=180*deg, alpha=5*deg)
hole3 = CircularCone(theta=0*deg, phi=0, alpha=30*deg)

detector = circle & ~(hole1 | hole2 | hole3)

integrator = Integrator(2*deg, detector)
integrator.PlotSamplingPoints()

cs = integrator.CrossSection(BVmodel,thetaIncident,incpol=[1,-1,0,0])

print("Cross section = %g µm^2" % cs)

Cross section = 0.000183612 µm^2

We can see what the cross section is as a function of sphere diameter.

diameters = np.exp(np.linspace(np.log(0.05),np.log(2),20)) #logarithmic spacing

refp = []
refs = []
spol = Polarization('s')
ppol = Polarization('p')
for d in diameters:
    BVmodel.setParameters(radius=d/2)
    refs.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=spol) )
    refp.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=ppol) )

plt.figure()
plt.loglog(diameters,refs,label="s")
plt.loglog(diameters,refp,label="p")
plt.xlabel("Diameter / $\mathrm{\mu m}$")
plt.ylabel("Cross section / $\mathrm{\mu m}^2$")
plt.title("Cross Section versus Diameter")
plt.legend()
plt.show()

Visualizing intensity

We can visualize the scatter distribution by using the PlotSamplingPoints with some arguments. The function DSC returns an array of differential scattering cross sections for each of the sampling points and we can pass that to the color argument of PlotSamplingPoints. Because the default symbol size for displaying the sample points leaves gaps between the symbols, we can expand their size to fill those gaps using the expand argument.

hemi = Hemisphere()
integrator = Integrator(1*deg, hemi)
dsc = integrator.DSC(BVmodel, thetaIncident, incpol=ppol)
integrator.PlotSamplingPoints(color=np.log(dsc), expand=3)

Polarization

This package is designed to simulate integrated scatter for directional incident irradiation. The incident polarization is therefore pretty straightward to define and has been illustrated above using Stokes vectors. The detector system, however, may also be polarization sensitive. We use the sensitivity argument of the SolidAngle to define that sensitivity. The Sensitivity function differs from the Polarization function, because a factor of two is usually needed when a detector is sensitive to a single polarization and not the other.

detector = CircularCone(theta=45*deg, phi=90*deg, alpha=30*deg, sensitivity=Sensitivity('s'))
integrator = Integrator(1*deg, detector)

diameters = np.exp(np.linspace(np.log(0.05), np.log(2),20)) #logarithmic spacing

ref = []
for d in diameters:
    BVmodel.setParameters(radius=d/2)
    ref.append( integrator.CrossSection(BVmodel, thetaIncident, incpol=Polarization('p')) )

plt.figure()
plt.loglog(diameters, ref)
plt.xlabel("Diameter / $\mathrm{\mu m}$")
plt.ylabel("Cross section / $\mathrm{\mu m}^2$")
plt.title("Cross Section versus Diameter")
plt.show()