{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using the Scatterer module"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook, we will be describing how to use the `Free_Space_Scatterer` class provided in the Scatterer module. This class provides access to the functions provided by the `Free_Space_Scatterer` class of SCATMECH. These include the models\n",
"* RayleighScatterer\n",
"* RayleighGansSphereScatterer\n",
"* MieScatterer\n",
"* CoatedMieScatterer\n",
"* MultilayerCoatedMieScatterer\n",
"* Optically_Active_Sphere_Scatterer\n",
"* TMatrix_Axisymmetric_Scatterer\n",
"\n",
"First, import the libraries we will need. We will be graphing some data, and using the NumPy package."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from pySCATMECH.scatterer import * "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mie scattering\n",
"\n",
"Let's start by looking at scattering by homogeneous spheres (Mie scatering). We start by creating a model and setting its parameters:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{None : 'MieScatterer',\n",
"'lambda' : '0.532',\n",
"'medium' : '(1,0)',\n",
"'radius' : '0.1',\n",
"'sphere' : '(1.59,0)'}\n"
]
}
],
"source": [
"parameters = {\n",
" 'lambda' : 0.532,\n",
" 'medium' : 1,\n",
" 'radius' : 0.1,\n",
" 'sphere' : 1.59}\n",
"\n",
"model = Free_Space_Scatterer(\"MieScatterer\",parameters)\n",
"\n",
"print(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can evaluate the Mueller matrix differential scattering cross section in any pair of directions, defined by two directional vectors. The basis set for the Mueller matrix is perpendicular and parallel to the plane of scattering, defined by the incident and scattering directions."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 2.30829762e-03 -6.16046580e-04 0.00000000e+00 0.00000000e+00]\n",
" [-6.16046580e-04 2.30829762e-03 0.00000000e+00 0.00000000e+00]\n",
" [ 0.00000000e+00 0.00000000e+00 2.22445889e-03 -2.25210041e-05]\n",
" [ 0.00000000e+00 0.00000000e+00 2.25210041e-05 2.22445889e-03]]\n"
]
}
],
"source": [
"theta = 45*deg\n",
"\n",
"vin = [0,0,1] \n",
"vout = [np.sin(theta),0,np.cos(theta)]\n",
"\n",
"m = model.DifferentialScatteringCrossSection(vin,vout)\n",
"print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot the scattering function:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def dsc(theta):\n",
" vin = [0,0,1]\n",
" vout = [np.sin(theta),0,np.cos(theta)]\n",
" return model.DifferentialScatteringCrossSection(vin,vout)\n",
" \n",
"angles = np.linspace(0,2*pi,180)\n",
"results = [dsc(angle) for angle in angles]\n",
"\n",
"spol = Polarization('S') # perpendicular to the scattering plane\n",
"ppol = Polarization('P') # parallel to the scattering plane\n",
"upol = Polarization('U') # unpolarized \n",
"usens = Sensitivity('U') # unpolarized sensitivity\n",
"\n",
"plt.figure()\n",
"plt.polar(angles, [usens @ r @ upol for r in results], label = \"unpol\")\n",
"plt.polar(angles, [usens @ r @ spol for r in results], label = \"perp\")\n",
"plt.polar(angles, [usens @ r @ ppol for r in results], label = \"par\")\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this next step, we look at a 0.04 µm diameter gold sphere as a function of wavelength. Here, we integrate the scatter and show the total scattering cross section as a function of wavelength. Notice the surface plasmon resonance at around 0.52 µm."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"parameters = {\n",
" 'lambda' : 0.532,\n",
" 'medium' : 1,\n",
" 'radius' : 0.03,\n",
" 'sphere' : 'gold'}\n",
"\n",
"model = Free_Space_Scatterer(\"MieScatterer\", parameters)\n",
"\n",
"def integrate_scatter(model, wavelength, steps):\n",
" \"Integrate `model` at `wavelength`, performing the integral in `steps` steps.\"\n",
" integ = 0\n",
" dtheta = pi / (steps - 1)\n",
" parameters['lambda'] = wavelength\n",
" model.setParameters(parameters)\n",
" vin = [0, 0, 1]\n",
" angles = np.arange(0, pi, dtheta)\n",
" # Integrate angles from 0 to pi...\n",
" for theta in angles:\n",
" vout = [np.sin(theta), 0, np.cos(theta)]\n",
" m = model.DifferentialScatteringCrossSection(vin, vout)\n",
" integ += m[0][0] * np.abs(np.sin(theta))\n",
" integ *= 2 * pi * dtheta\n",
" return integ\n",
"\n",
"wavelengths = np.linspace(0.380,0.780,200)\n",
"\n",
"plt.figure()\n",
"\n",
"# Diameter 0.04 µm...\n",
"D = 0.04\n",
"parameters['radius'] = D/2\n",
"model.setParameters(parameters)\n",
"\n",
"crosssection = [integrate_scatter(model,L,50) for L in wavelengths]\n",
"\n",
"plt.plot(wavelengths, crosssection, label = \"D = %g µm\"%D)\n",
"\n",
"plt.legend(frameon=False)\n",
"plt.xlabel(\"Wavelength [µm]\")\n",
"plt.ylabel(\"Cross section [µm$^2$]\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mie scattering with a coating"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `CoatedMieScatterer` model implements Mie theory extended to coated spheres. Let's see what the absorption, scattering, and extinction cross sections look like for different thicknesses of gold on glass spheres in water."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\r\n",
"\r\n",
"\r\n",
"\r\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"parameters = {\n",
" 'lambda' : 0.532,\n",
" 'medium' : 'water', \n",
" 'radius' : 0.03,\n",
" 'sphere' : 'glass',\n",
" 'coating' : 'gold',\n",
" 'thickness' : 0.005}\n",
"\n",
"model = Free_Space_Scatterer(\"CoatedMieScatterer\", parameters)\n",
"\n",
"def extinction(model, wavelength):\n",
" \"Calculate the extinction coefficient at a given `wavelength` for a given `model`.\"\n",
" parameters['lambda'] = wavelength\n",
" model.setParameters(parameters)\n",
" vin = [0, 0, 1]\n",
" return model.Extinction(vin)[0][0]\n",
"\n",
"def integrated_scatter(model,wavelength,steps):\n",
" \"Integrate the scatter over polar angle\"\n",
" integ = 0\n",
" dtheta = pi/(steps - 1)\n",
" parameters['lambda'] = wavelength\n",
" model.setParameters(parameters)\n",
" vin = [0, 0, 1]\n",
" for theta in np.linspace(0, pi, steps):\n",
" vout = [np.sin(theta), 0, np.cos(theta)]\n",
" m = model.DifferentialScatteringCrossSection(vin, vout)\n",
" integ += m[0][0]*np.abs(np.sin(theta))\n",
" integ *= 2*pi*dtheta\n",
" return integ\n",
"\n",
"wavelengths = np.linspace(0.500, 1.300, 200)\n",
"\n",
"fig = plt.figure(figsize=(10,5))\n",
"ax1 = fig.add_subplot(131)\n",
"ax2 = fig.add_subplot(132)\n",
"ax3 = fig.add_subplot(133)\n",
"\n",
"for t in [0.002,0.003,0.004,0.005,0.006,0.007,0.008]:\n",
" parameters['thickness'] = t\n",
" sca = np.array([integrated_scatter(model,L,50) for L in wavelengths])\n",
" ext = np.array([extinction(model,L) for L in wavelengths])\n",
" abs = ext-sca\n",
" ax1.plot(wavelengths, sca*1000, label = \"t = %g\"%t)\n",
" ax2.plot(wavelengths, ext*1000, label = \"t = %g\"%t)\n",
" ax3.plot(wavelengths, abs*1000, label = \"t = %g\"%t)\n",
" \n",
"ax1.legend()\n",
"ax1.set_title(\"Scattering\")\n",
"ax2.set_title(\"Extinction\")\n",
"ax3.set_title(\"Absorption\")\n",
"ax1.set_ylabel(\"Cross Section [µm$^2$/1000]\")\n",
"ax1.set_xlabel(\"Wavelength [µm]\")\n",
"ax2.set_xlabel(\"Wavelength [µm]\")\n",
"ax3.set_xlabel(\"Wavelength [µm]\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"scrolled": false
},
"source": [
"## Optically active sphere\n",
"\n",
"In the following, we look at the extinction matrix for a 2 µm diameter sphere that has a small index difference for right-hand-circular polarization (typical numbers for sucrose). The resulting matrix shows some optical rotation and some circular diattenuation."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 8.32522992e+00 0.00000000e+00 0.00000000e+00 7.94028248e-06]\n",
" [ 0.00000000e+00 8.32522992e+00 3.05860149e-06 -0.00000000e+00]\n",
" [ 0.00000000e+00 -3.05860149e-06 8.32522992e+00 0.00000000e+00]\n",
" [ 7.94028248e-06 0.00000000e+00 -0.00000000e+00 8.32522992e+00]]\n"
]
}
],
"source": [
"parameters = {\n",
" 'lambda' : 0.589,\n",
" 'medium' : 1,\n",
" 'radius' : 1,\n",
" 'materialright' : 1.56+ 6.9E-7,\n",
" 'materialleft' : 1.56}\n",
"\n",
"model = Free_Space_Scatterer(\"Optically_Active_Sphere_Scatterer\",parameters)\n",
"\n",
"def extinction(model,wavelength):\n",
" model.setParameters(wavelength=wavelength)\n",
" vin = [0,0,1]\n",
" return model.Extinction(vin)\n",
"\n",
"extmatrix = model.Extinction([0,0,1])\n",
"\n",
"print(extmatrix)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Ellipsoids"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the following, we calculate the extinction coefficient for ellipsoidal gold particles, where the major axis is 20 µm and the minor axis is 6 µm. We then calculate the Mueller matrix transmittance assuming a concentration-thickness product of 105 µm-2. When the radiation is incident perpendicular to the axis of the ellipsoids, the resulting Mueller matrix behaves very much like a polarizer. When the radiation is incident along the axis of the ellipsoids, the resulting Mueller matrix is just absorbing."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Mueller matrix transmittance perpendicular to ellipsoid axis:\n",
"[[ 4.90951752e-01 4.90951752e-01 8.56709626e-24 6.07946914e-24]\n",
" [ 4.90951752e-01 4.90951752e-01 8.55781735e-24 6.05189921e-24]\n",
" [-1.63365318e-22 -1.63342992e-22 -3.18191980e-06 7.75832676e-06]\n",
" [-1.47979786e-21 -1.47981003e-21 -7.75832676e-06 -3.18191980e-06]] \n",
"\n",
"Mueller matrix transmittance parallel to ellipsoid axis:\n",
"[[ 1.01786040e+00 -2.45123552e-07 -2.47343011e-17 6.68348925e-19]\n",
" [-2.45123552e-07 1.01786040e+00 4.12590470e-18 2.27875239e-19]\n",
" [-2.47343002e-17 -4.12589871e-18 1.01786040e+00 -1.83673418e-07]\n",
" [ 6.68344506e-19 -2.27876151e-19 1.83673418e-07 1.01786040e+00]] \n",
"\n"
]
}
],
"source": [
"shape = {None : 'Ellipsoid_Axisymmetric_Shape',\n",
" 'npoints' : '100',\n",
" 'vertical' : '0.010',\n",
" 'horizontal' : '0.003'\n",
" }\n",
"\n",
"tmatrix = {'lambda' : 0.8,\n",
" 'medium' : 1.4,\n",
" 'Shape' : shape,\n",
" 'particle' : 'gold',\n",
" 'lmax' : '30',\n",
" 'mmax' : '0'}\n",
"\n",
"model = Free_Space_Scatterer(\"TMatrix_Axisymmetric_Scatterer\",tmatrix)\n",
"\n",
"extCoeff100 = model.Extinction([1,0,0])\n",
"extCoeff001 = model.Extinction([0,0,1])\n",
"\n",
"print(\"Mueller matrix transmittance perpendicular to ellipsoid axis:\")\n",
"print(MuellerExp(-extCoeff100 * 10**5),'\\n')\n",
"print(\"Mueller matrix transmittance parallel to ellipsoid axis:\")\n",
"print(MuellerExp(-extCoeff001 * 10**5),'\\n')\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}