# 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
Sort By:
Posted 6 years ago
 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}] 
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?