If the single-scattering albedo equals $0$, the coefficient matrix will not be diagonalizable. It is still possible to express the coefficient matrix $A = G J G^{-1}$, but $J$ will be a Jordan matrix and we will need to use generalized eigenvectors to solve the system of ODEs. Currently, users may use single-scattering albedos that are close to $1$ to approximate conservative scattering and PythonicDISORT will warn if a single-scattering albedo is too close to $1$.
We would like PythonicDISORT to instead use the special case solution if the single-scattering albedo is too close to $1$. This may render the joblib parallelization obsolete.
