Message Boards Message Boards

GROUPS:

ContourPlot in Fast Spherical Harmonics on radiative transfer

Posted 6 years ago
5736 Views
|
4 Replies
|
2 Total Likes
|

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.

4 Replies

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 6 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.

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?

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.

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