""" How distribution width changes the scattering phase function ============================================================ Compare lognormal number distributions at the same median diameter. The right panel shows the unpolarized phase function per unit solid angle, normalized so its integral over the sphere is one. Intensities are averaged before normalization; averaging complex amplitudes or equally weighted normalized single-particle phase functions would give a different result. ``diameter_sampling`` controls the number of diameter quadrature nodes, not the number of physical particles. ``angular_sampling`` independently controls the resolution of the plotted scattering angles. Increase diameter sampling to check convergence of the ensemble phase function. This uses the non-interacting approximation: no interparticle electromagnetic coupling, multiple scattering, or interparticle interference. Refer to PackLab for correlation-based dependent-scattering approximations when correlations matter. The comparison isolates width at fixed number median; the arithmetic mean diameter and particle volume therefore also change with width. """ import matplotlib.pyplot as plt import numpy as np from PyMieSim import Measure, ParticleSizeDistribution, PlaneWave, PolarizationState, Simulation, Sphere, ureg MEDIAN_DIAMETER = 1.0 * ureg.micrometer WAVELENGTH = 532 * ureg.nanometer PARTICLE_INDEX = 1.59 MEDIUM_INDEX = 1.33 WIDTHS = (1.0, 1.05, 1.15, 1.30) DIAMETER_SAMPLING = 128 ANGULAR_SAMPLING = 901 def ensemble_phase_function(distribution, angles): """Return an incoherent, normalized unpolarized phase function and mean g. The two orthogonal polarization intensities are averaged for each sphere. All diameters share wavelength and medium, so the common amplitude scale cancels during normalization. Gauss-Legendre angular integration checks the normalization independently of the displayed angular grid. """ source = PlaneWave( wavelength=WAVELENGTH, amplitude=1 * ureg.volt / ureg.meter, polarization=PolarizationState(angle=0 * ureg.degree), ) cosine, angular_weights = np.polynomial.legendre.leggauss(512) evaluation_angles = np.concatenate([angles, np.arccos(cosine)]) * ureg.radian mean_intensity = np.zeros(len(evaluation_angles)) scattering_weights = [] asymmetries = [] for diameter, fraction in zip(distribution.diameters, distribution.number_fractions): simulation = Simulation( source=source, scatterer=Sphere(diameter=diameter, material=PARTICLE_INDEX, medium=MEDIUM_INDEX), ) s1, s2 = simulation.get_s1s2(evaluation_angles) intensity = (np.abs(s1.magnitude) ** 2 + np.abs(s2.magnitude) ** 2) / 2 mean_intensity += fraction * intensity scattering_weights.append(fraction * simulation.run(Measure.CSCA).magnitude) asymmetries.append(simulation.run(Measure.G).magnitude) integrated = mean_intensity[len(angles):] normalization = 2 * np.pi * np.dot(angular_weights, integrated) phase = mean_intensity[:len(angles)] / normalization g_from_angles = 2 * np.pi * np.dot(angular_weights, cosine * integrated) / normalization g_from_cross_sections = np.dot(scattering_weights, asymmetries) / np.sum(scattering_weights) np.testing.assert_allclose(g_from_angles, g_from_cross_sections, atol=2e-7, rtol=0) # The unscaled angular integral must agree with the Mie scattering cross # section up to the common wavenumber-squared factor for all diameters. wavenumber = 2 * np.pi * MEDIUM_INDEX / WAVELENGTH.to('meter').magnitude np.testing.assert_allclose(normalization / wavenumber ** 2, np.sum(scattering_weights), rtol=2e-7, atol=0) return phase, g_from_angles def make_figure(diameter_sampling=DIAMETER_SAMPLING, angular_sampling=ANGULAR_SAMPLING): """Plot the number distributions and their ensemble phase functions.""" angles = np.linspace(0, np.pi, angular_sampling) diameters_um = np.linspace(0.2, 2.5, 1000) colors = ('#253449', '#007f86', '#4676c4', '#d26435') phase_functions = [] with plt.rc_context({'font.size': 11, 'axes.spines.top': False, 'axes.spines.right': False}): figure, (size_axis, phase_axis) = plt.subplots( 1, 2, figsize=(13, 6.3), gridspec_kw={'width_ratios': [1, 1.85]}, ) figure.subplots_adjust(left=0.07, right=0.98, bottom=0.19, top=0.76, wspace=0.27) figure.suptitle('How particle-size width changes angular scattering', x=0.07, y=0.97, ha='left', fontsize=20, fontweight='bold', color='#253449') figure.text(0.07, 0.88, 'Same number median: 1 µm | Wavelength: 532 nm | Unpolarized light\n' 'Sphere refractive index: 1.59 · Medium refractive index: 1.33', color='#526071', linespacing=1.6) for width, color in zip(WIDTHS, colors): distribution = ParticleSizeDistribution.lognormal(MEDIAN_DIAMETER, width, sampling=diameter_sampling) phase, asymmetry = ensemble_phase_function(distribution, angles) phase_functions.append(phase) if width == 1: label = 'Single size' size_axis.axvline(1, color=color, linewidth=2, linestyle='--') size_axis.annotate('Single size', xy=(1, 7.5), xytext=(1.25, 8.6), color=color, arrowprops={'arrowstyle': '->', 'color': color}, fontsize=10) else: label = f'Geometric SD = {width:.2f}' sigma = np.log(width) density = np.exp(-0.5 * (np.log(diameters_um) / sigma) ** 2) / (diameters_um * sigma * np.sqrt(2 * np.pi)) size_axis.plot(diameters_um, density, color=color, linewidth=2) size_axis.fill_between(diameters_um, density, color=color, alpha=0.09) phase_axis.semilogy(np.rad2deg(angles), phase, color=color, linewidth=1.8 if width == 1 else 2.3, linestyle='--' if width == 1 else '-', label=f'{label} (g = {asymmetry:.3f})') size_axis.set(title='A Number-based size distributions', xlabel='Particle diameter [µm]', ylabel='Number probability density [µm⁻¹]', xlim=(0.2, 2.5), ylim=(0, 10)) phase_axis.set(title='B Scattering phase function', xlabel='Scattering angle [degrees]', ylabel='Normalized phase function [sr⁻¹]', xlim=(0, 180)) phase_axis.set_xticks(np.arange(0, 181, 30)) phase_axis.legend(loc='upper right', frameon=False, fontsize=10) phase_axis.grid(which='major', alpha=0.15) size_axis.grid(axis='y', alpha=0.15) phase_axis.text(0, -0.17, '0° = forward', transform=phase_axis.transAxes, fontsize=10, color='#526071') phase_axis.text(1, -0.17, '180° = backward', transform=phase_axis.transAxes, ha='right', fontsize=10, color='#526071') figure.text(0.07, 0.035, 'Non-interacting particles only. Intensities are averaged, then normalized over solid angle.\n' f'{diameter_sampling} diameter quadrature points per continuous distribution; ' f'{angular_sampling} plotted angles. The single-size case uses one diameter.', color='#526071', fontsize=9.5, linespacing=1.6) return figure, np.asarray(phase_functions) figure, phase_functions = make_figure() plt.show()