2
$\begingroup$

This code was written by Markus Roellig, and I'd like to mention that all credits should go to him. Please, see this link to open his original writing.

SphericalBesselI[0, 0] := SphericalBesselI[0, 0] = 1;
SphericalBesselI[l_, z_] := Sqrt[π/(2 z)] BesselI[l + 1/2, z];
σ[l_][g_] := g^l
h[l_][g_, ω_] := h[l][g, ω] = (2 l + 1) (1 - ω σ[l][g]);
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];
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]];
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]);
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}];

createSphericalHarmonics[g_, ω_, I0_, tauMax_, L_?OddQ] := 
 Module[
  {M = (L + 1)/2, angles, eigenvalues, kEigenValues, AList},
angles = -SetPrecision[getAngles[M], Infinity];
eigenvalues = 
    SetPrecision[calculateEigenvalues[L][g, ω], Infinity];
 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]];

 {
  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}]]
  }
 ]

By modifying the above code, Fast Spherical Harmonics radiative transfer, I am trying to plot intensity field, function of albedo (omega) and optical depth (tau) at theta angle = 0, using ContourPlot, but no luck yet. I can't figure out what I did wrong. Here is what I edit at the end of what your wrote.

omegafunction[omega_] := {meanIntensity, intensity} = createSphericalHarmonics[0.5, omega, 1., 10, 19];
ContourPlot[meanIntensity[tau]/meanIntensity[0], {tau, 0, 10}, {omega, 0, 1}, PlotRange -> {{0, 10}, All}]

Any help will be very appreciable. Thank you.

$\endgroup$
  • 3
    $\begingroup$ You should share the definition of createSphericalHarmonics and other constants too, so people can try to reproduce your problem. $\endgroup$ – MarcoB Feb 15 '16 at 23:49
  • $\begingroup$ Hi MarcoB, Thank you so much for your comment. I edited my question as you suggest. Hope I can get some comments from other experts. Thank you. $\endgroup$ – SungwooY Feb 16 '16 at 1:07
2
$\begingroup$

I think there is a typo: do you want to divide intensity by meanintensity?

Anyways, one way to do what you want is to tabulate your results and use ListContourPlot. Note that the computation fails for omega=1:

data = ParallelTable[Table[
{meanIntensity, intensity} = createSphericalHarmonics[0.5, omega, 1., 10, 19];
{tau,omega,Evaluate[intensity[tau, omega]/meanIntensity[0]]}, 
 {tau, 0, 1, .1}], {omega, 0, .9, 0.05}];
ListContourPlot[Flatten[data, 1]]

enter image description here

$\endgroup$
  • $\begingroup$ Dear Markus Roelling, It's an honor to receive a feedback from you who wrote the original code. Yes, it was a type and I fixed it.. This is the one I've been looking for. I appreciate your help. $\endgroup$ – SungwooY Feb 18 '16 at 3:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.