Fast Spherical Harmonics radiative transfer

Question

This is a rather specific question and I apologize for spamming you with some lengthy code. But it could be interesting for some reader and maybe you can help out, so please bear with me.

I am using the Spherical Harmonics Method to solve a radiative transfer problem in spherical symmetry, namely the penetration of diffuse UV radiation into interstellar clouds. It uses an expansion in Legendre polynomials to approximate the depth and angular dependence of the problem.

I have a working solution and I am interested in a some additional input:

1. Is there a more elegant/simpler way to set up solutions for the intensity for varying values of $\omega$ (scattering albedo) and g (mean scattering angle $\langle \cos(\theta)\rangle$)
2. My solution is monochromatic and still relatively slow. For real applications one needs to perform the same calculations on a wavelength grid numerous times. How to speed up the solution? Wouldn't it be neat if Mathematica could do the radiative transfer efficiently?
3. For higher approximation order, solving the Eigenvalue problem and getting the list of $A_m$ from the boundary conditions becomes extremely slow. Any ideas on how to improve?

Here is what I got so far. First some definitions:

(* definition of the modified spherical Bessel function of the first kind *)
SphericalBesselI[0, 0] := SphericalBesselI[0, 0] = 1;
SphericalBesselI[l_, z_] := Sqrt[π/(2 z)] BesselI[l + 1/2, z];
(* expansion coefficients for the Henyey-Greenstein phase function *)
σ[l_][g_] := g^l
(* diagonal matrix elements *)
h[l_][g_, ω_] := h[l][g, ω] = (2 l + 1) (1 - ω σ[l][g]);
(* low-order determinants *)
F[0][g_, ω_][k_] := F[0][g, ω][k] = -h[0][g, ω];
F[1][g_, ω_][k_] := F[1][g, ω][k] = h[0][g, ω] h[1][g, ω] - k^2;
F[l_][g_, ω_][k_] := F[l][g, ω][k] = -h[l][g, ω] F[l - 1][g, ω][k] -
  l^2 k^2 F[l - 2][g, ω][k];
(* helper function to access the i-th Eigenvalue k *)
kEV::mmatch = "Inconsistent list of Eigenvalues.";
kEV[list_List][m_Integer] := Module[{pos, neg},
   pos = Select[Sort@list, Positive];
   neg = Select[Sort@list, Negative];
   If[(Length[pos] == Length[neg]) && m <= Length[pos],
     If[Positive[m], pos[[m]], neg[[m]]],
     Message[kEV::mmatch]; $Failed]];
(* calculate the Eigenvalues *)
calculateEigenvalues[L_][g_, ω_] := 
   NSolve[F[L][g, ω][k] == 0, k][[All, 1, 2]];
R[0, m_][g_, ω_][k_] := R[0, m][g, ω][k] = 1;
R[1, m_][g_, ω_][k_] := R[1, m][g, ω][k] = (1 - ω)/k[m];
R[l_Integer, m_Integer][g_, ω_][k_] /; m < 0 := R[l, m][g, ω][k] = (-1)^l R[l, -m][g, ω][k]
R[l_Integer, m_][g_, ω_][k_] :=  R[l, m][g, ω][k] = 
     1/(l k[m]) (h[l - 1][g, ω] R[l - 1, m][g, ω][
      k] - (l - 1) k[m] R[l - 2, m][g, ω][k]);
(* Mark's boundary condition *)
getAngles[M_Integer] := 
 Sort@Select[
 List @@ (NRoots[LegendreP[2 M, x] == 0, x] /. Equal[_, x_] :> x), 
 Negative] 
B[i_Integer, m_Integer][g_?NumberQ, ω_?NumberQ][Lmax_, 
   taumax_?NumberQ, k_, angles_] := 
   Sum[(2 l + 1) R[l, m][g, ω][k] LegendreP[l, 
   angles[[i]]] SphericalBesselI[l, k[m] taumax], {l, 0, Lmax}];

Now the procedure to set up function to calculate the mean intensity (averaged over the full solid angle) as function of the optical depth $\tau$ and the specific intensity as function of optical depth $\tau$ and angle $\mu=\cos(\theta)$.

createSphericalHarmonics[g_, ω_, I0_, tauMax_, L_?OddQ] := 
 Module[
  {M = (L + 1)/2, angles, eigenvalues, kEigenValues, AList},
  (* angles for the boundary conditions from Mark's condition *)
  angles = -SetPrecision[getAngles[M], Infinity];
  (* calculate Eigenvalues from the recurrence relations F*)
  eigenvalues = 
    SetPrecision[calculateEigenvalues[L][g, ω], Infinity];
  (* calculate the A_k from the boundary conditions, i.e. 
     that the incoming intensity at the edge of the sphere and at
     angle_i equals I_ 0*)
  (* assuming isotropic illumination gives ConstantArray[I0, Mmax]. 
     Alternatively we could define an I_ 0[mu] as declination dependent I_ 0*)
  AList = LinearSolve[
     N[SetPrecision[
        Table[B[i, m][g, ω][L, tauMax, kEV[eigenvalues], angles], 
              {i, 1, M}, {m, 1, M}]
      , Infinity], 30], 
     N[SetPrecision[ConstantArray[I0, M], Infinity], 30]];

 (* the N[#, 30] should prevent the LinearSolve::luc error but it doesn't *)
 (* output are two InterpolationFunctions: {mean intensity as function of tau, intensity as function of angle and optical depth tau} *)
 {
  Function[tau, 
   Sum[AList[[m]] R[0, m][g, ω][
   kEV[eigenvalues]] SphericalBesselI[
   0, (tauMax - tau) kEV[eigenvalues][m]], {m, 1, M}]],
  Function[{tau, mu},
   Sum[(2 l + 1) LegendreP[l, mu] Sum[
   AList[[m]] R[l, m][g, ω][
     kEV[eigenvalues]] SphericalBesselI[
     l, (tauMax - tau) kEV[eigenvalues][m]], {m, 1, M}], {l, 0, L}]]
  }
 ]

To demonstrate the solution, let's assume $g=0.5$, albedo $\omega=0.6$, approximation order of the Legendre polynomials $L=19$ and incoming intensity $I_0=1$, and optical depth to cloud center $\tau_c=10$.

{meanIntensity, intensity} = createSphericalHarmonics[0.5, 0.6, 1., 10, 19];

PolarPlot[{1, intensity[0, Cos[mu]], intensity[1, Cos[mu]]}, {mu, 0, π/2}, 
     Epilog -> {Dashed, (Line[{{0, 0}, {Cos[ArcCos @ #], 
     Sin[ArcCos @ #]}}] & /@ (-getAngles[10]))}]

The wavy curve is the approximated intensity as function of angle $\mu$ at the surface of the cloud, the dashed lines show the angles that were used to specify the boundary conditions. The olive curve shows the angular dependence of the intensity at $\tau=1$.

The angular averaged, mean intensity:

Plot[Log10@(meanIntensity[tau]/meanIntensity[0]), {tau, 0, 10}, 
     PlotRange -> {{0, 10}, All}]

I am pretty sure there are plenty of bugs still in the code, but they are not the point of this post, but rather the performance and coding style.