|
1 |
| -############# |
| 1 | +# ## Spherical harmonic addition theorem |
2 | 2 | # This example confirms numerically that
|
| 3 | +# ```math |
| 4 | +# \frac{P_4(z\cdot y) - P_4(x\cdot y)}{z\cdot y - x\cdot y}, |
| 5 | +# ``` |
3 | 6 | #
|
4 |
| -# [P₄(z⋅y) - P₄(x⋅y)]/(z⋅y - x⋅y), |
5 |
| -# |
6 |
| -# is actually a degree-3 polynomial on 𝕊², where P₄ is the degree-4 |
7 |
| -# Legendre polynomial, and x,y,z ∈ 𝕊². |
8 |
| -# To verify, we sample the function on a 5×9 equiangular grid |
| 7 | +# is actually a degree-$3$ polynomial on $\mathbb{S}^2$, where $P_4$ is the degree-$4$ |
| 8 | +# Legendre polynomial, and $x,y,z \in \mathbb{S}^2$. |
| 9 | +# To verify, we sample the function on a $5\times9$ equiangular grid |
9 | 10 | # defined by:
|
10 |
| -# |
11 |
| -# θₙ = (n+1/2)π/N, for 0 ≤ n < N, and |
12 |
| -# |
13 |
| -# φₘ = 2π m/M, for 0 ≤ m < M; |
14 |
| -# |
| 11 | +# ```math |
| 12 | +# \theta_n = (n+\frac{1}{2})\pi/N,\quad{\rm for}\quad 0\le n < N,\quad{\rm and} |
| 13 | +# ``` |
| 14 | +# ```math |
| 15 | +# \varphi_m = 2\pi m/M,\quad{\rm for}\quad 0\le m < M; |
| 16 | +# ``` |
15 | 17 | # we convert the function samples to Fourier coefficients using
|
16 | 18 | # `plan_sph_analysis`; and finally, we transform
|
17 | 19 | # the Fourier coefficients to spherical harmonic coefficients using
|
18 | 20 | # `plan_sph2fourier`.
|
19 | 21 | #
|
20 | 22 | # In the basis of spherical harmonics, it is plain to see the
|
21 |
| -# addition theorem in action, since P₄(x⋅y) should only consist of |
22 |
| -# exact-degree-4 harmonics. |
| 23 | +# addition theorem in action, since $P_4(x\cdot y)$ should only consist of |
| 24 | +# exact-degree-$4$ harmonics. |
23 | 25 | #
|
24 |
| -# For the storage pattern of the arrays, please consult the documentation. |
25 |
| -############# |
| 26 | +# For the storage pattern of the arrays, please consult the |
| 27 | +# [documentation](https://MikaelSlevinsky.github.io/FastTransforms). |
26 | 28 |
|
27 | 29 | function threshold!(A::AbstractArray, ϵ)
|
28 | 30 | for i in eachindex(A)
|
|
33 | 35 |
|
34 | 36 | using FastTransforms, LinearAlgebra
|
35 | 37 |
|
36 |
| -# The colatitudinal grid (mod π): |
| 38 | +# The colatitudinal grid (mod $\pi$): |
37 | 39 | N = 5
|
38 | 40 | θ = (0.5:N-0.5)/N
|
39 | 41 |
|
40 |
| -# The longitudinal grid (mod π): |
| 42 | +# The longitudinal grid (mod $\pi$): |
41 | 43 | M = 2*N-1
|
42 | 44 | φ = (0:M-1)*2/M
|
43 | 45 |
|
44 |
| -# Arbitrarily, we place x at the North pole: |
| 46 | +# Arbitrarily, we place $x$ at the North pole: |
45 | 47 | x = [0,0,1]
|
46 | 48 |
|
47 | 49 | # Another vector is completely free:
|
48 | 50 | y = normalize([.123,.456,.789])
|
49 | 51 |
|
50 |
| -# Thus z ∈ 𝕊² is our variable vector, parameterized in spherical coordinates: |
| 52 | +# Thus $z \in \mathbb{S}^2$ is our variable vector, parameterized in spherical coordinates: |
51 | 53 | z = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]
|
52 | 54 |
|
53 |
| -# The degree-4 Legendre polynomial is: |
| 55 | +# The degree-$4$ Legendre polynomial is: |
54 | 56 | P4 = x -> (35*x^4-30*x^2+3)/8
|
55 | 57 |
|
56 | 58 | # On the tensor product grid, our function samples are:
|
57 | 59 | F = [(P4(z(θ,φ)⋅y) - P4(x⋅y))/(z(θ,φ)⋅y - x⋅y) for θ in θ, φ in φ]
|
58 | 60 |
|
| 61 | +# We precompute a spherical harmonic--Fourier plan: |
59 | 62 | P = plan_sph2fourier(F)
|
| 63 | + |
| 64 | +# And an FFTW Fourier analysis plan on $\mathbb{S}^2$: |
60 | 65 | PA = plan_sph_analysis(F)
|
61 | 66 |
|
62 |
| -# Its spherical harmonic coefficients demonstrate that it is degree-3: |
| 67 | +# Its spherical harmonic coefficients demonstrate that it is degree-$3$: |
63 | 68 | V = PA*F
|
64 | 69 | U3 = threshold!(P\V, 400*eps())
|
65 | 70 |
|
66 |
| -# Similarly, on the tensor product grid, the Legendre polynomial P₄(z⋅y) is: |
| 71 | +# Similarly, on the tensor product grid, the Legendre polynomial $P_4(z\cdot y)$ is: |
67 | 72 | F = [P4(z(θ,φ)⋅y) for θ in θ, φ in φ]
|
68 | 73 |
|
69 |
| -# Its spherical harmonic coefficients demonstrate that it is exact-degree-4: |
| 74 | +# Its spherical harmonic coefficients demonstrate that it is exact-degree-$4$: |
70 | 75 | V = PA*F
|
71 | 76 | U4 = threshold!(P\V, 3*eps())
|
72 | 77 |
|
73 |
| -nrm1 = norm(U4); |
| 78 | +# The $L^2(\mathbb{S}^2)$ norm of the function is: |
| 79 | +nrm1 = norm(U4) |
74 | 80 |
|
75 |
| -# Finally, the Legendre polynomial P₄(z⋅x) is aligned with the grid: |
| 81 | +# Finally, the Legendre polynomial $P_4(z\cdot x)$ is aligned with the grid: |
76 | 82 | F = [P4(z(θ,φ)⋅x) for θ in θ, φ in φ]
|
77 | 83 |
|
78 | 84 | # It only has one nonnegligible spherical harmonic coefficient.
|
79 | 85 | # Can you spot it?
|
80 | 86 | V = PA*F
|
81 | 87 | U4 = threshold!(P\V, 3*eps())
|
82 | 88 |
|
83 |
| -# That nonnegligible coefficient should be approximately √(2π/(4+1/2)), |
| 89 | +# That nonnegligible coefficient should be approximately `√(2π/(4+1/2))`, |
84 | 90 | # since the convention in this library is to orthonormalize.
|
85 | 91 |
|
86 |
| -nrm2 = norm(U4); |
| 92 | +nrm2 = norm(U4) |
87 | 93 |
|
88 |
| -# Note that the integrals of both functions P₄(z⋅y) and P₄(z⋅x) and their |
89 |
| -# L²(𝕊²) norms are the same because of rotational invariance. The integral of |
| 94 | +# Note that the integrals of both functions $P_4(z\cdot y)$ and $P_4(z\cdot x)$ and their |
| 95 | +# $L^2(\mathbb{S}^2)$ norms are the same because of rotational invariance. The integral of |
90 | 96 | # either is perhaps not interesting as it is mathematically zero, but the norms
|
91 | 97 | # of either should be approximately the same.
|
92 | 98 |
|
93 |
| -@show nrm1 |
94 |
| -@show nrm2 |
95 |
| -@show nrm1 ≈ nrm2 |
| 99 | +nrm1 ≈ nrm2 |
0 commit comments