4
$\begingroup$

I am plotting an anisotropy distribution function as color coding on the surface of a unit sphere with this code:

PAD[β_, θ_] :=
 1/(4 π) (1 + β LegendreP[2, Cos[θ]])

SphericalPlot3D[1, {θ, 0, π}, {ϕ, 0, 2 π}, 
 ColorFunction ->
  Function[{x, y, z, θ, ϕ, r},
   ColorData["Rainbow"][
    Rescale[PAD[β = 1, θ], {0, 0.5}, {0, π}]]], 
 Mesh -> False, Boxed -> False, Axes -> True, AxesLabel -> {x, y, z}, 
 ColorFunctionScaling -> False, PlotPoints -> 100,
 SphericalRegion -> True]

enter image description here

Thanks to this solution. Next, I would like to calculate the integral along the y axis and plot the result as a density map in the same color coding. It should look similar to this image (the outer ring, disregard the inner rings).

Embarrasingly, I fail to even find a proper starting point in this problem. I have to calculate the projection of my function on a sphere, scale the values according to the color coding of the 3D plot, transform to cartesian coordinates, integrate, density plot. While the last steps are simple, my mind is stuck on how to begin with the first two.

Update: I un-stuck my mind thanks to the discussion with Alexei and came to the following approach: I need a function for the radial part, and chose a Gaussian function

Gaussian[r_, r0_, s_] := Exp[-(r - r0)^2/(2 s^2)]

and then made a DensityPlot of a slice through the distribution with

DensityPlot[
 Gaussian[Sqrt[x^2 + y^2 + z^2], 1, 0.05] 
 PAD[β = 1, ArcTan[z, Sqrt[x^2 + y^2 + z^2]]] /. y -> 0,
 {x, -1.2, 1.2}, {z, -1.2, 1.2}, PlotPoints -> 150, 
 ColorFunction -> ColorData["Rainbow"], PlotRange -> All]

enter image description here

which works nicely. To get the projection along the y axis, I wrapped the plotted functions into NSum as in

NSum[<func>, {y, -1.2, 1.2, .25}]

and repeated the DensityPlot. This is taking forever and the results depend obviously a lot on the chosen increment of NSum. For larger values, I do not get a smooth distribution inside the ring, and smaller values take hours to complete. Does anybody have a better idea?

|improve this question|||||
$\endgroup$
  • $\begingroup$ What integral do you want to to find? Nothing is written about it. There should be some expression under the integral, should it not? $\endgroup$ – Alexei Boulbitch Oct 21 '15 at 9:44
  • $\begingroup$ Hi Alexei, thank you for your question. $\endgroup$ – Ole Oct 21 '15 at 12:42
  • $\begingroup$ I accidentally pressed enter, sorry. Next part of my answer: I am looking for the projection of the sphere onto the (x,z) plane, located e.g. at y=1 (in a graphical sense). In a numerical solution, I want to add up a lot of (x,z) slices through the sphere. Analytically this would be the projection $P(x,z) = \int_{-\infty}^\infty I(x,y,z) dy$, where I(x,y,z) is the plotted sphere with values corresponding to the color coding, transformed to cartesian coordinates. I want this for illustrative purposes only. $\endgroup$ – Ole Oct 21 '15 at 12:50
  • $\begingroup$ Hi, Ole, first there is the button "edit" below your question. With its help you may add and change your question any time. Second, your words: "I(x,y,z) is the plotted sphere with values corresponding to the color coding" are rather obscure. To integrate something you need to have a mathematical expression. In Mma the color coding may be done in several ways. For example, using the Hue function. Check Menu/Help/WolframDocumentation/ColorFunction and have a look at the very first example, or some other examples showing different ways. $\endgroup$ – Alexei Boulbitch Oct 21 '15 at 13:13
  • $\begingroup$ Continuation: Or do you, may be, have in mind the function Function[{x, y, z, \[Theta], \[Phi], r}, ColorData["Rainbow"][ Rescale[PAD[1, \[Theta]], {0, 0.5}, {0, \[Pi]}]]] defined in your question? $\endgroup$ – Alexei Boulbitch Oct 21 '15 at 13:14
1
$\begingroup$

Given the functions you have defined, I think you simply want to calculate the integral

Gaussian[Sqrt[x^2 + y^2 + z^2], 1, 0.05] PAD[β = 1, 
  ArcCos[z/Sqrt[x^2 + y^2 + z^2]]]

But I'm too impatient to wait on that to return an answer. So I just do it numerically. I define a function of x and y that takes the integral numerically,

projfunc[x_,z_]:=projfunc[x,z]=NIntegrate[Gaussian[Sqrt[x^2+y^2+z^2],1,0.05] 
                               PAD[β=1,ArcCos[z/Sqrt[x^2+y^2+z^2]]],{y,-∞,∞}]

Then I create a list to plot, which goes very quickly,

listdata = 
  Table[projfunc[x, z], {z, -1.2, 1.2, .05}, {x, -1.2, 1.2, .05}];

And I plot the data,

ListDensityPlot[listdata, InterpolationOrder -> 2, 
 ColorFunction -> ParulaCM, PlotLegends -> Automatic, 
 DataRange -> 1.2 {{-1, 1}, {-1, 1}}]

enter image description here

Of course, this is using a custom color function, so if you don't have that, just take out that option.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ This is exactly what I was looking for, thank you very much! Actually I looked at NIntegrate, but thought it could only integrate over all coordinates to yield one value. Thats the reason I went for NSum. I should have known better :-) Greetings from Kiel, did we unknowingly meet at the Femto12 in July? $\endgroup$ – Ole Oct 21 '15 at 15:45
  • $\begingroup$ @Ole Hallo, perhaps we did meet at the conference - there were a lot of really nice talks there. If you have any other issues plotting VMI data using Mathematica I'd be happy to help. $\endgroup$ – Jason B. Oct 28 '15 at 9:10

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.