Utilities
While not usually the star of the show, the following utilities can be quite helpful for actually using the rest of the code.
Complex powers
One common task we find when working with spherical functions is the computation of a range of integer powers of some complex number β so much so that it can be best to pre-compute the powers and cache their values. While a naive approach is generally quite accurate, and reasonably fast, we can do a little better with a specialized routine.
SphericalFunctions.complex_powers! β Methodcomplex_powers!(zpowers, z)Compute integer powers of z from z^0 through z^m, recursively, where m is one less than the length of the input zpowers vector.
Note that z is assumed to be normalized, with complex amplitude approximately 1.
See also: complex_powers
SphericalFunctions.complex_powers β Methodcomplex_powers(z, m)Compute integer powers of z from z^0 through z^m, recursively.
Note that z is assumed to be normalized, with complex amplitude approximately 1.
This algorithm is mostly due to Stoer and Bulirsch in "Introduction to Numerical Analysis" (page 24) β with a little help from de Moivre's formula, which is essentially exp(iΞΈ)βΏ = exp(inΞΈ), as well as my own alterations to deal with different behaviors in different quadrants.
There isn't usually a huge advantage to using this specialized function. If you just need a particular power, it will generally be far more efficient and just as accurate to compute either exp(iΞΈ)βΏ or exp(inΞΈ) explicitly. However, if you need all powers from 0 to m, this function is about 10 or 5 times faster than those options, respectively, for large m. Like those options, this function is numerically stable, in the sense that its errors are usually smaller than m times the error from machine-precision errors in the input argument β or at worst about 50% larger, which occurs as the phase approaches multiples of Ο/2.
See also: complex_powers!
Sizes of and indexing into $π$, $d$, and $Y$ data
By $Y$ data, we mean anything indexed like $Y_{\ell, m}$ modes; by $D$ data, we mean anything indexed like Wigner's $\mathfrak{D}$ matrices, or special subsets of them, like the $H$ matrices.
SphericalFunctions.WignerDindex β FunctionWignerDindex(β, mβ², m, mβ²βββ=β)Compute index into Wigner π matrix
See also WignerDrange and WignerDsize.
Notes
This assumes that the Wigner π matrix is arranged as
[
π(β, mβ², m)
for β in range(ββα΅’β, ββββ+1)
for mβ² in range(-min(β, mβ²βββ), min(β, mβ²βββ)+1)
for m in range(-β, β+1)
]SphericalFunctions.WignerDrange β FunctionWignerDrange(ββββ, mβ²βββ=ββββ)Create an array of (β, m', m) indices as in π array
See also WignerDsize and WignerDindex.
Notes
This assumes that the Wigner π matrix is arranged as
[
π(β, mβ², m)
for β in range(ββα΅’β, ββββ+1)
for mβ² in range(-min(β, mβ²βββ), min(β, mβ²βββ)+1)
for m in range(-β, β+1)
]SphericalFunctions.WignerDsize β MethodWignerDsize(ββββ, mβ²βββ=ββββ)Compute total size of Wigner π matrix
See also WignerDrange and WignerDindex.
Notes
This assumes that the Wigner π matrix is arranged as
[
π(β, mβ², m)
for β in ββα΅’β:ββββ
for mβ² in -min(β, mβ²βββ):min(β, mβ²βββ)
for m in -β:β
]SphericalFunctions.WignerHindex β FunctionWignerHindex(β, mβ², m, mβ²βββ)Index to "wedge" arrays.
See also WignerHsize and WignerHrange.
Notes
Here, it is assumed that only data with mβ₯|m'| are stored, and only corresponding values are passed. We also assume |m|β€β and |m'|β€β. Neither of these are checked. The wedge array that this function indexes is ordered as
[
H(β, mβ², m) for β in range(ββββ+1)
for mβ² in range(-min(β, mβ²βββ), min(β, mβ²βββ)+1)
for m in range(abs(mβ²), β+1)
]SphericalFunctions.WignerHrange β FunctionWignerHrange(ββββ, mβ²βββ=ββββ)Create an array of (β, m', m) indices as in H array
See also WignerHsize and WignerHindex
Notes
Here, it is assumed that only data with mβ₯|m'| are stored, and only corresponding values are passed. We also assume |m|β€β and |m'|β€β. Neither of these are checked. The wedge array that this function indexes is ordered as
[
H(β, mβ², m) for β in range(ββββ+1)
for mβ² in range(-min(β, mβ²βββ), min(β, mβ²βββ)+1)
for m in range(abs(mβ²), β+1)
]SphericalFunctions.WignerHsize β MethodWignerHsize(ββββ, mβ²βββ=ββββ)Total size of array of wedges of width mβ²βββ up to ββββ. If mβ²βββ is not given, it defaults to ββββ.
See also WignerHrange and WignerHindex.
Notes
Here, it is assumed that only data with mβ₯|mβ²| are stored, and only corresponding values are passed. We also assume |m|β€β and |mβ²|β€β. Neither of these are checked. The wedge array that this function indexes is ordered as
[
H(β, mβ², m) for β in 0:ββββ
for mβ² in -min(β, mβ²βββ):min(β, mβ²βββ)
for m in abs(mβ²):β
]SphericalFunctions.Yindex β FunctionYindex(β, m, ββα΅’β=0)Compute index into array of mode weights
Parameters
β : int Integer satisfying ββα΅’β <= β <= ββββ m : int Integer satisfying -β <= m <= β ββα΅’β : int, optional Integer satisfying 0 <= ββα΅’β. Defaults to 0.
Returns
i : int Index of a particular element of the mode-weight array as described below
See Also
Ysize : Total size of array of mode weights Yrange : Array of (β, m) indices corresponding to this array
Notes
This assumes that the modes are arranged (with fixed s value) as
[
Y(s, β, m)
for β in range(ββα΅’β, ββββ+1)
for m in range(-β, β+1)
]SphericalFunctions.Yrange β MethodYrange(ββα΅’β, ββββ)Create an array of (β, m) indices as in Y array
Parameters
ββα΅’β : int Integer satisfying 0 <= ββα΅’β <= ββββ ββββ : int Integer satisfying 0 <= ββα΅’β <= ββββ
Returns
i : int Total size of array of mode weights arranged as described below
See Also
Ysize : Total size of array of mode weights Yindex : Index of a particular element of the mode weight array
Notes
This assumes that the modes are arranged (with fixed s value) as
[
Y(s, β, m)
for β in range(ββα΅’β, ββββ+1)
for m in range(-β, β+1)
]SphericalFunctions.Ysize β MethodYsize(ββββ)Compute total size of array of mode weights
See Also
Yrange : Array of (β, m) indices corresponding to this array Yindex : Index of a particular element of the mode weight array
Notes
This assumes that the modes are arranged (with fixed s value) as
[
Y(s, β, m)
for β in range(ββα΅’β, ββββ+1)
for m in range(-β, β+1)
]SphericalFunctions._WignerHindex β MethodWignerHindex(β, mβ², m, mβ²βββ)
Helper function for WignerHindex, with more constraints.
This function assumes that m β₯ |mβ²|. The main WignerHindex function uses symmetries of the H array to account for cases that violate this assumption. (But note that both that function and this one assume that |m| β€ β and |mβ²| β€ β.)
SphericalFunctions.deduce_limits β Functiondeduce_limits(ysize, [βmin])Deduce the value of (βmin, βmax) that produces $Y$ arrays of the given size.
If βmin is not given, it is assumed to be 0. If it is set to nothing, the smallest possible value of βmin will be used. However, note that this is not a well-posed problem; multiple combinations of (βmin, βmax) can give rise to $Y$ arrays of the same size.
See also Ysize
SphericalFunctions.phi_theta β Methodphi_theta(NΟ, NΞΈ, [T=Float64])Construct (phi, theta) grid in order expected by this package.
See also theta_phi for the opposite ordering.
SphericalFunctions.theta_phi β Methodtheta_phi(NΞΈ, NΟ, [T=Float64])Construct (theta, phi) grid in spinsfast order.
Note that this order is different from the one assumed by this package; use phi_theta for the opposite ordering.
Combinatorics
Spherical functions frequently involve binomial coefficients and similar terms, with arguments proportional to $β$, which we aim to allow to be very large β of order 1,000 or more. Unfortunately, due to combinatorical explosions, this is frequently infeasible with naive methods. Here, we collect any specialized methods that help us beat the limits.
SphericalFunctions.sqrtbinomial β Methodsqrtbinomial(n, k, [T])Evaluate the square-root of the binomial coefficient binomial(n,k) for large coefficients.
Ordinarily, when n and k are standard Int arguments, the built-in binomial function will overflow around n=66, because it results in Ints. We need much larger values. This function, which is based on a related one in SpecialFunctions.jl, returns reasonably accurate results up to n β 1026 when k β n/2 (which is the case of interest in many applications in this package).
Computations are carried out (and returned) in type T, which defaults to Float64.