The objective of this notebook is the simulation of second order isotropic random fields defined on the three dimensional sphere using the spectral algorithm proposed in Ch. Lantuéjoul et al. (2019).
Any point on the sphere \(\mathbb{S}_2\) can be defined by two spherical coordinates, its colatitude \(\theta \in [0, \pi]\) and its longitude \(\varphi \in [0, 2\pi]\). If we consider the points of the sphere as points of \(\mathbb{R}^3\) the relationship between its Cartesian coordinates and its spherical coordinates is: \[ \left\{ \begin{array}{ccl} x & = & \sin (\theta) \, \cos (\varphi) \\ y & = & \sin (\theta) \, \sin (\varphi) \\ z & = & \cos (\theta) \end{array} \right. \]
The ingredients of the simulation method are:
\[ C(t) = \sum_{n = 0}^{+\infty} a_n P_n(\cos t) \text{ with } t \in [0, \pi]. \]
The spectral component \(Z(s) = 2\sqrt{2\pi} \, \Re(Y_{N, K}(s)\, e^{i\Phi})\) is a standardized random field with covariance function \(C\) with,
In the previous equations,
\(P_n(t)\), \(-1 \leq t \leq 1\) are the Legendre polynomials,
\(Y_{n, k}(\theta, \varphi)\) are the spherical harmonics, acting on the sphere exactly like trigonometric functions on the unit circle.
The Legendre polynomials are defined by the Rodrigues’formula:
\[ P_n(x) = \frac{(-1)^n}{2^n n!} \frac{d^n}{dx^n}(1-x^2)^n. \] They can be computed using the recursion formula:
\[ (n+1) P_{n+1}(x) = (2n+1)\,x\, P_{n}(x) - n P_{n-1}(x), \] with the initial values \(P_0(x) = 1\) and \(P_1(x) = x\).
The Legendre polynomials are orthogonal but not normalized:
\[ \int_{-1}^{+1} P_m(u)P_n(u)du = \frac{2}{2n+1}\delta_{m,n}. \]
For numerical stability, we use normalized polynomials \(\tilde{P}_l(x) = P_l(x)\times \sqrt{2n+1}\). In this case, the recursion formula is: \[ \tilde{P}_{n+1}(x) = \frac{\sqrt{(2n+3)(2n+1)}}{n+1}\,x\, \tilde{P}_{n}(x) - \frac{n}{n+1}\sqrt{\frac{2n+3}{2n-1}} \tilde{P}_{n-1}(x), \] with the initial values \(\tilde{P}_0(x) = 1\) and \(\tilde{P}_1(x) = x\sqrt{3}\).
The associated Legendre functions denoted \(P_l^m(x)\) with \(x \in [-1,1]\) and \(-l \leq m \leq l\) are defined by the equation: \[ P_l^m(x) = (-1)^m \, (1 - x^2)^{m/2} \frac{d^m}{dx^m}(P_l(x)), \] where \(P_l\) is the Legendre polynomial of degree \(l\) (\(P_l = P_l^0\)), and \(m\) its order.
We consider also normalized Legendre functions as,
\[ \tilde{P}_l^m(x) = P_l^m(x) \times \sqrt{(2l+1)\frac{(l-m)!}{(l+m)!}}, \] which can be compute using the two recursion formulae:
\[ \tilde{P}_{l+1}^{l+1}(x) = -\sqrt{\frac{2l+3}{2l+2}} \, (1-x^2)^{1/2} \, \tilde{P}_{l}^{l}(x) \] and
\[ \tilde{P}_{l+1}^{m}(x) = \sqrt{\frac{(2l+3)(2l+1)}{(l+m+1)(l-m+1)}} \, x \, \tilde{P}_{l}^{m}(x) - \sqrt{\frac{(l+m)(l-m)(2l+3)}{(l+m+1)(l-m+1)(2l-1)}} \, \tilde{P}_{l-1}^{m}(x) \] with the initial conditions, \(\tilde{P}_0^0 = 1\) and \(\tilde{P}_1^0 = x \sqrt{3}\).
Note: The recursion for the Legendre functions boils down to the recursion for the Legendre polynomials with \(l = 0\).
The spherical harmonics act on the sphere \(\mathbb{S}_2\) exactly like trigonometric functions on the unit circle, they are defined using the normalized Legendre functions:
\[ Y_{n, k}(\theta, \varphi) = \frac{1}{\sqrt{4\pi}} \tilde{P}_n^l (\cos \theta) e^{ik\varphi}. \] In addition, the value for negative order can is (\(k\geq0\)):
\[ Y_{n, -k}(\theta, \varphi) = (-1)^k \, \bar{Y}_{n, k}(\theta, \varphi). \]
The simulation of the isotropic random function on the sphere with a spectrum \(a_n\) is:
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display## Warning: 'rgl.init' failed, will use the null device.
## See '?rgl.useNULL' for ways to avoid this warning.defineDefaultSpace(ESpaceType_SN(),param=6371)## NULLnx = c(360, 180)
dx = nx/(nx-1)*pi/180
x0 = c(0,0)
grd = DbGrid_create(nx = nx, x0 = x0, dx = dx)
err = grd$setName("x1", "phi")
err = grd$setName("x2", "theta")Initialization of S2 for 3d viewer
viz_S2_in_3D = IniView_S2_in_3D(grd)
# display of the real part of the spherical harmonics using rgl
val = ut_sphericalHarmonicVec(n = 15, k = 5, theta = grd["theta"], phi = grd["phi"])
err = viz_S2_in_3D$display(val = Re(val), ncol = 50, palette = "rainbow")
rglwidget()