Message Boards Message Boards

ContourPlot in Fast Spherical Harmonics on radiative transfer

Posted 8 years ago

Hi, This code was written by Markus Roelling, 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}]]
  }
 ]

Then, I'd like 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, omega]/meanIntensity[0], {tau, 0, 10}, {omega, 0, 1}, PlotRange -> {{0, 10}, All}]

Any help will be very appreciable.

POSTED BY: Sungwoo Yang
4 Replies

Your function B wants omega as a number. What you feed it is a list omega = Range[0.1,0.9,0.2]. Perhaps you should specify the action of B when omega is a list.

POSTED BY: Gianluca Gorni

This is where I am confused. It seems that meanIntensity is just a function of one variable and intensity is a function of two variables. I guess what you are looking for is similar to the following simplified version:

f[x_]:=x+a
g[x_,y_]:=x+y+b

both f and g are functionals depending on the parameters a and b. Your plot should look like a family of curves or surfaces? Am I right?

POSTED BY: Shenghui Yang

In your case before the omegafunctlion is executed, meanIntensity would not be defined. A simplified experiment like below:

In[25]:= f[x_]:=a=3
In[26]:= Head[a]
Out[26]= Symbol
In[27]:= f[1]
Out[27]= 3
In[28]:= Head[a]
Out[28]= Integer

You can plot with the following code:

ContourPlot[ Evaluate[meanIntensity[tau, omega]/meanIntensity[0]], {tau, 0, 10}, {omega, 0, 1}, PlotRange -> {{0, 10}, All}]

sample1

POSTED BY: Shenghui Yang
Posted 8 years ago

Hi Shenghui, Thank you so much for your reply. I was able to get the plot you showed above, but it is not a function of omega. Here is what I typed as you suggested.

omega = 0.1;
{meanIntensity, intensity} = 
  createSphericalHarmonics[0.5, omega, 1., 5, 19];
ContourPlot[
 Evaluate[meanIntensity[tau, omega]/meanIntensity[0]], {tau, 0, 
  5}, {omega, 0, 1}, PlotRange -> {{0, 5}, All}]

But what I am looking for is plot which is function of omega. For example, if I edit below, the plot is not working.

omega = Range[0.1,0.9,0.2];
{meanIntensity, intensity} = 
  createSphericalHarmonics[0.5, omega, 1., 5, 19];
ContourPlot[
 Evaluate[meanIntensity[tau, omega]/meanIntensity[0]], {tau, 0, 
  5}, {omega, 0, 1}, PlotRange -> {{0, 5}, All}]

I appreciate your help.

POSTED BY: Sungwoo Yang
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract