diff --git a/.github/workflows/lint_and_test.yml b/.github/workflows/lint_and_test.yml new file mode 100644 index 0000000..c98ac0a --- /dev/null +++ b/.github/workflows/lint_and_test.yml @@ -0,0 +1,55 @@ +name: Lint, then test on all platforms + +on: + # empty "push:" will trigger CI on push to any branch + push: + pull_request: + branches: [ "develop", "master" ] + +jobs: + lint: + # only need to run on one platform since we're just looking at the code + runs-on: ubuntu-latest + defaults: + run: + shell: bash -el {0} + steps: + - uses: actions/checkout@v4 + - uses: astral-sh/ruff-action@v3 + with: + # check based on rulesets defined in pyproject.toml and fail if errors + args: "check" + # do another run to produce a linting report for pull requests. exit-zero + # treats all errors as warnings. This will flag codestyle problems in the + # PR + - run: ruff check --exit-zero --output-format=github + + test: + # linting must succeed for testing to run; this helps us rapidly flag code + # errors before going to testing + needs: lint + runs-on: ${{ matrix.os }} + # conda enviroment activation requires bash + defaults: + run: + shell: bash -el {0} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest, windows-latest] + python-version: ["3.13"] # , "3.10", "3.11", "3.12"] + max-parallel: 5 + + steps: + - uses: actions/checkout@v4 + - name: Set up miniforge with pymie environment + uses: conda-incubator/setup-miniconda@v3 + with: + environment-file: environment.yml + miniforge-version: latest + conda-remove-defaults: "true" + - name: Install plugin to annotate test results + run: pip install pytest-github-actions-annotate-failures + - name: Test with pytest + run: | + pytest diff --git a/README.md b/README.md index 209b283..5d4fdf5 100644 --- a/README.md +++ b/README.md @@ -3,9 +3,31 @@ Pure Python package for Mie scattering calculations. This package is used by, for example, the [structural-color python package](https://github.com/manoharan-lab/structural-color). +To install: + +```shell +git clone https://github.com/manoharan-lab/python-mie.git +python -m pip install ./python-mie +``` + +To remove: + +```shell +python -m pip uninstall pymie +``` + +You can then run the tests to make sure the package works properly: +```shell +cd python-mie +pytest -v +``` + The original code was developed by Jerome Fung in the [Manoharan Lab at Harvard -University](http://manoharan.seas.harvard.edu). This research was supported by -the National Science Foundation under grant numbers CBET-0747625, DMR-0820484, -DMR-1306410, and DMR-1420570. The code has since been updated. It now works in -Python 3, and it can handle quantities with dimensions (using +University](http://manoharan.seas.harvard.edu). The code has since been updated +by Victoria Hwang, Anna B. Stephenson, and Vinothan N. Manoharan. It works +in Python 3, and it can handle quantities with dimensions (using [pint](https://github.com/hgrecco/pint)). + +The development of this code was made possible with support from the National +Science Foundation under award numbers CBET-0747625, DMR-0820484, DMR-1306410, +DMR-1420570, and DMR-2011754. diff --git a/environment.yml b/environment.yml index f136943..c0755f6 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,7 @@ -# dependencies for structural color package +# dependencies for pymie package # # To use: -# conda env create -f .\environment.yml +# conda env create -f environment.yml # and then # conda activate pymie # @@ -12,7 +12,6 @@ name: pymie channels: - conda-forge - - defaults dependencies: - python>=3.11 - numpy diff --git a/pymie/__init__.py b/pymie/__init__.py index c5d70bc..57c627c 100644 --- a/pymie/__init__.py +++ b/pymie/__init__.py @@ -28,8 +28,8 @@ Github: https://github.com/hgrecco/pint Docs: https://pint.readthedocs.io/en/latest/ -.. moduleauthor :: Vinothan N. Manoharan -.. moduleauthor :: Sofia Magkiriadou . +.. moduleauthor:: Vinothan N. Manoharan +.. moduleauthor:: Sofia Magkiriadou . """ import numpy as np @@ -68,20 +68,22 @@ def index_ratio(n_particle, n_matrix): Parameters ---------- - n_particle: structcol.Quantity [dimensionless] - refractive index of particle at a particular wavelength + n_particle: structcol.Quantity [dimensionless] or ndarray thereof + refractive index of particle at particular wavelength(s) can be complex n_matrix: structcol.Quantity [dimensionless] refractive index of matrix at a particular wavelength Notes ----- - Returns a scalar rather than a Quantity because scipy special functions - don't seem to be able to handle pint Quantities + Nondimensionalizes from input arguments and strips units, returning a pure + ndarray (not a Quantity object) Returns ------- - float + ndarray or scalar (complex or float): + Return type depends on type of n_particle and n_matrix, and return + shape should be the same as n_particle """ return (n_particle/n_matrix).magnitude @@ -92,26 +94,42 @@ def size_parameter(wavelen, n_matrix, radius): Parameters ---------- - wavelen: structcol.Quantity [length] + wavelen: structcol.Quantity [length], array-like wavelength in vacuum - n_matrix: structcol.Quantity [dimensionless] + n_matrix: structcol.Quantity [dimensionless], array-like refractive index of matrix at wavelength=wavelen - radius: structcol.Quantity [length] + radius: structcol.Quantity [length], array-like radius of particle Notes ----- - Returns a scalar rather than a Quantity because scipy special functions - don't seem to be able to handle pint Quantities + Nondimensionalizes from input arguments and strips units, returning a pure + ndarray (not a Quantity object) Returns ------- - complex or float + ndarray or scalar (complex or float): + returns scalar if both wavelen and radius are scalars. If wavelen is an + array, returns shape [len(wavelen), 1]. If radius is an array, returns + shape [1, len(radius)]. If both are arrays, returns shape + [len(wavelen), len(radius)] + """ + # ensure size parameter calculation broadcasts correctly when both + # wavelength and radius are arrays + radius = np.broadcast_to(radius, (np.size(wavelen), np.size(radius))) + wavelen = np.reshape(wavelen, (np.size(wavelen), 1)) + # matrix index may be specifed as an array to account for dispersion + if isinstance(n_matrix, np.ndarray): + n_matrix = np.reshape(n_matrix, (np.size(wavelen), 1)) + sp = (2 * np.pi * n_matrix / wavelen * radius) + # must use to('dimensionless') in case the wavelength and radius are # specified in different units; pint doesn't automatically make # ratios such as 'nm'/'um' dimensionless - sp = (2 * np.pi * n_matrix / wavelen * radius) if isinstance(sp, Quantity): sp = sp.to('dimensionless').magnitude + if np.size(sp) == 1: + sp = sp.item() + return sp diff --git a/pymie/mie.py b/pymie/mie.py index 1e163cf..45b22d4 100644 --- a/pymie/mie.py +++ b/pymie/mie.py @@ -21,23 +21,28 @@ Notes ----- -Based on miescatlib.py in holopy. Also includes some functions from Jerome's -old miescat_1d.py library. +Based on miescatlib.py in HoloPy, written by Jerome Fung. Also includes some +functions from Jerome's old miescat_1d.py library and Jerome's multilayer +scattering code, copied from HoloPy on 12 Sept 2017. Numerical stability not guaranteed for large nstop, so be careful when calculating very large size parameters. A better-tested (and faster) version of this code is in the HoloPy package (http://manoharan.seas.harvard.edu/holopy). +Key reference for multilayer algorithm is [3]_ + References ---------- -[1] Bohren, C. F. and Huffman, D. R. ""Absorption and Scattering of Light by +[1] Bohren, C. F. and Huffman, D. R. "Absorption and Scattering of Light by Small Particles" (1983) [2] Wiscombe, W. J. "Improved Mie Scattering Algorithms" Applied Optics 19, no. 9 (1980): 1505. doi:10.1364/AO.19.001505 +[3] Yang, "Improved recursive algorithm for light scattering by a multilayered +sphere," Applied Optics 42, 1710-1720 (2003). -.. moduleauthor :: Jerome Fung -.. moduleauthor :: Vinothan N. Manoharan -.. moduleauthor :: Sofia Magkiriadou +.. moduleauthor:: Jerome Fung +.. moduleauthor:: Vinothan N. Manoharan +.. moduleauthor:: Sofia Magkiriadou """ import warnings @@ -47,7 +52,6 @@ from scipy.integrate import trapezoid from . import Quantity, index_ratio, mie_specfuncs -from . import multilayer_sphere_lib as msl from . import size_parameter, ureg from .mie_specfuncs import DEFAULT_EPS1, DEFAULT_EPS2 # default tolerances @@ -63,15 +67,17 @@ def calc_ang_dist(m, x, angles, mie = True, check = False): Parameters ---------- - m: complex particle relative refractive index, n_part/n_med - x: size parameter, x = ka = 2*pi*n_med/lambda * a (sphere radius a) - angles: ndarray(structcol.Quantity [dimensionless]) + m : complex or float, array-like + complex particle relative refractive index, n_part/n_med + x : complex or float, array-like + size parameter, x = ka = 2*pi*n_med/lambda * a (sphere radius a) + angles : ndarray(structcol.Quantity [dimensionless]) array of angles. Must be entered as a Quantity to allow specifying units (degrees or radians) explicitly - mie: Boolean (optional) + mie : Boolean (optional) if true (default) does full Mie calculation; if false, uses RG approximation - check: Boolean (optional) + check : Boolean (optional) if true, outputs scattering efficiencies Returns @@ -89,20 +95,11 @@ def calc_ang_dist(m, x, angles, mie = True, check = False): if isinstance(x, Quantity): x = x.to('').magnitude - #initialize arrays for holding ipar and iperp - ipar = np.array([]) - iperp = np.array([]) - if mie: # Mie scattering preliminaries nstop = _nstop(np.array(x).max()) - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - coeffs = msl.scatcoeffs_multi(m, x) - else: - coeffs = _scatcoeffs(m, x, nstop) + coeffs = _scatcoeffs(m, x, nstop) n = np.arange(nstop)+1. prefactor = (2*n+1.)/(n*(n+1.)) @@ -168,16 +165,9 @@ def calc_cross_sections(m, x, wavelen_media, eps1 = DEFAULT_EPS1, where I_0 is the incident intensity. See van de Hulst, p. 14. """ # This is adapted from mie.py in holopy - # TODO take arrays for m and x to describe a multilayer sphere and return - # multilayer scattering coefficients lmax = _nstop(np.array(x).max()) - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - albl = msl.scatcoeffs_multi(m, x, eps1=eps1, eps2=eps2) - else: - albl = _scatcoeffs(m, x, lmax, eps1=eps1, eps2=eps2) + albl = _scatcoeffs(m, x, lmax, eps1=eps1, eps2=eps2) cscat, cext, cback = tuple(np.abs(wavelen_media)**2 * c/2/np.pi for c in _cross_sections(albl[0], albl[1])) @@ -199,41 +189,37 @@ def calc_efficiencies(m, x): geometrical cross-section """ nstop = _nstop(np.array(x).max()) - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - coeffs = msl.scatcoeffs_multi(m, x) - else: - coeffs = _scatcoeffs(m, x, nstop) + coeffs = _scatcoeffs(m, x, nstop) cscat, cext, cback = _cross_sections(coeffs[0], coeffs[1]) - qscat = cscat * 2./np.abs(x)**2 - qext = cext * 2./np.abs(x)**2 - qback = cback * 1./np.abs(x)**2 + # for multilayer spheres, scale by the size parameter corresponding to + # outermost radius + x_outer = np.atleast_2d(x).max(axis=-1)[:, np.newaxis] + qscat = cscat[..., np.newaxis] * 2./np.abs(x_outer)**2 + qext = cext[..., np.newaxis] * 2./np.abs(x_outer)**2 + qback = cback[..., np.newaxis] * 1./np.abs(x_outer)**2 # in order: scattering, extinction and backscattering efficiency - return qscat, qext, qback + return qscat.squeeze(), qext.squeeze(), qback.squeeze() -def calc_g(m, x): +def calc_g(m, x, nstop=None): """ Asymmetry parameter """ - nstop = _nstop(np.array(x).max()) - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - coeffs = msl.scatcoeffs_multi(m, x) - else: - coeffs = _scatcoeffs(m, x, nstop) + if nstop is None: + nstop = _nstop(np.array(x).max()) + coeffs = _scatcoeffs(m, x, nstop) - cscat = _cross_sections(coeffs[0], coeffs[1])[0] * 2./np.array(x).max()**2 - g = ((4./(np.array(x).max()**2 * cscat)) + # for multilayer particle, need to scale by the x of the outermost layer + outer_x = np.array(x).max(axis=-1).squeeze() + cscat = _cross_sections(coeffs[0], coeffs[1])[0] * 2./outer_x**2 + g = ((4./(outer_x**2 * cscat)) * _asymmetry_parameter(coeffs[0], coeffs[1])) return g -@ureg.check(None, None, '[length]', ('[]','[]', '[]')) -def calc_integrated_cross_section(m, x, wavelen_media, theta_range): +@ureg.check(None, None, '[length]', '[]') +def calc_integrated_cross_section(m, x, wavelen_media, thetas): """ Calculate (dimensional) integrated cross section using quadrature @@ -243,22 +229,19 @@ def calc_integrated_cross_section(m, x, wavelen_media, theta_range): x: size parameter wavelen_media: structcol.Quantity [length] wavelength of incident light *in media* - theta_range: tuple of structcol.Quantity [dimensionless] - first two elements specify the range of polar angles over which to - integrate the scattering. Last element specifies the number of angles. + thetas: array of structcol.Quantity [dimensionless] + polar angles over which to integrate the scattering Returns ------- cross_section : float Dimensional integrated cross-section """ - theta_min = theta_range[0].to('rad').magnitude - theta_max = theta_range[1].to('rad').magnitude - angles = Quantity(np.linspace(theta_min, theta_max, theta_range[2]), 'rad') + angles = thetas.to('rad') form_factor = calc_ang_dist(m, x, angles) - integrand_par = form_factor[0]*np.sin(angles) - integrand_perp = form_factor[1]*np.sin(angles) + integrand_par = (form_factor[0]*np.sin(angles)).magnitude + integrand_perp = (form_factor[1]*np.sin(angles)).magnitude # scipy.integrate.trapezoid does not yet preserve units, so we will remove # the units before calling and put them back afterward. Can simplify code @@ -371,39 +354,39 @@ def calc_dwell_time(radius, n_medium, n_particle, wavelen, return dwell_time +# TODO: document and test the correctness of this function, or delete (it is +# not currently used anywhere) +@ureg.check('[length]', '[]', '[]', '[length]', None, None) def calc_reflectance(radius, n_medium, n_particle, wavelen, - min_angle=np.pi/2, num_angles=50, - eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): + min_angle=np.pi/2, num_angles=50): m = index_ratio(n_particle, n_medium) x = size_parameter(wavelen, n_medium, radius) wavelen_media = wavelen/n_medium - geometric_cross_sec = np.pi*radius**2 + # geometric cross section is defined by outermost radius + if np.isscalar(radius): + rmax = radius + else: + rmax = radius.max() + geometric_cross_sec = np.pi*rmax**2 thetas = Quantity(np.linspace(min_angle, np.pi, num_angles), 'rad') - # calculate reflectance cross section - if np.imag(x)>0: - angles = Quantity(np.linspace(min_angle, np.pi, num_angles), 'rad') - distance = radius.max() - k = 2*np.pi/wavelen_media + if np.any(np.imag(x) > 0): + distance = rmax + k = np.atleast_1d(2*np.pi/wavelen_media) (diff_cscat_par, - diff_cscat_perp) = diff_scat_intensity_complex_medium(m, - x, thetas, - k*distance) - + diff_cscat_perp) = diff_scat_intensity_complex_medium(m, x, thetas, + k*distance) refl_cscat = integrate_intensity_complex_medium(diff_cscat_par, - diff_cscat_perp, - distance, - angles, k)[0] + diff_cscat_perp, + distance, thetas, k)[0] else: - refl_cscat = calc_integrated_cross_section(m, x, wavelen_media, - (thetas[0], thetas[-1], - num_angles)) + thetas) - reflectance = refl_cscat/geometric_cross_sec/wavelen_media.magnitude**2 - reflectance = reflectance.magnitude + reflectance = ((refl_cscat/geometric_cross_sec).to('') + / wavelen_media**2).squeeze() return reflectance @@ -468,21 +451,151 @@ def _pis_and_taus(nstop, thetas): return pis[...,1:nstop+1], taus[...,1:nstop+1] def _scatcoeffs(m, x, nstop, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): + # index ratio should be specified as a 2D array with shape + # [num_values, num_layers] to calculate over a set of different values, + # such as wavelengths. If specified as a 1D array, shape is [num_layers]. + if np.atleast_2d(m).shape[-1] > 1: + return _scatcoeffs_multi(m, x) + + # Scattering coefficients for single-layer particles. # see B/H eqn 4.88 # implement criterion used by BHMIE plus a couple more orders to be safe # nmx = np.array([nstop, np.round(np.absolute(m*x))]).max() + 20 # Dnmx = mie_specfuncs.log_der_1(m*x, nmx, nstop) # above replaced with Lentz algorithm - Dnmx = mie_specfuncs.dn_1_down(m * x, nstop + 1, nstop, - mie_specfuncs.lentz_dn1(m * x, nstop + 1, + z = np.atleast_1d(m * x).squeeze() + Dnmx = mie_specfuncs.dn_1_down(z, nstop + 1, nstop, + mie_specfuncs.lentz_dn1(z, nstop + 1, eps1, eps2)) n = np.arange(nstop+1) + x = np.atleast_2d(x) psi, xi = mie_specfuncs.riccati_psi_xi(x, nstop) - psishift = np.concatenate((np.zeros(1), psi))[0:nstop+1] - xishift = np.concatenate((np.zeros(1), xi))[0:nstop+1] + + # insert zeroes at the beginning of second axis (order axis) + psishift = np.pad(psi, ((0,), (1,)))[:, 0:nstop+1] + xishift = np.pad(xi, ((0,), (1,)))[:, 0:nstop+1] an = ( (Dnmx/m + n/x)*psi - psishift ) / ( (Dnmx/m + n/x)*xi - xishift ) bn = ( (Dnmx*m + n/x)*psi - psishift ) / ( (Dnmx*m + n/x)*xi - xishift ) - return np.array([an[1:nstop+1], bn[1:nstop+1]]) # output begins at n=1 + + # coefficient array has shape [2, num_values, nstop] or [2, nstop] if + # only one value (only one wavelength, for example) + return np.array([an[..., 1:nstop+1], bn[..., 1:nstop+1]]).squeeze() + +def _scatcoeffs_multi(marray, xarray, nstop=None, eps1 = 1e-3, eps2 = 1e-16): + '''Calculate scattered field expansion coefficients (in the Mie formalism) + for a particle with an arbitrary number of spherically symmetric layers + with different refractive indices. + + Parameters + ---------- + marray : array_like, complex128 + array of layer indexes, innermost first. If specified as a 2D array, + axis=0 corresponds to the values over which to vectorize (e.g. + wavelengths), and axis=1 corresponds to the layer indexes + xarray : array_like, real + array of layer size parameters (k * outer radius), innermost first. If + specified as a 2D array, axis=0 corresponds to the values over which to + vectorize and axis=1 corresponds to the layer size parameters. + nstop : int + maximum order. If not specified, uses largest value of x to determine + eps1 : float, optional + underflow criterion for Lentz continued fraction for Dn1 + eps2 : float, optional + convergence criterion for Lentz continued fraction for Dn1 + + Returns + ------- + scat_coeffs : ndarray (complex) + Scattering coefficients + + ''' + # ensure correct data types and shapes + marray = np.atleast_2d(marray).astype(complex) + xarray = np.atleast_2d(xarray).astype(complex) + + # sanity check: marray and xarray must be same size + if marray.size != xarray.size: + raise ValueError("Arrays of layer indices and size parameters must " + "have the same dimensions") + + # need number of layers L + nlayers = marray.shape[-1] + + # calculate nstop based on largest radius + if nstop is None: + nstop = _nstop(xarray.max()) + + # initialize H_n^a and H_n^b in the core, see eqns. 12a and 13a + intl = mie_specfuncs.log_der_13(marray[..., 0] * xarray[..., 0], + nstop, eps1, eps2)[0] + hans = intl + hbns = intl + + # m_l x_{l-1} + z1 = marray[..., 1:] * xarray[..., :-1] + # # m_l x_l + z2 = marray[..., 1:] * xarray[..., 1:] + + # pre-calculate logarithmic derivatives for all layers + derz1s = mie_specfuncs.log_der_13(z1, nstop, eps1, eps2) + derz2s = mie_specfuncs.log_der_13(z2, nstop, eps1, eps2) + + # pre-calculate ratio Q_n^l for all layers + Qnl_arr = mie_specfuncs.Qratio(z1, z2, nstop, dns1 = derz1s, + dns2 = derz2s, eps1 = eps1, eps2 = eps2) + + # lay is l-1 (index on layers used by Yang) + for lay in np.arange(1, nlayers): + m = marray[..., lay] + mm1 = marray[..., lay-1] + + # calculate logarithmic derivatives D_n^1 and D_n^3 + dz1s0, dz1s1 = derz1s[0][:, lay-1], derz1s[1][:, lay-1] + dz2s0, dz2s1 = derz2s[0][:, lay-1], derz2s[1][:, lay-1] + + # calculate G1, G2, Gtilde1, Gtilde2 according to + # eqns 26-29 + # using H^a_n and H^b_n from previous layer + G1 = m[:, np.newaxis]*hans - mm1[:, np.newaxis]*dz1s0 + G2 = m[:, np.newaxis]*hans - mm1[:, np.newaxis]*dz1s1 + Gt1 = mm1[:, np.newaxis]*hbns - m[:, np.newaxis]*dz1s0 + Gt2 = mm1[:, np.newaxis]*hbns - m[:, np.newaxis]*dz1s1 + + # calculate ratio Q_n^l for this layer + Qnl = Qnl_arr[:, lay-1] + + # now calculate H^a_n and H^b_n in current layer + # see eqns 24 and 25 + hans = (G2*dz2s0 - Qnl*G1*dz2s1) / (G2 - Qnl*G1) + hbns = (Gt2*dz2s0 - Qnl*Gt1*dz2s1) / (Gt2 - Qnl*Gt1) + # repeat for next layer + + # Relate H^a and H^b in the outer layer to the Mie scat coeffs + # see Yang eqns 14 and 15 + # + # n = 0 to nstop + # (below we vectorize over the first dimension of xarray; we calculate the + # max x over layers for each value of the first dimension) + psiandxi = mie_specfuncs.riccati_psi_xi(xarray.max(axis=1)[:, np.newaxis], + nstop) + n = np.arange(nstop+1) + psi = psiandxi[0] + xi = psiandxi[1] + # this doesn't bother to calculate psi/xi_{-1} correctly, + # but OK since we're throwing out a_0, b_0 where it appears + psishift = np.insert(psi, 0, + np.zeros(psi.shape[:-1]), axis=-1)[..., 0:nstop+1] + xishift = np.insert(xi, 0, + np.zeros(xi.shape[:-1]), axis=-1)[..., 0:nstop+1] + mlast = marray[..., nlayers-1][:, np.newaxis] + xlast = xarray[..., nlayers-1][:, np.newaxis] + an = (((hans/mlast + n/xlast)*psi - psishift) + / ((hans/mlast + n/xlast)*xi - xishift)) + bn = (((hbns*mlast + n/xlast)*psi- psishift) + / ((hbns*mlast + n/xlast)*xi - xishift)) + + # output begins at n=1 + return np.array([an[..., 1:nstop+1], bn[..., 1:nstop+1]]).squeeze() def _internal_coeffs(m, x, n_max, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): ''' @@ -493,14 +606,25 @@ def _internal_coeffs(m, x, n_max, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): have different conventions (labeling of c_n and d_n and factors of m) for their internal coefficients. ''' - ratio = mie_specfuncs.R_psi(x, m * x, n_max, eps1, eps2) + if np.atleast_2d(x).shape[-1] > 1: + raise ValueError("Internal Mie coefficients cannot yet be calculated " + "for layered sphere") + + m = np.reshape(m, (np.size(m), 1)).astype(complex) + x = np.reshape(x, (np.size(x), 1)).astype(complex) + z = np.array(m * x) + + ratio = mie_specfuncs.R_psi(x, z, n_max, eps1, eps2) D1x, D3x = mie_specfuncs.log_der_13(x, n_max, eps1, eps2) - D1mx = mie_specfuncs.dn_1_down(m * x, n_max + 1, n_max, - mie_specfuncs.lentz_dn1(m * x, n_max + 1, + D1mx = mie_specfuncs.dn_1_down(z, n_max + 1, n_max, + mie_specfuncs.lentz_dn1(z, n_max + 1, eps1, eps2)) - cl = m * ratio * (D3x - D1x) / (D3x - m * D1mx) - dl = m * ratio * (D3x - D1x) / (m * D3x - D1mx) - return np.array([cl[1:], dl[1:]]) # start from l = 1 + cl = (m[..., np.newaxis] * ratio * (D3x - D1x) + / (D3x - m[..., np.newaxis] * D1mx)) + dl = (m[..., np.newaxis] * ratio * (D3x - D1x) + / (m[..., np.newaxis] * D3x - D1mx)) + # start from l = 1 + return np.array([cl[..., 1:], dl[..., 1:]]).squeeze() def _trans_coeffs(m, x, n_max, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): ''' @@ -598,7 +722,7 @@ def _nstop(x): # Criterion for calculating near-field properties with exact Mie solutions # (J. R. Allardice and E. C. Le Ru, Applied Optics, Vol. 53, No. 31 (2014). - return (np.round(np.absolute(x+11*x**(1./3.)+1))).astype('int') + return (np.round(np.absolute(x+11*x**(1./3.)+1))).squeeze().astype('int') def _asymmetry_parameter(al, bl): ''' @@ -607,13 +731,14 @@ def _asymmetry_parameter(al, bl): See discussion in Bohren & Huffman p. 120. The output of this function omits the prefactor of 4/(x^2 Q_sca). ''' - lmax = al.shape[0] + # axis -1 (last axis) is order axis + lmax = al.shape[-1] l = np.arange(lmax) + 1 selfterm = (l[:-1] * (l[:-1] + 2.) / (l[:-1] + 1.) * - np.real(al[:-1] * np.conj(al[1:]) + - bl[:-1] * np.conj(bl[1:]))).sum() + np.real(al[..., :-1] * np.conj(al[..., 1:]) + + bl[..., :-1] * np.conj(bl[..., 1:]))).sum(axis=-1) crossterm = ((2. * l + 1.)/(l * (l + 1)) * - np.real(al * np.conj(bl))).sum() + np.real(al * np.conj(bl))).sum(axis=-1) return selfterm + crossterm def _cross_sections(al, bl): @@ -625,13 +750,13 @@ def _cross_sections(al, bl): The output omits a scaling prefactor of 2 * pi / k^2 = lambda_media^2/2/pi. ''' - lmax = al.shape[0] + lmax = al.shape[-1] l = np.arange(lmax) + 1 prefactor = (2. * l + 1.) - cscat = (prefactor * (np.abs(al)**2 + np.abs(bl)**2)).sum() - cext = (prefactor * np.real(al + bl)).sum() + cscat = (prefactor * (np.abs(al)**2 + np.abs(bl)**2)).sum(axis=-1) + cext = (prefactor * np.real(al + bl)).sum(axis=-1) # see p. 122 and discussion in that section. The formula on p. 122 # calculates the backscattering cross-section according to the traditional @@ -639,7 +764,7 @@ def _cross_sections(al, bl): # jettison the factor of 4*pi to get values that correspond to the # differential scattering cross-section in the backscattering direction. alts = 2. * (np.arange(lmax) % 2) - 1 - cback = (np.abs((prefactor * alts * (al - bl)).sum())**2)/4.0/np.pi + cback = (np.abs((prefactor * alts * (al - bl)).sum(axis=-1))**2)/4.0/np.pi return cscat, cext, cback @@ -667,6 +792,10 @@ def _cross_sections_complex_medium_fu(al, bl, cl, dl, radius, n_particle, in an absorbing medium". Applied Optics, 40, 9 (2001). ''' + # ensure broadcasting will work correctly + num_wavelen = np.atleast_1d(wavelen).shape[0] + wavelen = np.reshape(wavelen, (num_wavelen, 1)) + # if the imaginary part of the medium index is close to 0, then use the # limit value of prefactor1 for the calculations if n_medium.imag.magnitude <= 1e-7: @@ -676,39 +805,42 @@ def _cross_sections_complex_medium_fu(al, bl, cl, dl, radius, n_particle, prefactor1 = eta**2 * wavelen / (2*np.pi*radius**2*n_medium.real* (1+(eta-1)*np.exp(eta))) - lmax = al.shape[0] + lmax = np.atleast_1d(al).shape[-1] l = np.arange(lmax) + 1 - prefactor2 = (2. * l + 1.) + prefactor2 = (2. * l + 1.)[np.newaxis, ...] # calculate the scattering efficiency _, xi = mie_specfuncs.riccati_psi_xi(x_medium, lmax) - xishift = np.concatenate((np.zeros(1), xi))[0:lmax+1] - xi = xi[1:] - xishift = xishift[1:] + xishift = np.insert(xi, 0, + np.zeros(xi.shape[:-1]), axis=-1)[..., 0:lmax+1] + xi = xi[..., 1:] + xishift = xishift[..., 1:] Bn = (np.abs(al)**2 * (xishift - l*xi/x_medium) * np.conj(xi) - np.abs(bl)**2 * xi * np.conj(xishift - l*xi/x_medium)) / (2*np.pi*n_medium/wavelen) - Qscat = prefactor1 * np.sum(prefactor2 * Bn.imag) + Qscat = prefactor1 * np.sum(prefactor2 * Bn.imag, axis=-1)[..., np.newaxis] # calculate the absorption and extinction efficiencies psi, _ = mie_specfuncs.riccati_psi_xi(x_scatterer, lmax) - psishift = np.concatenate((np.zeros(1), psi))[0:lmax+1] - psi = psi[1:] - psishift = psishift[1:] + psishift = np.insert(psi, 0, + np.zeros(xi.shape[:-1]), axis=-1)[..., 0:lmax+1] + psi = psi[..., 1:] + psishift = psishift[..., 1:] An = (np.abs(cl)**2 * psi * np.conj(psishift - l*psi/x_scatterer) - np.abs(dl)**2 * (psishift - l*psi/x_scatterer)* np.conj(psi)) / (2*np.pi*n_particle/wavelen) - Qabs = prefactor1 * np.sum(prefactor2 * An.imag) - Qext = prefactor1 * np.sum(prefactor2 * (An+Bn).imag) + Qabs = prefactor1 * np.sum(prefactor2 * An.imag, axis=-1)[..., np.newaxis] + Qext = (prefactor1 * + np.sum(prefactor2 * (An+Bn).imag, axis=-1)[..., np.newaxis]) # calculate the cross sections Cscat = Qscat *np.pi * radius**2 Cabs = Qabs *np.pi * radius**2 Cext = Qext *np.pi * radius**2 - return(Cscat, Cabs, Cext) + return(Cscat.squeeze(), Cabs.squeeze(), Cext.squeeze()) def _cross_sections_complex_medium_sudiarta(al, bl, x, radius): ''' @@ -726,59 +858,69 @@ def _cross_sections_complex_medium_sudiarta(al, bl, x, radius): (2001). ''' + # if multilayer, use outermost radius and size parameter corresponding to + # outermost radius radius = np.array(radius.magnitude).max() * radius.units - x = np.array(x).max() + x = np.array(x).max(axis=-1)[..., np.newaxis] k = x/radius - lmax = al.shape[0] + lmax = np.atleast_1d(al).shape[-1] l = np.arange(lmax) + 1 - prefactor = (2. * l + 1.) + prefactor = (2. * l + 1.)[np.newaxis, ...] # if the imaginary part of k is close to 0 (because the medium index is # close to 0), then use the limit value of factor for the calculations - if k.imag.magnitude <= 1e-8: - factor = 1/2 - else: - exponent = np.exp(2*radius*k.imag) - factor = (exponent/(2*radius*k.imag)+(1-exponent)/(2*radius*k.imag)**2) + # (see eq 10 of Sudiarta and Chylek for I_denom; the cross-section is + # calculated from W/I_denom) + factor_limit = 1/2 + exponent = np.exp(2*radius*k.imag) + with np.errstate(divide='ignore', invalid='ignore'): + factor = np.where(k.imag <= Quantity(1e-8, '1/nm'), factor_limit, + (exponent/(2*radius*k.imag) + + (1-exponent)/(2*radius*k.imag)**2)) I_denom = k.real * factor - _, xi = mie_specfuncs.riccati_psi_xi(x, lmax) - xishift = np.concatenate((np.zeros(1), xi))[0:lmax+1] - xi = xi[1:] - xishift = xishift[1:] + psi, xi = mie_specfuncs.riccati_psi_xi(x, lmax) + + xishift = np.insert(xi, 0, + np.zeros(xi.shape[:-1]), axis=-1)[..., 0:lmax+1] + xi = xi[..., 1:] + xishift = xishift[..., 1:] xideriv = xishift - l*xi/x - psi, _ = mie_specfuncs.riccati_psi_xi(x, lmax) - psishift = np.concatenate((np.zeros(1), psi))[0:lmax+1] - psi = psi[1:] - psishift = psishift[1:] + psishift = np.insert(psi, 0, + np.zeros(xi.shape[:-1]), axis=-1)[..., 0:lmax+1] + psi = psi[..., 1:] + psishift = psishift[..., 1:] psideriv = psishift - l*psi/x - # calculate the scattering cross section + # calculate the scattering cross section from eq 5 of Sudiarta and Chylek term1 = (-1j * np.abs(al)**2 *xideriv * np.conj(xi) + 1j* np.abs(bl)**2 * xi * np.conj(xideriv)) - numer1 = (np.sum(prefactor * term1) * k).real + numer1 = (np.sum(prefactor * term1, axis=-1)[..., np.newaxis] + * np.conj(k)).real Cscat = np.pi / np.abs(k)**2 * numer1 / I_denom - # calculate the absorption cross section + # calculate the absorption cross section from eq 7 of Sudiarta and Chylek term2 = (1j*np.conj(psi)*psideriv - 1j*psi*np.conj(psideriv) + 1j*bl*np.conj(psideriv)*xi + 1j*np.conj(bl)*psi*np.conj(xideriv) + 1j*np.abs(al)**2*xideriv*np.conj(xi) - 1j*np.abs(bl)**2*xi*np.conj(xideriv) - 1j*al*np.conj(psi)*xideriv - 1j*np.conj(al)*psideriv*np.conj(xi)) - numer2 = (np.sum(prefactor * term2) * k).real + numer2 = (np.sum(prefactor * term2, axis=-1)[..., np.newaxis] + * np.conj(k)).real Cabs = np.pi / np.abs(k)**2 * numer2 / I_denom - # calculate the extinction cross section + # calculate the extinction cross section from eq 8 of Sudiarta and Chylek term3 = (1j*np.conj(psi)*psideriv - 1j*psi*np.conj(psideriv) + 1j*bl*np.conj(psideriv)*xi + 1j*np.conj(bl)*psi*np.conj(xideriv) - 1j*al*np.conj(psi)*xideriv - 1j*np.conj(al)*psideriv*np.conj(xi)) - numer3 = (np.sum(prefactor * term3) * k).real + numer3 = (np.sum(prefactor * term3, axis=-1)[..., np.newaxis] + * np.conj(k)).real Cext = np.pi / np.abs(k)**2 * numer3 / I_denom - return(Cscat, Cabs, Cext) + return(Cscat.squeeze(), Cabs.squeeze(), Cext.squeeze()) def _scat_fields_complex_medium(m, x, thetas, kd, near_field=False): @@ -850,12 +992,7 @@ def _scat_fields_complex_medium(m, x, thetas, kd, near_field=False): nstop = _nstop(np.array(x).max()) n = np.arange(nstop)+1. - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - an, bn = msl.scatcoeffs_multi(m, x) - else: - an, bn = _scatcoeffs(m, x, nstop) + an, bn = _scatcoeffs(m, x, nstop) # calculate prefactor (omitting the incident electric field because it # cancels out when calculating the scattered intensity) @@ -869,8 +1006,8 @@ def _scat_fields_complex_medium(m, x, thetas, kd, near_field=False): # integrating to get the scattering cross section) # required for calculations with polarized light - th_shape = list(thetas.shape) - th_shape.append(len(n)) + # reshape to (num_values, num_angles, order) + th_shape = (kd.shape[0],) + thetas.shape + (len(n),) En = np.broadcast_to(En, th_shape) an = np.broadcast_to(an, th_shape) @@ -889,12 +1026,13 @@ def _scat_fields_complex_medium(m, x, thetas, kd, near_field=False): jn = spherical_jn(nstop_array, kd) yn = spherical_yn(nstop_array, kd) zn = jn + 1j*yn - zn = zn[1:] + zn = zn[..., 1:] _, xi = mie_specfuncs.riccati_psi_xi(kd, nstop) - xishift = np.concatenate((np.zeros(1), xi))[0:nstop+1] - xi = xi[1:] - xishift = xishift[1:] + # insert zeroes at the beginning of second axis (order axis) + xishift = np.pad(xi, ((0,), (1,)))[:, 0:nstop+1] + xi = xi[..., 1:] + xishift = xishift[..., 1:] bessel_deriv = xishift - n*xi/kd zn = np.broadcast_to(zn, th_shape) bessel_deriv = np.broadcast_to(bessel_deriv, th_shape) @@ -937,7 +1075,7 @@ def _scat_fields_complex_medium(m, x, thetas, kd, near_field=False): def diff_scat_intensity_complex_medium(m, x, thetas, kd, coordinate_system = 'scattering plane', phis = None, near_field=False, incident_vector=None): - ''' + """ Calculates the differential scattered intensity in an absorbing medium. User can choose whether to include near fields. @@ -970,27 +1108,33 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, Parameters ---------- - m: complex relative refractive index - x: size parameter using the medium's refractive index - thetas: array of scattering angles (Quantity in rad) - kd: k * distance, where k = 2*np.pi*n_matrix/wavelen, and distance is the + m : complex, array-like + complex particle relative refractive index, n_part/n_med + x : complex, array-like + size parameter, x = ka = 2*pi*n_med/lambda * a (sphere radius a) + thetas : array-like (Quantity [angle]) + scattering angles. Should be 2D, as output from np.meshgrid, if + coordinate_system = 'cartesian' + kd : float (Quantity, dimensionless) + k * distance, where k = 2*np.pi*n_matrix/wavelen, and distance is the distance away from the center of the particle. The standard far-field solutions are obtained when distance >> radius in a non-absorbing - medium. (Quantity, dimensionless) - coordinate_system: string + medium. + coordinate_system : string default value 'scattering plane' means scattering calculations will be carried out in the basis defined by basis vectors parallel and perpendicular to scattering plane. Variable also accepts value 'cartesian' which scattering calculations will be carried out in the basis defined by basis vectors x and y in the lab frame, with z as the direction of propagation. - phis: None or ndarray + phis : None or ndarray azimuthal angles for which to calculate the diff scat intensity. In the 'scattering plane' coordinate system, the scattering matrix does not depend on phi, so phi should be set to None. In the 'cartesian' coordinate system, the scattering matrix does depend on phi, so an - array of values should be provided. - near_field: boolean + array of values should be provided. For 'cartesian' both thetas and + phis should be 2D, as output from np.meshgrid. + near_field : boolean True to include the near-fields (default is False). Cannot be set to True while using coordinate_system='cartesian' because near field solutions are not implemented for cartesian coordinate system. Also @@ -1004,7 +1148,7 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, near fields and still integrate at the surface of the particle, we use the asymptotic form of the spherical Hankel function in the far field (p. 94 of Bohren and Huffman). - incident_vector: None or tuple + incident_vector : None or tuple vector describing the incident electric field. It is multiplied by the amplitude scattering matrix to find the vector scattering amplitude. If coordinate_system is 'scattering plane', then this vector should be in @@ -1026,16 +1170,18 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, Returns ------- - I components: tuple + I components : tuple of arrays tuple of the two orthogonal components of scattered intensity. If in cartesian coordinate system, each component is a function of theta and phi values. If in scattering plane coordinate system, each component - is an array of theta values (dimensionless).These intensities are + is an array of theta values (dimensionless). These intensities are technically 'unitless.' The intensities would get their units from the E_n term in the fields, which gets its units from an E_0 term, which is taken to be 1 here. To get an intensity with real units you would need to multiply these by |E_0|**2 where E_0 is the amplitude - of the incident wave at the origin. + of the incident wave at the origin. If multiple wavelengths are + specified, the shape of each array is (num_values, num_theta, + [num_phi]). Otherwise just (num_theta, [num_phi]). References ---------- @@ -1044,9 +1190,15 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, Q. Fu and W. Sun, "Mie theory for light scattering by a spherical particle in an absorbing medium". Applied Optics, 40, 9 (2001). - ''' + """ if isinstance(kd, Quantity): - kd = kd.to('') + kd = kd.to('').magnitude + + # ensure that broadcasting will work correctly + kd = np.atleast_1d(kd)[:, np.newaxis] + if coordinate_system == 'cartesian': + # add another axis to correspond to phi + kd = kd[:, np.newaxis] if near_field: if coordinate_system == 'scattering plane': @@ -1080,17 +1232,18 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, # get the final intensity in a non-absorbing medium (p. 113 of Bohren # and Huffman). factor = np.exp(-2*kd.imag) / ((kd.real)**2 + (kd.imag)**2) - I_1 = (np.abs(vec_scat_amp_1)**2)*factor.to('') # par or x - I_2 = (np.abs(vec_scat_amp_2)**2)*factor.to('') # perp or y + I_1 = (np.abs(vec_scat_amp_1)**2)*factor # par or x + I_2 = (np.abs(vec_scat_amp_2)**2)*factor # perp or y - return I_1.real, I_2.real # the intensities should be real + # the intensities should be real + return I_1.real.squeeze(), I_2.real.squeeze() def integrate_intensity_complex_medium(I_1, I_2, distance, thetas, k, phi_min=Quantity(0.0, 'rad'), phi_max=Quantity(2*np.pi, 'rad'), coordinate_system = 'scattering plane', phis = None): - ''' + """ Calculates the scattering cross section by integrating the differential scattered intensity at a distance of our choice in an absorbing medium. Choosing the right distance is essential in an absorbing medium because the @@ -1100,43 +1253,74 @@ def integrate_intensity_complex_medium(I_1, I_2, distance, thetas, k, Parameters ---------- - I_1, I_2: nd arrays + I_1, I_2 : array-like with shape (num_values, num_thetas, [num_phis]) differential scattered intensities, can be functions of theta or of theta and phi. If a function of theta and phi, the theta dimension MUST come first - distance: float (Quantity in [length]) + distance : float (Quantity in [length]) distance away from the scatterer - thetas: nd array (Quantity in rad) + thetas : array-like of Quantity in [angle], shape num_thetas scattering angles - k: wavevector given by 2 * pi * n_medium / wavelength - (Quantity in [1/length]) - phi_min: float (Quantity in rad). + k : array-like of Quantity in [1/length], shape num_values + wavevector given by 2 * pi * n_medium / wavelength + phi_min : float (Quantity in [angle]). minimum azimuthal angle, default set to 0 optional, only necessary if coordinate_system is 'scattering plane' - phi_max: float (Quantity in rad). + phi_max : float (Quantity in [angle]). maximum azimuthal angle, default set to 2pi optional, only necessary if coordinate_system is 'scattering plane' - phis: None or ndarray + phis : None or ndarray (Quantity in [angle], shape num_phis) azimuthal angles Returns ------- - sigma: float (in units of length**2) + sigma: array-like of Quantity in length**2 with shape num_values integrated cross section - sigma_1: float (in units of length**2) + sigma_1: array-like of Quantity in length**2 with shape num_values integrated cross section for first component of basis - sigma_2: float (in units of length**2) + sigma_2: array-like of Quantity in length**2 with shape num_values integrated cross section for second component of basis - dsigma_1: ndarray (in units of length**2) + dsigma_1: array-like (Quantity in length**2), shape num_values, num_angles differential cross section for first component of basis - dsigma_2: ndarray (in units of length**2) + dsigma_2: array-like (Quantity in length**2), shape num_values, num_angles differential cross section for second component of basis - ''' + """ # convert to radians from whatever units the user specifies if isinstance(thetas, Quantity): thetas = thetas.to('rad').magnitude + # check that if phis is specified, both thetas and phis are given as 2D + # arrays + if phis is not None: + # convert to radians + if isinstance(phis, Quantity): + phis = phis.to('rad').magnitude + if phis.ndim == 1: + phis, thetas = np.meshgrid(phis, thetas) + + # reshape arrays for broadcasting. We use length of k (which should be + # num_values) to determine whether I_1 and I_2 are specified as + # (num_values, num_angles) or just (num_angles) + k = np.atleast_1d(k) + num_values = k.shape[0] + num_thetas = thetas.shape[0] + if phis is not None: + num_phis = phis.shape[-1] + k = k.reshape((num_values, 1, 1)) + I_1 = I_1.reshape((num_values, num_thetas, num_phis)) + I_2 = I_2.reshape((num_values, num_thetas, num_phis)) + thetas = thetas[np.newaxis, ...] + # phis has only two dimensions because by the time we use it, we have + # already integrated over theta + phis = phis[0, :] + phis = phis.reshape(1, num_phis) + else: + I_1 = I_1.reshape((num_values, num_thetas)) + I_2 = I_2.reshape((num_values, num_thetas)) + k = k.reshape((num_values, 1)) + thetas = thetas.reshape((1, num_thetas)) + # this line converts the unitless intensities to cross section # Multiply by distance (= to radius of particle in montecarlo.py) because # this is the integration factor over solid angles (see eq. 4.58 in @@ -1147,12 +1331,12 @@ def integrate_intensity_complex_medium(I_1, I_2, distance, thetas, k, dsigma_1 = I_1 * distance**2 dsigma_2 = I_2 * distance**2 - if coordinate_system == 'scattering plane': + if coordinate_system == "scattering plane": if phis is not None: - warnings.warn('''azimuthal angles specified for scattering plane - calculations. Scattering plane calculations do not - depend on azimuthal angle, so specified values will - be ignored''') + warnings.warn("azimuthal angles specified for scattering plane " + "calculations. Scattering plane calculations do not " + "depend on azimuthal angle, so specified values " + "will be ignored") # convert to radians phi_min = phi_min.to('rad').magnitude @@ -1186,31 +1370,25 @@ def integrate_intensity_complex_medium(I_1, I_2, distance, thetas, k, sigma_2 = (integral_perp * (phi_max/2 - np.sin(2*phi_max)/4 - phi_min/2 + np.sin(2*phi_min)/4)) - elif coordinate_system == 'cartesian': + elif coordinate_system == "cartesian": if phis is None: - raise ValueError('phis set to None, but azimuthal angle must be \ - specified for scattering calculations in \ - cartesian coordinate system') - - # convert to radians - if isinstance(phis, Quantity): - phis = phis.to('rad').magnitude - - # Integrate over theta and phi - thetas_bc = thetas.reshape((len(thetas),1)) # reshape for broadcasting + raise ValueError("phis set to None, but azimuthal angle must be " + "specified for scattering calculations in " + "cartesian coordinate system") # strip units from integrand if isinstance(dsigma_1, Quantity): - integrand_1 = dsigma_1.magnitude * np.abs(np.sin(thetas_bc)) + integrand_1 = dsigma_1.magnitude * np.abs(np.sin(thetas)) else: - integrand_1 = dsigma_1 * np.abs(np.sin(thetas_bc)) + integrand_1 = dsigma_1 * np.abs(np.sin(thetas)) if isinstance(dsigma_2, Quantity): - integrand_2 = dsigma_2.magnitude * np.abs(np.sin(thetas_bc)) + integrand_2 = dsigma_2.magnitude * np.abs(np.sin(thetas)) else: - integrand_2 = dsigma_2 * np.abs(np.sin(thetas_bc)) + integrand_2 = dsigma_2 * np.abs(np.sin(thetas)) - sigma_1 = trapezoid(trapezoid(integrand_1, x=thetas, axis=0), x=phis) - sigma_2 = trapezoid(trapezoid(integrand_2, x=thetas, axis=0), x=phis) + # Integrate over theta and phi + sigma_1 = trapezoid(trapezoid(integrand_1, x=thetas, axis=1), x=phis) + sigma_2 = trapezoid(trapezoid(integrand_2, x=thetas, axis=1), x=phis) # restore units to integral if isinstance(dsigma_1, Quantity): @@ -1218,25 +1396,33 @@ def integrate_intensity_complex_medium(I_1, I_2, distance, thetas, k, if isinstance(dsigma_2, Quantity): sigma_2 = Quantity(sigma_2, dsigma_2.units) else: - raise ValueError('The coordinate system specified has not yet been \ - implemented. Change to \'cartesian\' or \'scattering plane\'') + raise ValueError("The coordinate system specified has not yet been " + "implemented. Change to \'cartesian\' or " + "\'scattering plane\'") # multiply by factor that accounts for attenuation in the incident light # (see Sudiarta and Chylek (2001), eq 10). # if the imaginary part of k is close to 0 (because the medium index is # close to 0), then use the limit value of factor for the calculations - if k.imag <= Quantity(1e-8, '1/nm'): - factor = 2 - else: - exponent = np.exp(2*distance*k.imag) - factor = 1 / (exponent / (2*distance*k.imag)+ - (1 - exponent) / (2*distance*k.imag)**2) + exponent = np.exp(2*distance*k.imag) + factor_limit = 2 + # ignore division by zero in case k.imag=0; we'll replace the nans with the + # limit value of factor anyway + with np.errstate(divide='ignore', invalid='ignore'): + factor = np.where(k.imag <= Quantity(1e-8, '1/nm'), factor_limit, + 1 / (exponent / (2*distance*k.imag) + + (1 - exponent) / (2*distance*k.imag)**2)) + + # prepare for broadcasting (this will add trailing axes of size 1) + sigma_1 = sigma_1.reshape(factor.shape) + sigma_2 = sigma_2.reshape(factor.shape) # calculate the averaged sigma sigma = (sigma_1 + sigma_2)/2 * factor - return(sigma, sigma_1*factor, sigma_2*factor, dsigma_1*factor/2, - dsigma_2*factor/2) + return(sigma.squeeze(), (sigma_1*factor).squeeze(), + (sigma_2*factor).squeeze(), (dsigma_1*factor/2).squeeze(), + (dsigma_2*factor/2).squeeze()) def diff_abs_intensity_complex_medium(m, x, thetas, ktd): ''' @@ -1254,17 +1440,21 @@ def diff_abs_intensity_complex_medium(m, x, thetas, ktd): Parameters ---------- - m: complex relative refractive index - x: size parameter using the medium's refractive index - thetas: array of scattering angles (Quantity in rad) - ktd: kt * distance, where kt = 2*np.pi*n_particle/wavelen, and distance is + m : array-like + complex relative refractive index + x : array-like + size parameter using the medium's refractive index + thetas: Quantity[angle], array-like + array of scattering angles (Quantity in rad or degrees) + ktd: Quantity[dimensionless], array-like + kt * distance, where kt = 2*np.pi*n_particle/wavelen, and distance is the distance away from the center of the particle. The far-field - solution is obtained when distance >> radius. (Quantity, dimensionless) + solution is obtained when distance >> radius. Returns ------- - I_par, I_perp: differential absorption intensities for an array of theta - (dimensionless). + I_par, I_perp : Quantity[dimensionless], array-like + differential absorption intensities for an array of theta Reference --------- @@ -1273,8 +1463,10 @@ def diff_abs_intensity_complex_medium(m, x, thetas, ktd): ''' # convert units from whatever units the user specifies - thetas = thetas.to('rad').magnitude - ktd = ktd.to('').magnitude + if isinstance(thetas, Quantity): + thetas = thetas.to('rad').magnitude + if isinstance(ktd, Quantity): + ktd = ktd.to('').magnitude # calculate mie coefficients nstop = _nstop(np.array(x).max()) @@ -1370,7 +1562,7 @@ def amplitude_scattering_matrix(m, x, thetas, particle x: float, or array size parameter, array if multilayer particle - thetas: nd array + thetas: array theta angles coordinate_system: string default value 'scattering plane' means scattering calculations will be @@ -1379,7 +1571,7 @@ def amplitude_scattering_matrix(m, x, thetas, 'cartesian' which scattering calculations will be carried out in the basis defined by basis vectors x and y in the lab frame, with z as the direction of propagation. - phis: None or ndarray + phis: None or array azimuthal angles for which to calculate the scattering matrix. In the 'scattering plane' coordinate system, the scattering matrix does not depend on phi, so phi should be set to None. In the 'cartesian' @@ -1388,29 +1580,29 @@ def amplitude_scattering_matrix(m, x, thetas, Returns: -------- - S1, S2, S3, S4: tuple of nd arrays - amplitude scattering matrix elements for all theta. S2 and S1 have the - same shape as theta. + S1, S2, S3, S4: tuple of arrays + amplitude scattering matrix elements for all values (e.g. wavelengths) + and theta. Shapes of all arrays are (num_values, num_theta, [num_phi]) """ # calculate n-array nstop = _nstop(np.array(x).max()) n = np.arange(nstop)+1. prefactor = (2*n+1)/(n*(n+1)) + if isinstance(thetas, Quantity): + thetas = thetas.to('rad').magnitude + if isinstance(phis, Quantity): + phis = phis.to('rad').magnitude + # calculate mie coefficients - # if the index ratio m is an array with more than 1 element, it's a - # multilayer particle - if len(np.atleast_1d(m)) > 1: - coeffs = msl.scatcoeffs_multi(m, x) - else: - coeffs = _scatcoeffs(m, x, nstop) + coeffs = _scatcoeffs(m, x, nstop) # calculate amplitude scattering matrix in 'scattering plane' coordinate # system S2_sp, S1_sp = _amplitude_scattering_matrix(nstop, prefactor, coeffs, thetas) - S3_sp = 0 - S4_sp = 0 + S3_sp = np.zeros_like(S1_sp) + S4_sp = np.zeros_like(S1_sp) if coordinate_system == 'cartesian': # raise error if no phis are specified @@ -1428,7 +1620,6 @@ def amplitude_scattering_matrix(m, x, thetas, S2_xy = S2_sp*(cosphi)**2 + S1_sp*(sinphi)**2 S3_xy = S2_sp*sinphi*cosphi - S1_sp*sinphi*cosphi S4_xy = S2_sp*cosphi*sinphi - S1_sp*cosphi*sinphi - return S1_xy, S2_xy, S3_xy, S4_xy elif coordinate_system == 'scattering plane': if phis is not None: @@ -1485,11 +1676,11 @@ def vector_scattering_amplitude(m, x, thetas, incident_vector = None, Parameters: ---------- - m: float + m: float or array-like index ratio between the particle and sample - x: float + x: float or array-like size parameter - thetas: nd array + thetas: array-like scattering angles incident_vector: None or tuple vector describing the incident electric field. It is multiplied by the @@ -1513,13 +1704,13 @@ def vector_scattering_amplitude(m, x, thetas, incident_vector = None, coordinate_system: string describes the coordinate system. Can be either 'scattering plane' or 'cartesian' - phis: nd array or None + phis: ndarray or None azimuthal angles Returns: -------- - vector scattering amplitude: tuple + vector scattering amplitude: tuple of arrays (num_values, num_angles) tuple describing the vector scattering amplitude in the specified coordinate system. not normalized. ''' @@ -1547,21 +1738,47 @@ def vector_scattering_amplitude(m, x, thetas, incident_vector = None, def _amplitude_scattering_matrix(n_stop, prefactor, coeffs, thetas): - # amplitude scattering matrix from Mie coefficients + """Amplitude scattering matrix from Mie coefficients + + """ pis, taus = _pis_and_taus(n_stop, thetas) + + # to broadcast correctly over the dimensions of coeffs (which may be + # wavelength or other variable), we need to add leading dimensions to the + # pis and taus, which have shape [num_angles, order]. Result should have + # pis, taus shape: [1, ..., 1, num_angles, order] + # Similarly, we need to insert dimensions in coeffs corresponding to the + # angles in pis and taus. Result should have + # coeffs[0].shape: [num_values, ..., 1, order] + num_leading_dims = len(coeffs[0].shape[:-1]) + num_insert_dims = len(pis.shape[:-1]) + new_coeffs_shape = (coeffs.shape[:-1] + num_insert_dims*(1,) + + (coeffs.shape[-1],)) + pis = pis.reshape(num_leading_dims*(1,) + pis.shape) + taus = taus.reshape(num_leading_dims*(1,) + taus.shape) + coeffs = coeffs.reshape(new_coeffs_shape) + + # result should have shape [num_values, ..., num_angles] S1 = np.sum(prefactor*(coeffs[0]*pis + coeffs[1]*taus), axis=-1) S2 = np.sum(prefactor*(coeffs[0]*taus + coeffs[1]*pis), axis=-1) return S2, S1 def _amplitude_scattering_matrix_RG(prefactor, x, thetas): - # amplitude scattering matrix from Rayleigh-Gans approximation - u = 2 * x * np.sin(thetas/2.) - S1 = prefactor * (3./u**3) * (np.sin(u) - u*np.cos(u)) - S2 = prefactor * (3./u**3) * (np.sin(u) - u*np.cos(u)) * np.cos(thetas) - return S2, S1 + """Amplitude scattering matrix from Rayleigh-Gans approximation + + """ + if np.atleast_2d(x).shape[-1] > 1: + raise ValueError("Rayleigh-Gans approximation cannot be used for " + "layered spheres") + u = 2 * x * np.sin(thetas/2.) + # for theta=0 the limit is 1*prefactor; the following will avoid a divide + # by zero error by dividing everywhere that u!=0 and returning 1 where u=0 + p = np.divide(np.sin(u) - u*np.cos(u), u**3, out=np.ones_like(u), + where=u!=0) -# TODO: copy multilayer code from multilayer_sphere_lib.py in holopy and -# integrate with functions for calculating scattering cross sections and form -# factor. + # result should have shape [num_values, ..., num_angles] + S1 = prefactor * 3 * p + S2 = S1 * np.cos(thetas) + return S2, S1 diff --git a/pymie/mie_specfuncs.py b/pymie/mie_specfuncs.py index 8071dfb..2007e03 100644 --- a/pymie/mie_specfuncs.py +++ b/pymie/mie_specfuncs.py @@ -42,8 +42,7 @@ .. moduleauthor:: Vinothan N. Manoharan """ import numpy as np -from numpy import arange, array, exp, imag, real, sin, zeros -from scipy.special import riccati_jn, riccati_yn, spherical_jn, spherical_yn +from scipy.special import spherical_jn, spherical_yn # default tolerances DEFAULT_EPS1 = 1e-3 @@ -52,66 +51,95 @@ def riccati_psi_xi(x, nstop): """ Construct riccati hankel function of 1st kind by linear combination of - RB's based on j_n and y_n + RB's based on j_n and y_n. + + Notes: we use scipy's spherical_jn/spherical_yn, which are vectorized, + rather than riccati_jn/riccati_yn, which are not. Also we use the general + formula which is valid for both real and complex x """ - if np.imag(x) != 0.: - # if x is complex, calculate spherical bessel functions and compute the - # complex riccati-bessel solutions - nstop_array = np.arange(0,nstop+1) - psin = spherical_jn(nstop_array, x)*x - xin = psin + 1j*spherical_yn(nstop_array, x)*x - rbh = array([psin, xin]) - else: - x = x.real - psin = riccati_jn(nstop, x) - # scipy sign on y_n consistent with B/H - xin = psin[0] + 1j*riccati_yn(nstop, x)[0] - rbh = array([psin[0], xin]) + + # Calculate spherical bessel functions and compute the complex + # riccati-bessel solutions + nstop_array = np.arange(0,nstop+1) + psin = spherical_jn(nstop_array, x)*x + xin = psin + 1j*spherical_yn(nstop_array, x)*x + rbh = np.array([psin, xin]) return rbh def lentz_dn1(z, n, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): - """ - Calculate logarithmic derivative D_n(z) of the Riccati-Bessel - function for a single value of n using the Lentz (1976) - continued fraction method. + """Calculate logarithmic derivative D_n(z) of the Riccati-Bessel function + for a single value of n using the Lentz (1976) continued fraction method + [1]. + + This function is used to seed the down-recursion algorithm in + `dn_1_down()`, which calculates the other logarithmic derivatives needed to + calculate the Mie coefficients. Notes ----- Implements check/workaround for ill-conditioning described under "Algorithm - Improvement" in Lentz (1976); see also Wiscombe/NCAR Mie report. + Improvement" in Lentz (1976) [1]_; see also Wiscombe/NCAR Mie report [2]_. + Function is vectorized across z (uses array operations to handle multiple + values of z in parallel) but still requires a loop to test for convergence Parameters ---------- - z: complex argument - n: order of the logarithmic derivative - eps1: value of continued fraction numerator or denominator - triggering ill-conditioning workaround. Recommend - 1e-3. - eps2: converge when additional products in continued fraction - differ by less than eps2 from 1. Recommend 1e-16. + z : array-like (complex) + argument of the logarithmic derivative + n : integer + order of the logarithmic derivative + eps1 : float + value of continued fraction numerator or denominator triggering + ill-conditioning workaround. Recommend 1e-3. + eps2 : float + converge when additional products in continued fraction differ by less + than eps2 from 1. Recommend 1e-16. Returns ------- - value of D_n(z) + array-like : shape of z + value of D_n(z) + + References + ---------- + [1] W. J. Lentz, "Generating Bessel functions in Mie scattering + calculations using continued fractions," Applied Optics 15(3), 668-671 + (1976). + + [2] w. Wiscombe, "Mie Scattering Calculations: Advances in Technique and + Fast, Vector-speed Computer Codes," University Corporation for Atmospheric + Research. https://doi.org/10.5065/D6ZP4414 (Original work published 1979) + """ def a_i(i): - return (-1.)**(i + 1) * 2. * (n + i - 0.5) / z + # elements of continued fraction expansion + # Note: Wiscombe uses a different value for a_1 than the one given by + # this formula, but Lentz uses the same formula for all a_i + val = (-1.)**(i + 1) * 2. * (n + i - 0.5) / z + return np.atleast_1d(val) + # calculate numerator and denominator and a_2 numerator = a_i(2) + 1. / a_i(1) denominator = a_i(2) - product = a_i(1) * numerator / denominator - ratio = product - + ratio = product.copy() ctr = 3 - while abs(product.real - 1) > eps2 or abs(product.imag) > eps2: + while ((np.abs(product.real - 1) > eps2).any() + or (np.abs(product.imag) > eps2).any()): + # find the indices corresponding to the values that have not converged + # yet, and iterate only those (otherwise further iteration of values + # that have converged could lead to instabilities) + ind = np.where(np.logical_or(np.abs(product.real - 1) > eps2, + np.abs(product.imag) > eps2)) ai = a_i(ctr) - numerator = ai + 1. / numerator - denominator = ai + 1. / denominator - - if abs(numerator / ai) < eps1 or abs(denominator / ai) < eps1: + numerator[ind] = ai[ind] + 1. / numerator[ind] + denominator[ind] = ai[ind] + 1. / denominator[ind] + # TODO check the ill-conditioning code and update to include indices. + # Requires unit tests to be written. + if ((np.abs(numerator / ai) < eps1).any() + or (np.abs(denominator / ai) < eps1).any()): # ill conditioning xi1 = 1. + a_i(ctr + 1) * numerator xi2 = 1. + a_i(ctr + 1) * denominator @@ -120,9 +148,10 @@ def a_i(i): numerator = a_i(ctr + 2) + numerator / xi1 denominator = a_i(ctr + 2) + denominator / xi2 ctr = ctr + 2 - product = numerator / denominator - ratio = ratio * product + product[ind] = numerator[ind] / denominator[ind] + ratio[ind] = ratio[ind] * product[ind] ctr = ctr + 1 + return ratio - n / z def dn_1_down(z, nmx, nstop, start_val): @@ -141,12 +170,18 @@ def dn_1_down(z, nmx, nstop, start_val): ----- psi_n(z) is related to the spherical Bessel function j_n(z). ''' - dn = zeros(nmx+1, dtype = 'complex128') - dn[nmx] = start_val + z = np.array(z) + start_val = np.atleast_1d(start_val) + dn = np.zeros(start_val.shape + (nmx+1,), dtype=complex) + dn[..., nmx] = start_val - for i in np.arange(nmx-1, -1, -1): - dn[i] = (i+1.)/z - 1.0/(dn[i+1] + (i+1.)/z) - return dn[0:nstop+1] + # pre-calculate 1/z so we save a little time in the loop + i_range = np.arange(nmx-1, -1, -1) + one_over_z = 1/z + for i in i_range: + dn[..., i] = ((i+1)*one_over_z + - 1.0/(dn[..., i+1] + ((i+1)*one_over_z))) + return dn[..., 0:nstop+1] def log_der_13(z, nstop, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): @@ -164,7 +199,7 @@ def log_der_13(z, nstop, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): eps1: underflow criterion for Lentz continued fraction for Dn1 eps2: convergence criterion for Lentz continued fraction for Dn1 ''' - z = np.complex128(z) # convert to double precision + z = np.atleast_1d(z).astype(complex) # Calculate Dn_1 (based on \psi(z)) using downward recursion. # See Mackowski eqn. 62 @@ -174,33 +209,53 @@ def log_der_13(z, nstop, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): # Calculate Dn_3 (based on \xi) by up recurrence # initialize - dn3 = zeros(nstop+1, dtype = 'complex128') - psixi = zeros(nstop+1, dtype = 'complex128') - dn3[0] = 1.j - psixi[0] = -1j*exp(1.j*z)*sin(z) - for dindex in arange(1, nstop+1): + dn3 = np.zeros(z.shape + (nstop+1,), dtype = complex) + psixi = np.zeros(z.shape + (nstop+1,), dtype = complex) + dn3[..., 0] = 1.j + psixi[..., 0] = -1j*np.exp(1.j*z)*np.sin(z) + for dindex in np.arange(1, nstop+1): # Mackowski eqn 63 - psixi[dindex] = psixi[dindex-1] * ( (dindex/z) - dn1[dindex-1]) * ( - (dindex/z) - dn3[dindex-1]) + psixi[..., dindex] = (psixi[..., dindex-1] + * ( (dindex/z) - dn1[..., dindex-1]) + * ( (dindex/z) - dn3[..., dindex-1])) # Mackowski eqn 64 - dn3[dindex] = dn1[dindex] + 1j/psixi[dindex] + dn3[..., dindex] = dn1[..., dindex] + 1j/psixi[..., dindex] return dn1, dn3 # calculate ratio of RB's defined in Yang eqn. 23 by up recursion relation def Qratio(z1, z2, nstop, dns1 = None, dns2 = None, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): - ''' - Calculate ratio of Riccati-Bessel functions defined in Yang eq. 23 - by up recursion. + """Calculate ratio of Riccati-Bessel functions defined in Yang eq. 23 by + up recursion. + Notes + ----- Logarithmic derivatives calculated automatically if not specified. - ''' - # convert z1 and z2 to 128 bit complex to prevent division problems - z1 = np.complex128(z1) - z2 = np.complex128(z2) - if dns1 is None: + Inputs z1 and z2 should be 2d complex arrays with shape [num_values, + num_layers], where num_values could be the number of wavelengths or other + variable. + + Parameters + ---------- + z1 : array-like with shape [num_values, num_layers] + m for layer * x for previous layer + z2 : array-like with shape [num_values, num_layers] + m for layer * x for layer + nstop : integer + maximum order of computation + eps1 : float + underflow criterion to pass to Lentz continued fraction calculation + eps2 : float + convergence criterion to pass to Lentz continued fraction calculation + + Returns + ------- + Qnl : array-like with shape [num_values, num_layers, order] + Q_n^l for all values (e.g. wavelengths) and layers in z + """ + if (dns1 is None) and (dns2 is None): logdersz1 = log_der_13(z1, nstop, eps1, eps2) logdersz2 = log_der_13(z2, nstop, eps1, eps2) d1z1 = logdersz1[0] @@ -213,19 +268,36 @@ def Qratio(z1, z2, nstop, dns1 = None, dns2 = None, d1z2 = dns2[0] d3z2 = dns2[1] - qns = zeros(nstop+1, dtype = 'complex128') - # initialize according to Yang eqn. 34 - a1 = real(z1) - a2 = real(z2) - b1 = imag(z1) - b2 = imag(z2) - qns[0] = exp(-2.*(b2-b1)) * (exp(-1j*2.*a1)-exp(-2.*b1)) / (exp(-1j*2.*a2) - - exp(-2.*b2)) - # Loop to do upwards recursion in eqn. 33 - for i in arange(1, nstop+1): - qns[i] = qns[i-1]* (((d3z1[i] + i/z1) * (d1z2[i] + i/z2)) - / ((d3z2[i] + i/z2) * (d1z1[i] + i/z1))) + a1 = np.real(z1) + a2 = np.real(z2) + b1 = np.imag(z1) + b2 = np.imag(z2) + qns0 = (np.exp(-2.*(b2-b1)) * (np.exp(-1j*2.*a1)-np.exp(-2.*b1)) + / (np.exp(-1j*2.*a2) - np.exp(-2.*b2))) + # shape is [num_values, num_layers, order] + qns0 = qns0[:, :, np.newaxis] + + # Vectorized loop (using np.cumprod) to do upwards recursion in eqn. 33 + irange = np.arange(1, nstop+1) + # shape is [num_values, num_layers, order] + i_over_z1 = irange[np.newaxis, np.newaxis, :]/z1[:, :, np.newaxis] + i_over_z2 = irange[np.newaxis, np.newaxis, :]/z2[:, :, np.newaxis] + prod = ((d3z1[..., 1:] + i_over_z1) * (d1z2[..., 1:] + i_over_z2) + / ((d3z2[..., 1:] + i_over_z2) * (d1z1[..., 1:] + i_over_z1))) + qns = np.concatenate((qns0, qns0 * np.cumprod(prod, axis=-1)), axis=-1) + + # equivalent non-vectorized loop is below + # + # qns = np.zeros(z1.shape + (nstop+1,), dtype=complex) + # qns[..., 0] = (np.exp(-2.*(b2-b1)) * (np.exp(-1j*2.*a1)-np.exp(-2.*b1)) + # / (np.exp(-1j*2.*a2) - np.exp(-2.*b2))) + # for i in np.arange(1, nstop+1): + # qns[..., i] = qns[..., i-1]* (((d3z1[..., i] + i/z1) + # * (d1z2[..., i] + i/z2)) + # / ((d3z2[..., i] + i/z2) + # * (d1z1[..., i] + i/z1))) + return qns def R_psi(z1, z2, nmax, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): @@ -233,13 +305,31 @@ def R_psi(z1, z2, nmax, eps1 = DEFAULT_EPS1, eps2 = DEFAULT_EPS2): Calculate ratio of Riccati-Bessel function psi: psi(z1)/psi(z2). See Mackowski eqns. 65-66. + + z1, z2 are complex arrays with shape [num_values, 1] ''' - output = zeros(nmax + 1, dtype = 'complex128') - output[0] = sin(z1) / sin(z2) + # Vectorized loop (using np.cumprod) to do up recursion + output_0 = (np.sin(z1) / np.sin(z2))[:, :, np.newaxis] dnz1 = dn_1_down(z1, nmax + 1, nmax, lentz_dn1(z1, nmax + 1, eps1, eps2)) dnz2 = dn_1_down(z2, nmax + 1, nmax, lentz_dn1(z2, nmax + 1, eps1, eps2)) - # use up recursion - for i in arange(1, nmax + 1): - output[i] = output[i - 1] * (dnz2[i] + i / z2) / (dnz1[i] + i / z1) - return output + irange = np.arange(1, nmax+1) + # shape is [num_values, num_layers, order] + i_over_z1 = irange[np.newaxis, np.newaxis, :]/z1[:, :, np.newaxis] + i_over_z2 = irange[np.newaxis, np.newaxis, :]/z2[:, :, np.newaxis] + prod = (dnz2[..., 1:] + i_over_z2) / (dnz1[..., 1:] + i_over_z1) + output_vec = np.concatenate((output_0, + output_0 * np.cumprod(prod, axis=-1)), + axis=-1) + + # equivalent non-vectorized loop is below + # + # output = np.zeros(z1.shape + (nmax + 1,), dtype=complex) + # output[..., 0] = np.sin(z1) / np.sin(z2) + # dnz1 = dn_1_down(z1, nmax + 1, nmax, lentz_dn1(z1, nmax + 1, eps1, eps2)) + # dnz2 = dn_1_down(z2, nmax + 1, nmax, lentz_dn1(z2, nmax + 1, eps1, eps2)) + # for i in np.arange(1, nmax + 1): + # output[..., i] = output[..., i-1] * ((dnz2[..., i] + i / z2) + # / (dnz1[..., i] + i / z1)) + + return output_vec diff --git a/pymie/multilayer_sphere_lib.py b/pymie/multilayer_sphere_lib.py deleted file mode 100755 index 989a7e7..0000000 --- a/pymie/multilayer_sphere_lib.py +++ /dev/null @@ -1,126 +0,0 @@ -# Copyright 2011-2016, Vinothan N. Manoharan, Thomas G. Dimiduk, -# Rebecca W. Perry, Jerome Fung, Ryan McGorty, Anna Wang, Solomon Barkley -# -# This file is part of the python-mie python package. -# -# This package is free software: you can redistribute it and/or modify it under -# the terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. -# -# This package is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -# details. -# -# You should have received a copy of the GNU General Public License along with -# this package. If not, see . - -''' -multilayer_sphere_lib.py - -Author: -Jerome Fung (fung@physics.harvard.edu) - -Copied from holopy on 12 Sept 2017. Functions to calculate the scattering from -a spherically symmetric particle with an arbitrary number of layers with -different refractive indices. - -Key reference for multilayer algorithm: -Yang, "Improved recursive algorithm for light scattering by a multilayered -sphere," Applied Optics 42, 1710-1720, (1993). - -''' - -import numpy as np - -try: - from . import mie - from .mie_specfuncs import Qratio, log_der_13, riccati_psi_xi -except ImportError: - pass - -def scatcoeffs_multi(marray, xarray, eps1 = 1e-3, eps2 = 1e-16): - ''' - Calculate scattered field expansion coefficients (in the Mie formalism) - for a particle with an arbitrary number of spherically symmetric layers. - - Parameters - ---------- - marray : array_like, complex128 - array of layer indices, innermost first - xarray : array_like, real - array of layer size parameters (k * outer radius), innermost first - eps1 : float, optional - underflow criterion for Lentz continued fraction for Dn1 - eps2 : float, optional - convergence criterion for Lentz continued fraction for Dn1 - - Returns - ------- - scat_coeffs : ndarray (complex) - Scattering coefficients - ''' - # ensure correct data types - marray = np.array(marray, dtype = 'complex128') - xarray = np.array(xarray, dtype = 'complex128') - - # sanity check: marray and xarray must be same size - if marray.size != xarray.size: - raise ValueError('Arrays of layer indices \ - and size parameters must be the same length!') - - # need number of layers L - nlayers = marray.size - - # calculate nstop based on outermost radius - nstop = mie._nstop(xarray.max()) - - # initialize H_n^a and H_n^b in the core, see eqns. 12a and 13a - intl = log_der_13(marray[0]*xarray[0], nstop, eps1, eps2)[0] - hans = intl - hbns = intl - - # lay is l-1 (index on layers used by Yang) - for lay in np.arange(1, nlayers): - z1 = marray[lay]*xarray[lay-1] # m_l x_{l-1} - z2 = marray[lay]*xarray[lay] # m_l x_l - - # calculate logarithmic derivatives D_n^1 and D_n^3 - derz1s = log_der_13(z1, nstop, eps1, eps2) - derz2s = log_der_13(z2, nstop, eps1, eps2) - - # calculate G1, G2, Gtilde1, Gtilde2 according to - # eqns 26-29 - # using H^a_n and H^b_n from previous layer - G1 = marray[lay]*hans - marray[lay-1]*derz1s[0] - G2 = marray[lay]*hans - marray[lay-1]*derz1s[1] - Gt1 = marray[lay-1]*hbns - marray[lay]*derz1s[0] - Gt2 = marray[lay-1]*hbns - marray[lay]*derz1s[1] - - # calculate ratio Q_n^l for this layer - Qnl = Qratio(z1, z2, nstop, dns1 = derz1s, dns2 = derz2s, eps1 = eps1, - eps2 = eps2) - - # now calculate H^a_n and H^b_n in current layer - # see eqns 24 and 25 - hans = (G2*derz2s[0] - Qnl*G1*derz2s[1]) / (G2 - Qnl*G1) - hbns = (Gt2*derz2s[0] - Qnl*Gt1*derz2s[1]) / (Gt2 - Qnl*Gt1) - # repeat for next layer - - # Relate H^a and H^b in the outer layer to the Mie scat coeffs - # see Yang eqns 14 and 15 - psiandxi = riccati_psi_xi(xarray.max(), nstop) # n = 0 to nstop - n = np.arange(nstop+1) - psi = psiandxi[0] - xi = psiandxi[1] - # this doesn't bother to calculate psi/xi_{-1} correctly, - # but OK since we're throwing out a_0, b_0 where it appears - psishift = np.concatenate((np.zeros(1), psi))[0:nstop+1] - xishift = np.concatenate((np.zeros(1), xi))[0:nstop+1] - - an = ((hans/marray[nlayers-1] + n/xarray[nlayers-1])*psi - psishift) / ( - (hans/marray[nlayers-1] + n/xarray[nlayers-1])*xi - xishift) - bn = ((hbns*marray[nlayers-1] + n/xarray[nlayers-1])*psi - psishift) / ( - (hbns*marray[nlayers-1] + n/xarray[nlayers-1])*xi - xishift) - return np.array([an[1:nstop+1], bn[1:nstop+1]]) # output begins at n=1 diff --git a/pymie/tests/test_mie.py b/pymie/tests/test_mie.py index 77f42e9..c936b52 100644 --- a/pymie/tests/test_mie.py +++ b/pymie/tests/test_mie.py @@ -122,7 +122,7 @@ def test_efficiencies(): 1.02022022710453, 0.51835427781473, 0.331000402174976]) - wavelen = Quantity('658.0 nm') + # wavelen = Quantity('658.0 nm') n_matrix = Quantity(1.00, '') n_particle = Quantity(1.59 + 1e-4 * 1.0j, '') m = index_ratio(n_particle, n_matrix) @@ -149,7 +149,7 @@ def test_efficiencies(): def test_absorbing_materials(): # test calculations for gold, which has a high imaginary refractive index - wavelen = Quantity('658.0 nm') + # wavelen = Quantity('658.0 nm') n_matrix = Quantity(1.00, '') n_particle = Quantity(0.1425812 + 3.6813284 * 1.0j, '') m = index_ratio(n_particle, n_matrix) @@ -583,7 +583,7 @@ def test_multilayer_complex_medium(): # Hankel equations (because they simplify when the fields are multiplied by # their conjugates to get the intensity) matches old result before simplifying cscat_imag_old = 6275.240019849266 - assert_almost_equal(cscat_imag_old, cscat_imag.magnitude, decimal=11) + assert_almost_equal(cscat_imag_old, cscat_imag.magnitude, decimal=9) assert_array_almost_equal(cscat_real.magnitude, cscat_imag.magnitude, decimal=3) @@ -632,8 +632,8 @@ def test_vector_scattering_amplitude_2d_theta_cartesian(): as_vec_x = S2_sp*np.cos(phis_2d)**2 + S1_sp*np.sin(phis_2d)**2 as_vec_y = S2_sp*np.cos(phis_2d)*np.sin(phis_2d) - S1_sp*np.cos(phis_2d)*np.sin(phis_2d) - assert_almost_equal(as_vec_x0.magnitude, as_vec_x.magnitude) - assert_almost_equal(as_vec_y0.magnitude, as_vec_y.magnitude) + assert_almost_equal(as_vec_x0, as_vec_x.magnitude) + assert_almost_equal(as_vec_y0, as_vec_y.magnitude) def test_diff_scat_intensity_complex_medium_cartesian(): ''' @@ -677,7 +677,7 @@ def test_diff_scat_intensity_complex_medium_cartesian(): I_par_perp_mag = np.sqrt(I_par**2 + I_perp**2) # check that the magnitudes are equal - assert_array_almost_equal(I_xy_mag.magnitude, I_par_perp_mag.magnitude, decimal=16) + assert_array_almost_equal(I_xy_mag, I_par_perp_mag, decimal=16) def test_integrate_intensity_complex_medium_cartesian(): ''' @@ -750,7 +750,7 @@ def test_value_errors(): I_x, I_y = mie.diff_scat_intensity_complex_medium(m, x, thetas_2d, kd, coordinate_system = 'weird', phis = phis_2d, near_field=True) - + with pytest.raises(ValueError): # try to calculate new I_x, I_y = mie.diff_scat_intensity_complex_medium(m, x, thetas_2d, kd, coordinate_system = 'cartesian', phis = phis_2d, @@ -767,13 +767,13 @@ def test_value_errors(): # integrate the differential scattered intensities cscat_xy = mie.integrate_intensity_complex_medium(I_x, I_y, distance, thetas, k, coordinate_system = 'cartesian')[0] - + with pytest.raises(ValueError): cscat_weird = mie.integrate_intensity_complex_medium(I_x, I_y, distance, thetas, k, coordinate_system = 'weird')[0] - + with pytest.raises(ValueError): as_vec_weird = mie.vector_scattering_amplitude(m, x, thetas_2d, coordinate_system = 'weird', phis = phis_2d) - + with pytest.raises(ValueError): as_vec_xy = mie.vector_scattering_amplitude(m, x, thetas_2d, coordinate_system = 'cartesian') diff --git a/pymie/tests/test_mie_vectorized.py b/pymie/tests/test_mie_vectorized.py new file mode 100644 index 0000000..a218f55 --- /dev/null +++ b/pymie/tests/test_mie_vectorized.py @@ -0,0 +1,966 @@ +# Copyright 2016, Vinothan N. Manoharan +# +# This file is part of the python-mie python package. +# +# This package is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This package is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this package. If not, see . +""" +Tests vectorization behavior of the mie module + +.. moduleauthor:: Vinothan N. Manoharan +""" + +from .. import Quantity, index_ratio, size_parameter, np, mie +from .. import mie_specfuncs +from numpy.testing import assert_allclose, assert_equal +import pytest + +def mx(num_wavelen, num_layer, start_wavelen=400, end_wavelen=800, + start_radius = 100, end_radius = 1000, + start_n_particle = 1.33+0j, end_n_particle = 1.59+0j, + n_matrix = 1.0, return_all=False): + """Convenience function to set up various combinations of wavelength and + layer inputs to test functions. + + """ + if num_layer > 1: + radius = Quantity(np.linspace(start_radius, end_radius, num_layer), + 'nm') + n_particle = np.linspace(start_n_particle, end_n_particle, num_layer) + n_particle = n_particle[np.newaxis, :] + else: + radius = Quantity(start_radius, 'nm') + n_particle = np.atleast_1d(start_n_particle)[np.newaxis, :] + if num_wavelen > 1: + wavelen = Quantity(np.linspace(start_wavelen, end_wavelen, num_wavelen), + 'nm') + # let index be the same at all wavelengths + n_particle = Quantity(np.ones((num_wavelen, num_layer))*n_particle, '') + else: + wavelen = Quantity(start_wavelen, 'nm') + n_particle = Quantity(n_particle, '') + + n_matrix = Quantity(n_matrix, '') + m = index_ratio(n_particle, n_matrix) + x = size_parameter(wavelen, n_matrix, radius) + + if return_all: + return m, x, wavelen, radius, n_particle, n_matrix + else: + return m, x + + +def calc_coeffs(m, x): + """Utility function to calculate scattering coefficients""" + nstop = mie._nstop(x.max()) + coeffs = mie._scatcoeffs(m, x, nstop) + + return nstop, coeffs + + +def test_parameter_shapes(): + """Test to make sure vectorized size_parameter() and index_ratio() have the + right shapes""" + + num_wavelen = 8 + num_layer = 3 + wavelen = Quantity(np.linspace(400, 800, num_wavelen), 'nm') + radius = Quantity(np.linspace(0.85, 1.0, num_layer), 'um') + n_matrix = Quantity(1.00, '') + # let index be the same at all wavelengths, but different at each layer + n_particle = np.linspace(1.33, 1.59, num_layer) + n_particle = np.repeat(np.array([n_particle]), num_wavelen, axis=0) + n_particle = Quantity(n_particle, '') + + # multiple wavelengths, multiple layers. m and x should have shape + # [num_wavelen, num_layer]. + expected_shape = (num_wavelen, num_layer) + # The following should be true by construction of n_particle, but we test + # anyway to make sure that index_ratio() doesn't change shape + m = index_ratio(n_particle, n_matrix) + assert m.shape == expected_shape + # x should have shape [num_wavelen, num_layer] + x = size_parameter(wavelen, n_matrix, radius) + assert x.shape == expected_shape + + # one wavelength, multiple layers; index specified as 1D array. Should + # return a 1D index ratio and a 2D size parameter + wavelen = Quantity(400, 'nm') + num_layer = 6 + radius = Quantity(np.linspace(0.85, 1.0, num_layer), 'um') + n_particle = Quantity(np.linspace(1.33, 1.59, num_layer), '') + m = index_ratio(n_particle, n_matrix) + assert m.shape == (num_layer, ) + x = size_parameter(wavelen, n_matrix, radius) + assert x.shape == (1, num_layer) + + # one wavelength, multiple layers; index specified as 2D array with shape + # [1, num_layers]. Should return a 2D index ratio and a 2D size parameter + wavelen = Quantity(400, 'nm') + num_layer = 6 + radius = Quantity(np.linspace(0.85, 1.0, num_layer), 'um') + n_particle = Quantity(np.linspace(1.33, 1.59, num_layer)[np.newaxis,:], '') + m = index_ratio(n_particle, n_matrix) + assert m.shape == (1, num_layer) + x = size_parameter(wavelen, n_matrix, radius) + assert x.shape == (1, num_layer) + + # multiple wavelengths, one layer; index specified as a 2D array with shape + # [num_wavelen, 1]. Should return a 2D index ratio and 2D size parameter + num_wavelen = 8 + wavelen = Quantity(np.linspace(400, 800, num_wavelen), 'nm') + radius = Quantity(0.85, 'um') + n_particle = Quantity(np.linspace(1.33, 1.59, num_wavelen)[:,np.newaxis], + '') + m = index_ratio(n_particle, n_matrix) + assert m.shape == (num_wavelen, 1) + x = size_parameter(wavelen, n_matrix, radius) + assert x.shape == (num_wavelen, 1) + + +# all functions in this class are tested for various combinations of +# wavelengths and layers +@pytest.mark.parametrize("num_wavelen,num_layer", + [(1, 1), (10, 1), (1, 5), (10, 5)]) +class TestVectorizedSpecialFuncs(): + """Tests that simplifying/removing loops from Mie special functions + produces same results as using loops. Special functions and corresponding + tests of vectorization are as follows: + + riccati_psi_xi() : + not tested in this class, but tested implicitly in + `test_vectorized_scatcoeffs()` + lentz_dn1() : + tested by `test_lentz_dn1()` + dn_1_down() : + tested by `test_dn_1_down()` + log_der_13() : + not tested explicitly but tested implicitly in `test_Qratio()` and in + testing of the multilayer scattering coefficients and internal + scattering coefficients in `TestVectorizedUserFunctions` + Qratio() : + tested by `test_Qratio()` + R_psi() : + tested by `test_R_psi()` + + """ + mxargs = {"start_wavelen": 400, + "end_wavelen": 800, + "start_radius": 100, + "end_radius": 1000, + "start_n_particle": 1.59 + 0.001j, + "end_n_particle": 1.33 + 0.005j, + "n_matrix": Quantity(1.00, '')} + + # vectorizing functions may lead to small differences from loops, due to + # floating point precision. We set 10^-14 as a relative tolerance for + # differences, given that floating point errors tend to accumulate with + # multiple operations. + rtol = 1e-14 + + def test_lentz_dn1(self, num_wavelen, num_layer): + """Test whether Lentz continued fraction approximation for nth order + logarithmic derivative parallelizes properly with both wavelengths and + layers + + """ + m, x = mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs) + + # quick check on shapes + if np.isscalar(m): + assert (num_wavelen, num_layer) == (1, 1) + else: + assert m.shape == (num_wavelen, num_layer) + if np.isscalar(x): + assert (num_wavelen, num_layer) == (1, 1) + else: + assert x.shape == (num_wavelen, num_layer) + + # vectorized computation + nstop = mie._nstop(np.array(x).max()) + n = nstop + 1 + z = m * x + lentz_vec = mie_specfuncs.lentz_dn1(z, n) + + # looped computation + lentz = np.zeros((num_wavelen, num_layer), dtype=complex) + for i in range(num_wavelen): + for j in range(num_layer): + z = np.atleast_2d(m)[i, j] * np.atleast_2d(x)[i, j] + lentz[i, j] = mie_specfuncs.lentz_dn1(z, n).item() + + # vectorized and loop results should be exactly the same because the + # iteration count should be determined individually for each z, even in + # the vectorized version + assert_equal(lentz_vec, lentz) + + # check result against Lentz (1976) equation 9, which gives the ratio + # of Bessel functions of order nu = 9.5 at x = 1. Converting nu to n + # gives n = 9, and then noting that A_n = -n/z + ratio of Bessel + # functions, we add n/z to the result: + expected_ratio = 18.95228198 + assert_allclose(mie_specfuncs.lentz_dn1(1.0, 9) + 9, expected_ratio) + + def test_dn_1_down(self, num_wavelen, num_layer): + """Tests that down-recurrence for logarithmic derivatives can be + vectorized over wavelengths and layers. + + """ + m, x = mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs) + nstop = mie._nstop(np.array(x).max()) + nmx = nstop + 1 + + z = m * x + start_val = mie_specfuncs.lentz_dn1(z, nmx) + + # loop version of dn_1_down + dn = np.zeros(start_val.shape + (nmx+1,), dtype=complex) + dn[..., nmx] = start_val + for i in np.arange(nmx-1, -1, -1): + dn[..., i] = (i+1.)/z - 1.0/(dn[..., i+1] + (i+1.)/z) + dn = dn[..., 0:nstop+1] + + # vectorized version + dn_vec = mie_specfuncs.dn_1_down(z, nmx, nstop, start_val) + + # currently these differ at the 10^-15 level for some + # elements. The following test should pass: + assert_allclose(dn_vec, dn, rtol=self.rtol) + + def test_Qratio(self, num_wavelen, num_layer): + """Tests that vectorized version of Qratio (without loop for up + recursion) works the same as loop version. This test runs some of the + same code that is used to calculate the scattering coefficients for + multilayer particles. + + """ + def Qratio_loop(z1, z2, nstop, dns1, dns2): + # non-vectorized (loop-based) version of Qratio calculation, from + # previous version of pymie + d1z1 = dns1[0] + d3z1 = dns1[1] + d1z2 = dns2[0] + d3z2 = dns2[1] + + # initialize according to Yang eqn. 34 + a1 = np.real(z1) + a2 = np.real(z2) + b1 = np.imag(z1) + b2 = np.imag(z2) + qns = np.zeros(z1.shape + (nstop+1,), dtype=complex) + qns[..., 0] = (np.exp(-2.*(b2-b1)) * + (np.exp(-1j*2.*a1)-np.exp(-2.*b1)) + / (np.exp(-1j*2.*a2) - np.exp(-2.*b2))) + for i in np.arange(1, nstop+1): + qns[..., i] = qns[..., i-1]* (((d3z1[..., i] + i/z1) + * (d1z2[..., i] + i/z2)) + / ((d3z2[..., i] + i/z2) + * (d1z1[..., i] + i/z1))) + return qns + + m, x = mx(num_wavelen=num_wavelen, num_layer=num_layer, + start_wavelen=self.mxargs["start_wavelen"], + end_wavelen=self.mxargs["end_wavelen"], start_radius=850, + end_radius=1000, + start_n_particle=1.33, + end_n_particle=1.59) + + marray = np.atleast_1d(m).astype(complex) + xarray = np.atleast_1d(x).astype(complex) + + nstop = mie._nstop(xarray.max()) + # m_l x_{l-1} + z1 = marray[..., 1:] * xarray[..., :-1] + # # m_l x_l + z2 = marray[..., 1:] * xarray[..., 1:] + + # pre-calculate logarithmic derivatives for all layers + derz1s = mie_specfuncs.log_der_13(z1, nstop) + derz2s = mie_specfuncs.log_der_13(z2, nstop) + + # vectorized calculation of Q_n^l for all layers + Qnl_vec = mie_specfuncs.Qratio(z1, z2, nstop, dns1 = derz1s, + dns2 = derz2s) + # non-vectorized (loop-based) calculation of Q_n^l + Qnl_loop = Qratio_loop(z1, z2, nstop, derz1s, derz2s) + + assert_allclose(Qnl_vec.real, Qnl_loop.real, rtol=self.rtol) + # Use a different test to look at imaginary elements because the + # differences are pretty close to zero but can vary a lot in their + # magnitudes. + assert_allclose(np.abs(Qnl_vec.imag - Qnl_loop.imag), 0, atol=1e-10) + + def test_R_psi(self, num_wavelen, num_layer): + """Tests that the up-recurrence in the calculation of the ratio of + Riccati-Bessel functions can be done without a loop + + """ + # The R_psi calculation shows up in the calculation of internal + # coefficients, which is only valid (for now) for non-multilayer + # spheres. Nonetheless, the R_psi calculation can be vectorized over + # layers, and we test that it vectorizes properly over both wavelength + # and layers here + m, x = mx(num_wavelen, num_layer) + z1 = x + z2 = m * x + nstop = mie._nstop(np.array(x).max()) + nmax = nstop + 1 + + # note that inputs to R_psi must be 2-dimensional + if np.isscalar(z1): + z1 = z1 * np.ones((1,1)) + z2 = z2 * np.ones((1,1)) + + # vectorized version + output_vec = mie_specfuncs.R_psi(z1, z2, nmax) + + # loop version of R_psi (from a previous version of pymie) + output = np.zeros(z1.shape + (nmax + 1,), dtype=complex) + output[..., 0] = np.sin(z1) / np.sin(z2) + dnz1 = mie_specfuncs.dn_1_down(z1, nmax + 1, nmax, + mie_specfuncs.lentz_dn1(z1, nmax + 1)) + dnz2 = mie_specfuncs.dn_1_down(z2, nmax + 1, nmax, + mie_specfuncs.lentz_dn1(z2, nmax + 1)) + for i in np.arange(1, nmax + 1): + output[..., i] = output[..., i-1] * ((dnz2[..., i] + i / z2) + / (dnz1[..., i] + i / z1)) + + assert_allclose(output_vec.imag, output.imag, rtol=self.rtol) + assert_allclose(output_vec.real, output.real, rtol=self.rtol) + +class TestVectorizedInternalFunctions(): + """Test vectorization of the internal Mie calculation functions (the ones + starting with an underscore) over wavelength and layers. These tests check + primarily that the functions return the same values for array arguments as + they do when we loop over the arrays. They do not check for correctness + of the results. + + Internal functions and corresponding tests of vectorization are as follows: + + _pis_and_taus() : + not tested explicitly here, but tested implicitly in + `test_vectorized_calc_ang_dist()`. Vectorization over angles is tested + in `test_mie.py::test_pis_taus()` + _scatcoeffs() : + tested by `test_vectorized_scatcoeffs()` + _scatcoeffs_multi() : + tested by `test_vectorized_scatcoeffs()` + _internal_coeffs() : + tested by `test_vectorized_internal_coeffs()` + _trans_coeffs() : + * vectorization not yet tested + _time_coeffs() : + * vectorization not yet tested + _W0() : + * vectorization not yet tested + _nstop() : + tested by `test_vectorized_nstop()` + _asymmetry_parameter() : + not tested explicitly here, but tested implicitly in + `test_vectorized_asymmetry_parameter()`, which tests the user-facing + function for calculating asymmetry parameters + _cross_sections() : + not tested explicitly here, but tested implicitly in + `test_vectorized_cross_sections()`, which tests the user-facing + function for calculating cross-sections + _cross_sections_complex_medium_fu() : + tested by `test_vectorized_cross_sections_complex_medium()` + _cross_sections_complex_medium_sudiarta() : + tested by `test_vectorized_cross_sections_complex_medium()` + _scat_fields_complex_medium() : + * vectorization not yet tested + diff_scat_intensity_complex_medium() : + tested by `test_vectorized_angular_functions()` + integrate_intensity_complex_medium() : + tested by `test_vectorized_angular_functions()` + diff_abs_intensity_complex_medium() : + * vectorization not yet tested + amplitude_scattering_matrix() : + tested by `test_vectorized_angular_functions()` + vector_scattering_amplitude() : + tested by `test_vectorized_angular_functions()` + _amplitude_scattering_matrix() : + not tested explicitly here, but tested implicitly in + `test_vectorized_calc_ang_dist()` + _amplitude_scattering_matrix_RG() : + not tested explicitly here, but tested implicitly in + `test_vectorized_calc_ang_dist()` + + """ + mxargs = {"start_wavelen": 400, + "end_wavelen": 800, + "start_radius": 100, + "end_radius": 1000, + "start_n_particle": 1.59 + 0.001j, + "end_n_particle": 1.33 + 0.005j, + "n_matrix": Quantity(1.00, '')} + + num_theta = 20 + thetas = Quantity(np.linspace(0, np.pi, num_theta), 'rad') + + @pytest.mark.parametrize("num_wavelen", [1, 10, 100]) + def test_vectorized_nstop(self, num_wavelen): + # Just checks that the shape of nstop is correct + # (should scale with number of wavelengths) + m, x = mx(num_wavelen=num_wavelen, num_layer=1, **self.mxargs) + nstop = mie._nstop(x) + assert np.atleast_1d(nstop).shape == (num_wavelen,) + + @pytest.mark.parametrize("num_wavelen,num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_scatcoeffs(self, num_wavelen, num_layer): + """Tests that mie._scatcoeffs() and mie._scatcoeffs_multi() vectorize + properly. + + """ + m, x = mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs) + nstop, coeffs = calc_coeffs(m, x) + + # if multilayer, check that _scatcoeffs is actually calling the + # multilayer code + if num_layer > 1: + coeffs_direct = mie._scatcoeffs_multi(m, x) + assert_equal(coeffs, coeffs_direct) + + # make sure shape is correct + if num_wavelen == 1: + expected_shape = (2, nstop) + assert coeffs.shape == expected_shape + # no further test since no loop required in this case + else: + expected_shape = (2, num_wavelen, nstop) + assert coeffs.shape == expected_shape + + # we should get same value from loop + coeffs_loop = np.zeros(expected_shape, dtype=complex) + for i in range(m.shape[0]): + if num_layer == 1: + c = mie._scatcoeffs(m[i], x[i], nstop) + else: + # need to specify nstop here; otherwise we will get a + # different number of scattering coefficients for each + # wavelength, since _scatcoeffs_multi() picks the largest x + # for each wavelength. + c = mie._scatcoeffs_multi(m[i], x[i], nstop) + coeffs_loop[:, i] = c + assert_equal(coeffs, coeffs_loop) + + @pytest.mark.parametrize("num_wavelen,num_layer", + [(1, 1), (10, 1), (1, 5), (10, 5)]) + def test_vectorized_internal_coeffs(self, num_wavelen, num_layer): + """Tests that mie._internal_coeffs() vectorizes properly + + """ + m, x = mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs) + m = np.atleast_1d(m) + x = np.atleast_1d(x) + nstop = mie._nstop(x.max()) + + # should not work for a layered sphere + if np.atleast_2d(m).shape[-1] > 1: + with pytest.raises(ValueError, match="Internal Mie coefficients"): + coeffs = mie._internal_coeffs(m, x, nstop) + + else: + coeffs = mie._internal_coeffs(m, x, nstop) + + # make sure shape is correct + if num_wavelen == 1: + expected_shape = (2, nstop) + assert coeffs.shape == expected_shape + # no further test since no loop required in this case + else: + expected_shape = (2, num_wavelen, nstop) + assert coeffs.shape == expected_shape + + # we should get same values from loop + coeffs_loop = np.zeros(expected_shape, dtype=complex) + for i in range(m.shape[0]): + c = mie._internal_coeffs(m[i], x[i], nstop) + coeffs_loop[:, i] = c + + assert_equal(coeffs, coeffs_loop) + + @pytest.mark.parametrize("n_matrix", + [1.33, pytest.param(1.33+0.001j), + pytest.param(1.33+0.1j)]) + @pytest.mark.parametrize("num_wavelen,num_layer", + [(1, 1), (10, 1), (1, 5), (10, 5)]) + def test_vectorized_cross_sections_complex_medium(self, num_wavelen, + num_layer, n_matrix): + """ + Test the vectorized versions of the (Fu, Sun) and (Sudiarta, Chylek) + solutions for the cross-sections in absorbing medium: + mie._cross_sections_complex_medium_fu() + mie._cross_sections_complex_medium_sudiarta() + Note that these functions should work properly only for non-layered + particles, but we test that the Sudiarta function runs anyway for + layered spheres (it will take the largest radius). + """ + mxargs = self.mxargs.copy() + mxargs['n_matrix'] = n_matrix + m, x, wavelen, radius, n_particle, n_matrix = \ + mx(num_wavelen=num_wavelen, num_layer=num_layer, **mxargs, + return_all=True) + + m = np.atleast_1d(m) + x = np.atleast_1d(x) + nstop = mie._nstop(x.max()) + + # Calculate coefficients. For the Fu solution, the internal + # coefficiences cl and dl are needed for the calculation of the + # absorption and extinction cross sections. + al, bl = mie._scatcoeffs(m, x, nstop) + if num_layer == 1: + x_scat = size_parameter(wavelen, n_particle, radius) + cl, dl = mie._internal_coeffs(m, x, nstop) + c_fu = mie._cross_sections_complex_medium_fu(al, bl, cl, dl, + radius, n_particle, + n_matrix, x_scat, x, + wavelen) + + c_sudiarta = mie._cross_sections_complex_medium_sudiarta(al, bl, + x, radius) + + + c_sud_sca = np.zeros(num_wavelen) + c_sud_abs = np.zeros_like(c_sud_sca) + c_sud_ext = np.zeros_like(c_sud_sca) + c_fu_sca = np.zeros(num_wavelen) + c_fu_abs = np.zeros_like(c_fu_sca) + c_fu_ext = np.zeros_like(c_fu_sca) + wavelen = np.atleast_1d(wavelen) + al = np.reshape(al, (num_wavelen, al.shape[-1])) + bl = np.reshape(bl, (num_wavelen, bl.shape[-1])) + if num_layer == 1: + cl = np.reshape(cl, (num_wavelen, cl.shape[-1])) + dl = np.reshape(dl, (num_wavelen, dl.shape[-1])) + for i in range(num_wavelen): + c_sud_loop = \ + mie._cross_sections_complex_medium_sudiarta(al[i], bl[i], x[i], + radius) + c_sud_sca[i] = c_sud_loop[0].magnitude + c_sud_abs[i] = c_sud_loop[1].magnitude + c_sud_ext[i] = c_sud_loop[2].magnitude + + x_scat = size_parameter(wavelen[i], n_particle[i], radius) + if num_layer == 1: + c_fu_loop = \ + mie._cross_sections_complex_medium_fu(al[i], bl[i], cl[i], + dl[i], radius, + n_particle[i], + n_matrix, + x_scat, x[i], + wavelen[i]) + c_fu_sca[i] = c_fu_loop[0].magnitude + c_fu_abs[i] = c_fu_loop[1].magnitude + c_fu_ext[i] = c_fu_loop[2].magnitude + + + # the loop calculations might differ at the 1e-15 level for Sudiarta + rtol = 1e-15 + assert_allclose(c_sudiarta[0].magnitude, c_sud_sca, rtol=rtol) + assert_allclose(c_sudiarta[1].magnitude, c_sud_abs, rtol=rtol) + assert_allclose(c_sudiarta[2].magnitude, c_sud_ext, rtol=rtol) + if num_layer == 1: + assert_equal(c_fu[0].magnitude, c_fu_sca) + assert_equal(c_fu[1].magnitude, c_fu_abs) + assert_equal(c_fu[2].magnitude, c_fu_ext) + + # also check that the Fu and Sudiarta cross sections agree for both + # real and complex matrix indices + for i in range(len(c_fu)): + assert_allclose(c_fu[i].magnitude, c_sudiarta[i].magnitude, + rtol=1e-10) + + @pytest.mark.parametrize("num_wavelen,num_layer", + [(1, 1), (10, 1), (1, 5), (10, 5)]) + @pytest.mark.parametrize("coordinate_system", ["scattering plane", + "cartesian"]) + def test_vectorized_angular_functions(self, num_wavelen, + num_layer, + coordinate_system): + """Tests that mie.vector_scattering_amplitude(), + mie.amplitude_scattering_matrix(), + diff_scat_intensity_complex_medium(), and + integrate_intensity_complex_medium() vectorize properly + + TODO: test vectorized near-field calculation in + diff_scat_intensity_complex_medium() + + """ + m, x, wavelen, radius, n_particle, n_matrix = \ + mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs, + return_all=True) + + if coordinate_system == "scattering plane": + phis = None + thetas = self.thetas + else: + phis = Quantity(np.linspace(0, 2*np.pi, self.num_theta), '') + thetas, phis = np.meshgrid(self.thetas, phis) + + vsa = mie.vector_scattering_amplitude(m, x, thetas, + coordinate_system = + coordinate_system, phis = phis) + + mat = mie.amplitude_scattering_matrix(m, x, thetas, + coordinate_system = + coordinate_system, + phis = phis) + + # choose distance reasonably close to the particle for differential + # scattering calculations + k = 2*np.pi*n_matrix/wavelen + d = 10*np.atleast_1d(radius)[-1] + kd = np.atleast_1d(k*d) + i12 = mie.diff_scat_intensity_complex_medium(m, x, thetas, + kd, + coordinate_system = + coordinate_system, + phis = phis) + + integral = mie.integrate_intensity_complex_medium(*i12, d, + thetas, k, phi_min = Quantity(0.0, 'rad'), + phi_max = Quantity(2*np.pi, 'rad'), + coordinate_system = coordinate_system, phis = phis) + + # check that shapes of all the computed quantities are correct + for element in vsa + mat + i12: + if num_wavelen > 1: + assert element.shape == (num_wavelen, ) + thetas.shape + else: + assert element.shape == thetas.shape + + # check that vectorized calculations match looped calculations over + # scalars + amp0 = np.zeros((num_wavelen, ) + thetas.shape, dtype=complex) + amp1 = np.zeros_like(amp0) + S1 = np.zeros((num_wavelen, ) + thetas.shape, dtype=complex) + S2 = np.zeros_like(S1) + S3 = np.zeros_like(S1) + S4 = np.zeros_like(S1) + i1 = np.zeros_like(S1) + i2 = np.zeros_like(S1) + + sigma = np.zeros(num_wavelen) + sigma_1 = np.zeros_like(sigma) + sigma_2 = np.zeros_like(sigma) + dsigma_1 = np.zeros((num_wavelen, ) + thetas.shape) + dsigma_2 = np.zeros_like(dsigma_1) + + m = np.atleast_1d(m) + x = np.atleast_1d(x) + k = np.atleast_1d(k) + for i in range(num_wavelen): + mat_loop = mie.amplitude_scattering_matrix(m[i], x[i], thetas, + coordinate_system = + coordinate_system, + phis = phis) + + vsa_loop = mie.vector_scattering_amplitude(m[i], x[i], thetas, + coordinate_system = + coordinate_system, + phis = phis) + i_loop = mie.diff_scat_intensity_complex_medium(m[i], x[i], + thetas, + kd[i], + coordinate_system = + coordinate_system, + phis = phis) + + integral_loop = mie.integrate_intensity_complex_medium(*i_loop, d, + thetas, k[i], phi_min = Quantity(0.0, 'rad'), + phi_max = Quantity(2*np.pi, 'rad'), + coordinate_system = coordinate_system, phis = phis) + + S1[i], S2[i], S3[i], S4[i] = mat_loop + amp0[i], amp1[i] = vsa_loop + i1[i], i2[i] = i_loop + # can't assign using tuple notation because pint will confuse numpy + # by trying to assign units to each element, so we have to take + # the magnitudes here + units = integral_loop[0].units + sigma[i] = integral_loop[0].magnitude + sigma_1[i] = integral_loop[1].magnitude + sigma_2[i] = integral_loop[2].magnitude + dsigma_1[i] = integral_loop[3].magnitude + dsigma_2[i] = integral_loop[4].magnitude + + assert_equal(mat[0], S1.squeeze()) + assert_equal(mat[1], S2.squeeze()) + assert_equal(mat[2], S3.squeeze()) + assert_equal(mat[3], S4.squeeze()) + + assert_equal(vsa[0], amp0.squeeze()) + assert_equal(vsa[1], amp1.squeeze()) + + assert_equal(i12[0], i1.squeeze()) + assert_equal(i12[1], i2.squeeze()) + + assert_equal(integral[0].magnitude, sigma.squeeze()) + assert_equal(integral[1].magnitude, sigma_1.squeeze()) + assert_equal(integral[2].magnitude, sigma_2.squeeze()) + assert_equal(integral[3].magnitude, dsigma_1.squeeze()) + assert_equal(integral[4].magnitude, dsigma_2.squeeze()) + + +class TestVectorizedUserFunctions(): + """Test vectorization of the user-facing Mie calculation functions over + wavelength for solid (one layer) spheres. These tests check primarily that + the functions return the same values for array arguments as they do when + we loop over the arrays. They do not check for correctness of the + results. + + User functions and corresponding tests of vectorization are as follows: + + calc_ang_dist() : + tested by test_vectorized_calc_ang_dist() + calc_cross_sections() : + tested by test_vectorized_cross_sections() + calc_efficiencies() : + tested by test_vectorized_calc_efficiencies() + calc_g() : + tested by test_vectorized_asymmetry_parameter() + calc_integrated_cross_section() : + tested by test_vectorized_cross_sections() + calc_energy() : + * vectorization not yet tested + calc_dwell_time() : + * vectorization not yet tested + calc_reflectance() : + tested by test_vectorized_calc_reflectance() + + """ + mxargs = {"start_wavelen": 400, + "end_wavelen": 800, + "start_radius": 100, + "end_radius": 1000, + "start_n_particle": 1.59 + 0.001j, + "end_n_particle": 1.33 + 0.005j, + "n_matrix": Quantity(1.00, '')} + + num_angle = 19 + angles = Quantity(np.linspace(0, 180., num_angle), 'deg') + + @pytest.mark.parametrize("num_wavelen, num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_asymmetry_parameter(self, num_wavelen, num_layer): + """Tests that mie.calc_g() vectorizes properly. Also implicitly checks + that mie._asymmetry_parameter() vectorizes properly + + """ + m, x = mx(num_wavelen, num_layer, **self.mxargs) + # make sure shape is [num_wavelen] + g = mie.calc_g(m,x) + if num_wavelen == 1: + expected_shape = () + assert g.shape == expected_shape + # no further test needed since no loop is required in this case + else: + expected_shape = (num_wavelen,) + assert g.shape == expected_shape + + # we should get same values from loop. Need to set nstop to the + # same value as used in the vectorized calculation. + g_loop = np.zeros(expected_shape, dtype=float) + nstop = mie._nstop(x.max()) + for i in range(num_wavelen): + g_loop[i] = mie.calc_g(m[i], x[i], nstop=nstop) + assert_equal(g, g_loop) + + @pytest.mark.parametrize("num_wavelen, num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_cross_sections(self, num_wavelen, num_layer): + """Tests that mie.calc_cross_sections() and + mie.calc_integrated_cross_sections() vectorize properly. Also + implicitly checks that _cross_sections() vectorizes properly + + """ + m, x = mx(num_wavelen, num_layer, **self.mxargs) + # wavelength in medium + wavelen = Quantity(np.linspace(self.mxargs["start_wavelen"], + self.mxargs["end_wavelen"], + num_wavelen), + "nm") + wavelen_med = wavelen/self.mxargs["n_matrix"] + cscat, cext, cback, cabs, asym = mie.calc_cross_sections(m, x, + wavelen_med) + + # test shape + expected_shape = (num_wavelen,) + for cs in [cscat, cext, cback, cabs, asym]: + assert cs.shape == expected_shape + + # we should get same values from loop + cscat_loop = np.zeros(expected_shape, dtype=float) + cext_loop = np.zeros(expected_shape, dtype=float) + cback_loop = np.zeros(expected_shape, dtype=float) + cabs_loop = np.zeros(expected_shape, dtype=float) + asym_loop = np.zeros(expected_shape, dtype=float) + for i in range(num_wavelen): + cs = mie.calc_cross_sections(m[i], x[i], wavelen[i]) + cscat_loop[i], cext_loop[i], cback_loop[i], \ + cabs_loop[i], asym_loop[i] = (c.magnitude for c in cs) + assert_equal(cscat.magnitude, cscat_loop) + assert_equal(cext.magnitude, cext_loop) + assert_equal(cback.magnitude, cback_loop) + assert_equal(cabs.magnitude, cabs_loop) + assert_equal(asym.magnitude, asym_loop) + + # check that numerical integration works too, and that it gives a + # similar value for the total cross section + num_angles = 100 + thetas = Quantity(np.linspace(0, np.pi, num_angles), 'rad') + c_integrated = mie.calc_integrated_cross_section(m, x, wavelen_med, + thetas) + + c_integrated = c_integrated.to(cscat.units) + # small grid gives large integration error, but should still be within + # 10% + assert_allclose(c_integrated.magnitude, cscat.magnitude, rtol=1e-1) + assert c_integrated.shape == (num_wavelen,) + + @pytest.mark.parametrize("num_wavelen, num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_calc_efficiencies(self, num_wavelen, num_layer): + """Tests that mie.calc_efficiencies() vectorizes properly with + wavelength, including with multi-layered particles. + + """ + m, x = mx(num_wavelen, num_layer, **self.mxargs) + qscat, qext, qback = mie.calc_efficiencies(m, x) + + # test shape + if num_wavelen == 1: + expected_shape = () + # no further test because no loop is needed in this case + else: + expected_shape = (num_wavelen,) + for q in [qscat, qext, qback]: + assert q.shape == expected_shape + + if num_wavelen > 1: + # we should get same values from loop + qscat_loop = np.zeros(expected_shape, dtype=float) + qext_loop = np.zeros(expected_shape, dtype=float) + qback_loop = np.zeros(expected_shape, dtype=float) + for i in range(num_wavelen): + qs = mie.calc_efficiencies(m[i], x[i]) + qscat_loop[i], qext_loop[i], qback_loop[i] = (q for q in qs) + assert_equal(qscat, qscat_loop) + assert_equal(qext, qext_loop) + assert_equal(qback, qback_loop) + + @pytest.mark.parametrize("num_wavelen, num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_calc_ang_dist(self, num_wavelen, num_layer): + """Tests that mie.calc_ang_dist() vectorizes properly. Also implicitly + checks that _amplitude_scattering_matrix() and + _amplitude_scattering_matrix_RG() vectorize properly. Also checks for + correctness of RG calculations by comparing against Mie calculations + for small refractive index difference. + + """ + m, x = mx(num_wavelen, num_layer, **self.mxargs) + form_factor = mie.calc_ang_dist(m, x, self.angles) + if num_wavelen == 1: + expected_shape = (self.num_angle,) + for pol in form_factor: + assert pol.shape == expected_shape + # no further test required since there is only one wavelength + return + else: + expected_shape = (num_wavelen, self.num_angle) + for pol in form_factor: + assert pol.shape == expected_shape + + # we should get same values from loop + ipar_loop = np.zeros(expected_shape, dtype=float) + iperp_loop = np.zeros(expected_shape, dtype=float) + for i in range(num_wavelen): + ipar, iperp = mie.calc_ang_dist(m[i], x[i], self.angles) + ipar_loop[i] = ipar + iperp_loop[i] = iperp + assert_equal(form_factor[0], ipar_loop) + assert_equal(form_factor[1], iperp_loop) + + # check vectorization for Rayleigh-Gans approximation + if num_layer > 1: + with pytest.raises(ValueError, + match="Rayleigh-Gans approximation cannot"): + form_factor_RG = mie.calc_ang_dist(m, x, self.angles, + mie=False) + return + + form_factor_RG = mie.calc_ang_dist(m, x, self.angles, + mie=False) + + expected_shape = (num_wavelen, self.num_angle) + for pol in form_factor_RG: + assert pol.shape == expected_shape + + # also check that we recover approximately the same result for RG as we + # do for Mie in the limit of low refractive index + radius = Quantity('0.85 um') + n_matrix = Quantity(1.00, '') + # let index be the same at all wavelengths. We look at small index + # contrast (1 + 1e-8) to be in the RG regime. If we go smaller we run + # into numerical issues + n_particle = Quantity(np.ones(num_wavelen)*(1 + 1e-8), '') + wavelen = Quantity(np.linspace(self.mxargs["start_wavelen"], + self.mxargs["end_wavelen"], + num_wavelen), + 'nm') + m = index_ratio(n_particle, n_matrix)[:, np.newaxis] + x = size_parameter(wavelen, n_matrix, radius) + num_angle = 1000 + # 0 degree scattering may give differences between RG and Mie, so we + # compare at a few degrees and higher; also we do a lot of angles to + # capture the sharp dips in the form factor + angles = Quantity(np.linspace(10, 180., num_angle), 'deg') + form_factor_RG = mie.calc_ang_dist(m, x, angles, mie=False) + form_factor_mie = mie.calc_ang_dist(m, x, angles) + + # Since we are comparing small numbers at the dips of the form factor, + # the Mie and RG solutions may have a relative difference of up to a + # few percent. But the absolute difference should be very small + # (smaller than the default atol for this test). + assert_allclose(form_factor_RG, form_factor_mie, rtol=1e-1) + + @pytest.mark.parametrize("n_medium", + [1.33, pytest.param(1.33+0.1j)]) + @pytest.mark.parametrize("num_wavelen, num_layer", + [(10, 1), (1, 5), (10, 5)]) + def test_vectorized_calc_reflectance(self, n_medium, num_wavelen, + num_layer): + """Tests that vectorized calc_reflectance() returns same result as + loop. + + """ + m, x, wavelen, radius, n_particle, n_matrix = \ + mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs, + return_all=True) + n_medium = Quantity(n_medium, '') + refl = mie.calc_reflectance(radius, n_medium, n_particle, wavelen) + + refl_loop = np.zeros(num_wavelen, dtype=complex) + wavelen = np.atleast_1d(wavelen) + n_particle = np.atleast_1d(n_particle) + for i in range(num_wavelen): + reflectance = mie.calc_reflectance(radius, n_medium, n_particle[i], + wavelen[i]).magnitude + refl_loop[i] = reflectance + + assert_allclose(refl.magnitude, refl_loop, rtol=1e-14) + assert refl.units == 1/wavelen.units**2 diff --git a/pymie/tests/test_multilayer.py b/pymie/tests/test_multilayer.py index 4b32f1c..d2c58e5 100755 --- a/pymie/tests/test_multilayer.py +++ b/pymie/tests/test_multilayer.py @@ -22,7 +22,6 @@ """ import os from .. import Quantity, size_parameter, np, mie -from .. import multilayer_sphere_lib as msl from numpy.testing import assert_array_almost_equal, assert_allclose import yaml @@ -44,7 +43,7 @@ def test_scatcoeffs_multi(): marray = [1.15, 1.15] # layer index ratios, innermost first multi_radius = Quantity(np.array([100.0, 100.0]),'nm') xarray = size_parameter(wavelen, n_sample, multi_radius) - coeffs_multi = msl.scatcoeffs_multi(marray, xarray) + coeffs_multi = mie._scatcoeffs_multi(marray, xarray) assert_array_almost_equal(coeffs, coeffs_multi) @@ -53,7 +52,7 @@ def test_scatcoeffs_multi(): marray2 = [1.15, 1.15, 1.15] # layer index ratios, innermost first multi_radius2 = Quantity(np.array([100.0, 100.0, 100.0]),'nm') xarray2 = size_parameter(wavelen, n_sample, multi_radius2) - coeffs_multi2 = msl.scatcoeffs_multi(marray2, xarray2) + coeffs_multi2 = mie._scatcoeffs_multi(marray2, xarray2) assert_array_almost_equal(coeffs, coeffs_multi2) @@ -67,8 +66,8 @@ def test_scatcoeffs_multi_absorbing_particle(): multi_radius = Quantity(np.array([100.0, 110.0]),'nm') xarray = size_parameter(wavelen, n_sample, multi_radius) - coeffs_multi_real = msl.scatcoeffs_multi(marray_real, xarray) - coeffs_multi_imag = msl.scatcoeffs_multi(marray_imag, xarray) + coeffs_multi_real = mie._scatcoeffs_multi(marray_real, xarray) + coeffs_multi_imag = mie._scatcoeffs_multi(marray_imag, xarray) assert_array_almost_equal(coeffs_multi_real, coeffs_multi_imag) @@ -80,13 +79,18 @@ def test_sooty_particles(): We will use the data in [Yang2003]_ Table 3 on p. 1717, cases 2, 3, and 4 as our gold standard. ''' + # size parameter corresponding to outer radius x_L = 100 + + # index ratios m_med = 1.33 m_abs = 2. + 1.j + + # volume fraction f_v = 0.1 def efficiencies_from_scat_units(m, x): - asbs = msl.scatcoeffs_multi(m, x) + asbs = mie._scatcoeffs_multi(m, x) qs = np.array(mie._cross_sections(*asbs)) * 2 / x_L**2 # there is a factor of 2 conventional difference between # "backscattering" and "radar backscattering" efficiencies. @@ -118,5 +122,14 @@ def efficiencies_from_scat_units(m, x): rtol = 1e-3) assert_allclose(efficiencies_from_scat_units(m_as, x_as), gold[1], rtol = 2e-5) - assert_allclose(efficiencies_from_scat_units(m_sm, x_sm), gold[2], - rtol = 1e-3) + sooty_efficiencies = efficiencies_from_scat_units(m_sm, x_sm) + assert_allclose(sooty_efficiencies, gold[2], rtol = 1e-3) + + # also ensure that we get the same efficiencies from calc_efficiencies() + # for the 900-layer particle + qscat, qext, qback = mie.calc_efficiencies(m_sm, x_sm) + # Rearrange to match order in gold file. We multiply qback by 4 pi to get + # the radar backscattering efficiency (factor of 2 is included in + # calc_efficiencies()) + efficiencies = np.array([qext, qscat, qback*4*np.pi]) + assert_allclose(efficiencies, sooty_efficiencies) diff --git a/pyproject.toml b/pyproject.toml index 8228754..5942298 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,28 @@ +[build-system] +requires = ["setuptools"] +build-backend = "setuptools.build_meta" + +[project] +name = "pymie" +version = "0.2" +description = "Pure Python package for Mie scattering calculations" +readme = "README.md" +authors = [ + { name = "Manoharan Lab, Harvard University", email = "vnm@seas.harvard.edu" }, +] +dependencies = [ + "numpy", + "pint", + "scipy", + # pytest and pyyaml are needed for running unit tests + "pytest", + "pyyaml", +] +license-files = ["LICENSE",] + +[project.urls] +Homepage = "https://github.com/manoharan-lab/python-mie" + # convert all warnings to errors [tool.pytest.ini_options] filterwarnings = [ @@ -23,4 +48,7 @@ select = [ "**/tests/*" = ["E501", "F841"] # Ignore "ambiguous variable name" when we use "l" as a variable in Mie # calculations -"pymie/mie.py" = ["E741"] \ No newline at end of file +"pymie/mie.py" = ["E741"] + +[tool.setuptools.packages.find] +include = ["pymie"] diff --git a/setup.py b/setup.py deleted file mode 100644 index 2827b89..0000000 --- a/setup.py +++ /dev/null @@ -1,10 +0,0 @@ -from setuptools import setup - -setup(name='pymie', - version='0.1', - description='Pure Python package for Mie scattering calculations.', - url='https://github.com/manoharan-lab/python-mie', - author='Manoharan Lab, Harvard University', - author_email='vnm@seas.harvard.edu', - packages=['pymie'], - install_requires=['pint', 'numpy', 'scipy'])