26
$\begingroup$

Context

A couple of years ago I posted this question for an efficient code to generate an n-D Gaussian random field (sometimes called processes in other fields of research), which has applications in cosmology.

For instance, the function GaussianRandomField would work as follows in 2D

 u = GaussianRandomField[] //Chop// GaussianFilter[#, 1] &; Image[u] // ImageAdjust

a 2D GRF

Question

I am interested in an efficient way to proceed on the sphere.

Motivation

I would eventually like to make maps like this map of the Cosmic Microwave Background seen by the Planck satellite.

Mathematica graphics

using this code.

Possible difficulties

It most likely involves fast Spherical Harmonic transforms of relatively high order? One possibility might be to link to the Healpix library, but hopefully it would be an overkill.

Feeble attempt

Extracting the transformation code from the above mentioned page

invmollweide[{x_, y_}] := With[{theta = ArcSin[y]}, {Pi x/(2 Cos[theta]), ArcSin[(2 theta + Sin[2 theta])/Pi]}]

I can generate a map:

 Clear[alms]; lmax = 48; Do[alms[l, m] = (Random[NormalDistribution[]] + I Random[NormalDistribution[]]) /Sqrt[(l + 2) (l + 1)]; alms[l, -m] = (-1)^m Conjugate@alms[l, m];, {l, 0, lmax}, {m, 0, l}]; Do[alms[l, 0] = (Random[NormalDistribution[]])/Sqrt[(l + 2) (l + 1)]; , {l, 0, lmax}]; Clear[field]; field[θ_, ϕ_] = Sum[alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l}]; fieldN = Compile[{θ, ϕ}, field[θ, ϕ] // Evaluate]; dat = ParallelTable[fieldN[θ, ϕ], {θ, 0, Pi, Pi/128.}, {ϕ, 0.,2 Pi, 2 Pi/256.}] // N;

and then plot it: 

im = Re[dat] // Image // ImageAdjust// Colorize[#, ColorFunction -> "LightTemperatureMap"] &; mol=ImageTransformation[im, invmollweide, DataRange -> {{-Pi, Pi}, {-Pi/2, Pi/2}}, PlotRange -> {{-2, 2}, {-1, 1}},Padding-> White]

Mathematica graphics

But

  1. The production of the map is rather slow…

  2. The code breaks down e.g. at lmax=64

Mathematica graphics

probably because of the accuracy of the spherical harmonics (whereas in astronomy people routinely use lmax = 4096 or more).

It's a pity because Mathematica allows for some cool visualization. 

SphericalPlot3D[1, {u, 0, Pi}, {v, 0, 2 Pi}, PlotPoints -> 50, MaxRecursion -> 0, Mesh -> True, TextureCoordinateFunction -> ({#5, 1 - #4} &), PlotStyle -> Directive[Texture[im], Specularity[White, 50]], Lighting -> "Neutral", Boxed -> False, Axes -> False]

Mathematica graphics Mathematica graphics

See also this

Possible Extension: random vector field on sphere

One can also define vector fields on the sphere as follows

fc[phi_] := Block[{theta}, If[Abs[phi] == Pi/2, phi, theta /. FindRoot[2 theta + Sin[2 theta] == Pi Sin[phi], {theta, phi}]]]; cart[{lambda_, phi_}] := With[{theta = fc[phi]}, {2/Pi*lambda Cos[theta], Sin[theta]}]

Then, if we define VectorSphericalHarmonicV

Clear[ϵ];(*Polarization vector*) ϵ[λ_] = Switch[λ, -1, {1, -I, 0}/Sqrt[2], 0, {0, 0, 1}, 1, {1, I, 0}/Sqrt[2]]; Clear[VectorSphericalHarmonicV]; VectorSphericalHarmonicV[ℓ_, J_, M_, θ_, ϕ_] /; J >= 0 && ℓ >= 0 && Abs[J - ℓ] <= 1 && Abs[M] <= J := Sum[If[Abs[M - λ] <= ℓ, ClebschGordan[{ℓ, M - λ}, {1, λ}, {J, M}], 0]* SphericalHarmonicY[ℓ, M - λ, θ, ϕ]*ϵ[λ], {λ, -1, 1}] pp = Sum[(Random[NormalDistribution[]] + I Random[NormalDistribution[]]) Rest@ VectorSphericalHarmonicV[i, i - 1, i - 1, θ, ϕ], {i,2, 8}];

We can write 

 gr2 = StreamPlot[pp//Re, {ϕ, -Pi, Pi}, {θ, -Pi/2, Pi/2}, AspectRatio -> 1/2, Frame -> False, StreamColorFunction -> "ThermometerColors", StreamPoints -> 250]; gr2 = gr2 /. Arrow[pts_] :> Arrow[(cart /@ pts)] /. Point[pts_] :> Point[cart[pts]] // Show[#, PlotRange -> {{-2, 2}, {-1, 1}}] &;

and

Graphics[{Inset[mol, {-2, -1}, {0, 0}, {4, 2}], First[gr2]}, PlotRange -> {{-2, 2}, {-1, 1}}]

Mathematica graphics

I guess once we have a fast spherical harmonic code on the sphere it will be trivial to generalize it to vector fields.

Update

This package might be of relevance. I quote

Wolfram Mathematica is a powerful and convenient software package used by many cosmologists everywhere. However, since it is not always the most efficient for low-level computing, most popular algorithms for numerical computations in cosmology have been written in C or Fortran, since these languages are typically much better suited for the task at hand. This package bundles the functionality of some of these algorithms in a Mathematica package, which makes them easier to use and avoids the need to learn C or Fortran.

Improve this question
$\endgroup$
4
  • $\begingroup$ @shrx thanks for the bounty! $\endgroup$ Commented Sep 15, 2015 at 20:27
  • $\begingroup$ are you sufficiently satisfied with either of the answers so I can award the bounty? $\endgroup$ Commented Sep 21, 2015 at 8:54
  • $\begingroup$ @shrx I guess I would rather leave it to you to decide? $\endgroup$ Commented Sep 21, 2015 at 9:20
  • $\begingroup$ To be honest neither answers your question sufficiently in my opinion, so I wouldn't award it. $\endgroup$ Commented Sep 21, 2015 at 10:21

2 Answers 2

Reset to default
6
+25
$\begingroup$

GaussianRandomField is a special case. More generally, what is required is fast code for the (inverse) Spherical Harmonic Transform (SHT), which will work for any coefficients $a_{l,m}$. SHTns is a high performance library for Spherical Harmonic Transform written in C and so should be straightforward to link in using MathLink. It would be very nice if this code was built-in to a future version of Mathematica.

To see why fast code is required, for $l_{max} = 2^{12}$ direct computation requires evaluation of around $2^{24}$ spherical harmonics on your grid of $\theta$ and $\phi$ values (another factor of $2^{15}$). There are recursive tricks, of course, but fast methods rely on the FFT.

Regarding your example code, I do not see why you have (-1)^m in your definition for alms[l, -m]. To obtain an explicitly real function you should just have alms[l, -m] = Conjugate@alms[l, m].

Further to your visualisation on the sphere, you can also select your view-centre, e.g. Galactic Coordinates.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for your indications. Regarding the code, if ` alms[l, -m] = Conjugate@alms[l, m].` then the map appears pi periodic horizontaly. $\endgroup$ Commented Sep 18, 2015 at 19:18
  • $\begingroup$ One could also use the library libsharp $\endgroup$ Commented Sep 21, 2015 at 7:41
  • 1
    $\begingroup$ As stated above, a possible solution is to link Mathematica to a high performance external SHT library with MathLink. Among them there is the C GNU-GPL library libsharp (github.com/dagss/libsharp, aanda.org/10.1051/0004-6361/201321494) which has the advantage of supporting several different sphere pixelisations (Gauss-Legendre, ECP, HEALPix, ...), and SHT of arbitrary spins. $\endgroup$ Commented Sep 21, 2015 at 9:34
5
$\begingroup$

A part from efficiency, I noticed that with the definition below, lmax >= 64 works:

Clear[field]; field[θ_, ϕ_] := Chop@ Total[ Table[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ] , 2 ]; nn = 4.; dat = ParallelTable[ field[θ, ϕ], {θ, 0, Pi, Pi/nn}, {ϕ, 0., 2 Pi, 2 Pi/nn/2} ]; Round[Re[dat], 0.01]

{ {0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51},

{0.02, -1.54, -1.48, -0.47, -2.29, 0.61, 2.05, 1.82, 0.02},

{0.24, -1.61, 0.44, 0.32, -0.55, 0.65, -1.12, -0.08, 0.24},

{-0.34, 0.63, -0.54, 2.08, -0.72, 1.09, -2.09, -1.3, -0.34},

{0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} }

which looks better and produces a reasonable plot.

So, I tested also:

Clear[field]; field[θ_, ϕ_] := Sum[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ];

and this too seems to me to work.

Consequently, it appears to me that Compile over Sum is doing some inappropriate manipulation.

EDIT

I could perhaps track down the problem a little better: it seems to be the SphericalHarmonicY at angles θ below Pi/4 and above Pi - Pi/4, which finally provide extremely high values in dat and hides the details of the other values (at least for lmax = 64).

For example the table:

dat = ParallelTable[ fieldN[θ, ϕ], {θ, Pi/4., Pi - Pi/4., Pi/128.}, {ϕ, 0, 2 Pi, 2 Pi/256.} ];

seems to me to produce a reasonable chart.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.