Add method to compute scattering phase matrix

[nollety]

Changes

I've added a method named mie_phase_matrix that computes the scattering phase matrix. The method has the same signature than mie_S1_S2, namely:

def mie_phase_matrix(m, x, mu, nom="albedo"):
    ...

The phase scattering matrix is derived from the scattering amplitude functions according to equations 5.2.105-6 in K. N. Liou (2002) - An Introduction to Atmospheric Radiation, Second Edition.

I am not familiar with numba, so I just duplicated mie_phase_matrix in the modules miepython and miepython_nojit, but I assume mie_phase_matrix could be optimised in miepython?

Usage

For example, the following:

import miepython
import numpy as np

miepython.mie_phase_matrix(
    m=1.5,
    x=5.0,
    mu=np.linspace(-1, 1, 1000),
)

returns an array with shape (4, 4, 1000), the first axis corresponding to the matrix row index, the second axis corresponding to the matrix column index and the last axis corresponding to the scattering angle' cosine (mu). If mu is a scalar, it returns a (4, 4) array, i.e., the last axis is squeezed. The (0, 0) element of the scattering phase matrix provides the scattering phase function.

Checklist

• [x] Tests passed (python=3.9)

To-do

• [ ] Update changelog
• [ ] Add documentation
• [ ] Add tests

[nollety]

@scottprahl I must do a bit of research to find test data to validate this implementation. I will come back to this pull request as soon as I find the test data. In the meantime, I have marked this pull request as draft.

[nollety]

I tried to measure the performance of the jit version of miepython versus the non-jit version with respect to mie_phase_matrix, similar to what is done here and I found the jit improvement to be 7 to 371X (versus 31 to 501X for the scattering phase function) on my machine. These numbers don't look fantastic, I assume because my implementation is not optimised?

[scottprahl]

Your interface looks good.

Please make sure it returns a simple 4x4 matrix when a single scalar mu value is passed, e.g.,

 miepython.mie_phase_matrix(1.5, 2.0, 0.5)

Finding a few test cases would be fantastic. Obviously, you can test basic scattering matrix symmetry.

[scottprahl]

As far as jit/no-jit speeds go, I did the simplest possible numba acceleration. It is highly likely that other speed optimizations are possible.

Since it makes no sense to jit the mie_phase_matrix() function, your implementation looks sane. The speeds are what they are.

[nollety]

Please make sure it returns a simple 4x4 matrix when a single scalar mu value is passed, e.g.,

 miepython.mie_phase_matrix(1.5, 2.0, 0.5)

Ah yes, I had forgotten about that point, I'll refactor accordingly.

[scottprahl]

I guess I should have responded here and closed the other issue. Any progress?