[riclarsson]

AI message by request:

Scattering on rtepack types; T-matrix in C++ with azimuthally random orientation; polarization and accuracy fixes

Besides moving the scattering containers to rtepack types (bit-identical, no physics change), this PR changes
scattering physics in four places and fixes several bugs that gave wrong numbers without any error.

Physics

ARTS's Stokes convention for particle scattering (behaviour change). Mie and the T-matrix returned scattering
matrices in Mishchenko's (and Bohren & Huffman's) convention, whose U has the opposite sign to ARTS's. As a result,
F34 had the wrong sign, and so did Z13, Z23, Z31, Z32, Z34 and Z43 of fixed orientation. V produced by scattering of
linearly polarized radiation therefore had the wrong sign relative to the rest of ARTS (Zeeman, surfaces). Both are
now in ARTS's convention: the T-matrix phase matrix is the Mueller matrix of its amplitude matrix, and F34 changes
sign in Mie, the T-matrix and their habits.

• ARTS's Mie F34 now equals Evans and Stephens (1991), Table 3. Before, it was exactly -1 times that.
• When scattering is the only source of V, I and Q are unchanged, since they depend on F34 only through even powers.
  When Zeeman or surfaces also produce V, I and Q change at the order of that coupling.

T-matrix in C++. Mishchenko's solver is ported, and the Fortran is removed. The port computes in double precision
throughout (the original stored the T-matrix in single precision). It evaluates the amplitude at the angles given,
and it computes the orientation average with exact Clebsch-Gordan coefficients (wigxjpf).

Azimuthally random orientation (new). Particles with a distribution of tilts from the vertical and random azimuth,
such as aligned ice or falling drops. They give the full phase matrix, the extinction (K11, K12, K34) and the
absorption (I, Q) on ARTS's grids, as native ARO habits. The azimuthal average and the absorption integral use the
number of nodes that the band limits of the truncated T-matrix make exact (derived in concept.tmatrix). They are
checked against Mie, against random orientation (uniform tilts) and against twice the nodes.

Generalized Chebyshev particles (new). Raindrop shapes r0 (1 + c0 + sum c_n cos n theta) with the user's
coefficients. The original fixed one drop's coefficients in its code and could not compute random orientation for
them.

Bugs fixed

• Random orientation was wrong for large particles. The original's Clebsch-Gordan recursion loses accuracy as the
  truncation n1 grows: errors of 1e-14 at n1 = 14, 1.4e-6 at n1 = 100, and about 1e8 at n1 = 240, the limit ARTS had
  raised it to. The exact coefficients remove this.
• Bulk ARO properties of spectral TRO particles were approximate. A Legendre series was put on a grid of 2L + 2
  angles and then interpolated linearly in the scattering angle, about 1.4 % off for a Rayleigh-size T-matrix habit.
  The series is now evaluated at each direction pair's exact scattering angle. This affects every spectral TRO habit
  used by an ARO consumer (MC, RT4, VDISORT).
• Polarization bases near the poles and near forward/backward scattering. rtepack treated two directions as
  parallel below a cross product of 1e-12, about 4500 eps:
  • rays within 1e-12 rad of the zenith or nadir took a fixed basis instead of their own meridian;
  • near-forward and near-backward pairs used the wrong scattering plane, about 0.99 off at 1e-13 rad;
  • the plane of incidence at normal incidence had the same problem.
    The bound is now 16 eps. The plane is also projected back across the incident ray; without that projection,
    rounding distorted the phase matrix by about (eps / |k_in x k_out|)^2 even for a regular F.
• T-matrix amplitude near the particle axis. The original moved angles 1e-7 rad off the singular directions and
  used the pole limits within 1.4e-5 rad of the axis, an error of about 2e-5 there. Its arccos angle lost sqrt(eps).
  Angles now come from atan2 of the components, and the limits are used only on the axis itself.
• Power-law size distribution. The radius root was searched in (1e-5, a - 1e-5), and an endpoint was returned when
  the root lay outside, which gave the wrong distribution for large variances (b = 10 has r1 ~ 3e-9 a). It is now
  found on the whole interval, and a root that double precision cannot resolve is an error.
• Silent non-convergence. The original used an unconverged T-matrix when its quadrature reached its array limits,
  and it ended the expansion series at a fixed threshold whatever the requested accuracy. The series is now
  complete, and non-convergence is an error, also for a T-matrix that is not finite. Such matrices were handed to
  LAPACK, where the outcome depended on the BLAS library.
• matrix_exp was wrong. Its Padé denominator had the wrong sign for even powers, 7 % off at |A| ~ 1, because
  nonstd::pow(-1, k) returned -1 for every k. nonstd::pow returned -|x|^v for negative x. It is now defined for
  x >= 0 only, and integer powers use exact forms. No production code called matrix_exp.
• Smaller fixes:
  • ExtinctionMatrixData<ARO>::extract_stokes_coeffs had never compiled;
  • the Fourier-mode extinction and absorption of bulk ARO properties could not be returned to Python.

For owners of legacy scattering data

ARTS2-format data (and probably the ARTS-SSDB) appear to use Mishchenko's convention: the regenerated MC fixture
matches the ARTS2 XML only in that convention. ARTS3 reads them without conversion, so their F34, and the U rows and
columns of their ARO phase matrices, have the opposite sign to ARTS's convention. This PR does not change the legacy
readers. Each database's convention, including how it defines the azimuth, should be checked and converted on read.

Evidence

The T-matrix is checked against:

• the independent Mie code, in all 16 elements of the laboratory-frame phase matrix;
• every number printed by the original programs (double and extended precision);
• the original Fortran's outputs over 124 cases, recorded before its removal (spheroids, cylinders, Chebyshev and
  generalized Chebyshev particles, all size distributions);
• the analytic random-orientation average.

[riclarsson]

#1194 is affected by this merge. A design change to make muelmat{} return an empty instead of diagonal matrix was changed so that the interpolation routine could be made to work. There are a few places in the surface scattering PR that made use of the previous default behavior.