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
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
.
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]
But
The production of the map is rather slow…
The code breaks down e.g. at
lmax=64
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]
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}}]
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.