Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : /*!
7 : * \brief Holds functions that converts spatial tensors between scalar-Ylm
8 : * and tensor-Ylm basis, and functions that filter by doing such a conversion.
9 : *
10 : * \details We expand a spatial tensor of arbitrary rank
11 : * in terms of Cartesian components as follows:
12 : * \begin{align}
13 : * {\mathbf T} &= \sum_{\ell,m,\tilde{A}} {}_0 Y_{\ell m}
14 : * T^{\tilde A}_{\ell m}\mathbf{e}_{\tilde A},
15 : * \end{align}
16 : * where ${}_0 Y_{\ell m}$ are the usual scalar spherical harmonics,
17 : * $T^{\tilde A}_{\ell m}$ are complex-valued expansion coefficients,
18 : * $\tilde A$ is a multi-index
19 : * $\tilde{A}=(\tilde{a}_1,\tilde{a}_2,\tilde{a}_3,\ldots)$, where each
20 : * of the $\tilde{a}_i$ take on values from 0 to 2, referring to
21 : * either $\mathbf{e}_x$, $\mathbf{e}_y$, or $\mathbf{e}_z$.
22 : * The sum over $\tilde A$ goes over all $\tilde A$ of the same rank.
23 : *
24 : * Similarly, we can expand the same spatial
25 : * tensor in terms of a complex tetrad:
26 : * \begin{align}
27 : * {\mathbf T} &= \sum_{\ell,m,A} {}_{s(\!A\!)}Y_{\ell m}
28 : * T^A_{\ell m} \mathbf{e}_A,
29 : * \end{align}
30 : * where ${}_sY_{\ell m}$ are the spin-weighted harmonics,
31 : * $A$ is a multi-index $A=(a_1,a_2,a_3,\ldots)$, where each of the
32 : * $a_i$ refer to either $\mathbf{l}$,$\mathbf{m}$, or $\mathbf{\bar{m}}$,
33 : * and $s(\!A\!)$ is the spin weight, which is a function of
34 : * $A$: each $\mathbf{l}$ in $A$ adds the value zero, each $\mathbf{m}$
35 : * adds the value -1, and each $\mathbf{\bar{m}}$ adds the value +1.
36 : * (These values are the spin weights of the complex conjugates
37 : * of the basis vectors.)
38 : *
39 : * To allow for a compact notation where we don't need to write different
40 : * equations for different combinations of basis vectors,
41 : * we define two auxiliary quantities $k_j$ and $m_j$
42 : * that are functions of the Cartesian basis vectors:
43 : * \begin{align}
44 : * k_j(\mathbf{e}_x) &= -j,\\
45 : * k_j(\mathbf{e}_y) &= i,\\
46 : * k_j(\mathbf{e}_z) &= 1/\sqrt{2},\\
47 : * m_j(\mathbf{e}_x) &= j,\\
48 : * m_j(\mathbf{e}_y) &= j,\\
49 : * m_j(\mathbf{e}_z) &= 0,
50 : * \end{align}
51 : * where the $i$ above is the imaginary unit $i=\sqrt{-1}$.
52 : * These quantities will be used in explicit formulas for transformations
53 : * between expansion coefficients.
54 : *
55 : * ### Spherepack coefficients
56 : *
57 : * Instead of using the actual spherical-harmonic or
58 : * spin-weighted spherical-harmonic coefficients $T^{B}_{\ell' m'}$ and
59 : * $T^{\tilde A}_{\ell'm'}$, we often work with
60 : * coefficients computed by Spherepack, which have a different normalization
61 : * and phase than the standard coefficients. While Spherepack internally
62 : * stores real matrices it calls $a_{lm}$ and $b_{lm}$ (see Spherepack
63 : * documentation for details), we simplify things by defining
64 : * complex Spherepack coefficients
65 : * ${\breve T}^{B}_{\ell' m'}$ and ${\breve T}^{\tilde A}_{\ell'm'}$ that
66 : * have Spherepack normalization. These are related to the standard coefficients
67 : * by
68 : * \begin{align}
69 : * {\breve T}^{A}_{\ell m} &= (-1)^m \sqrt{\frac{2}{\pi}} T^{B}_{\ell m},\\
70 : * {\breve T}^{\tilde A}_{\ell m} &= (-1)^m \sqrt{\frac{2}{\pi}}
71 : * T^{\tilde A}_{\ell m}.
72 : * \end{align}
73 : * Some of the formulas below are slightly different when acting on
74 : * Spherepack coefficients compared to standard coefficients.
75 : *
76 : * ## Transformations
77 : *
78 : * We can define the following transformations between the expansion
79 : * coefficients $T^A_{\ell m}$ and $T^{\tilde A}_{\ell m}$:
80 : * \begin{align}
81 : * T^{\tilde B}_{\ell' m'} &=
82 : * \sum_{\ell,m,A} C_{\ell' m' A}^{\ell m\tilde{B}} T^A_{\ell m},\\
83 : * T^{B}_{\ell' m'} &= \sum_{\ell m \tilde{A}}
84 : * C_{\ell' m'\tilde{A}}^{\ell mB}
85 : * T^{\tilde A}_{\ell m},
86 : * \end{align}
87 : * where the above transformations define the coefficients
88 : * $C_{\ell' m' A}^{\ell m\tilde{B}}$ and $C_{\ell' m'\tilde{A}}^{\ell mB}$.
89 : * Analytic expressions for these coefficients are derived in Klinger
90 : * and Scheel (in prep) and will be used here to compute the coefficients.
91 : *
92 : * For real-valued tensors, the coefficients $T^A_{\ell m}$ and
93 : * $T^{\tilde A}_{\ell m}$ for negative $m$ can be computed from the
94 : * complex conjugates of the same coefficients with positive $m$. Therefore
95 : * we store and loop over only nonnegative $m$. This means we can write
96 : * \begin{align}
97 : * T^{\tilde B}_{\ell m} &= \sum_{A, \ell',m'\geq 0}
98 : * \left[
99 : * C_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}
100 : * +{\hat C}_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}{}^\star
101 : * \right]\label{eq:S2C},\\
102 : * T^{B}_{\ell m} &= \sum_{\tilde{A},\ell',m'\geq 0}
103 : * \left[
104 : * C_{\ell m\tilde{A}}^{\ell' m' B}
105 : * T^{\tilde A}_{\ell' m'}
106 : * +{\hat C}_{\ell m\tilde{A}}^{\ell' m' B}
107 : * T^{\tilde A}_{\ell' m'}{}^\star
108 : * \right]\label{eq:C2S},
109 : * \end{align}
110 : * where
111 : * \begin{align}
112 : * {\hat C}_{\ell m A}^{\ell' m'\tilde{B}} &=
113 : * (1-\delta_{0m'})C_{\ell m A^\star}^{\ell' -m'\tilde{B}}
114 : * (-1)^{S(A)+m'},\\
115 : * {\hat C}_{\ell m\tilde{A}}^{\ell' m' B} &=
116 : * (1-\delta_{0m'})C_{\ell m \tilde{A}}^{\ell' -m'B}
117 : * (-1)^{m'}.
118 : * \end{align}
119 : *
120 : * The functions FillCartToSphere and FillSphereToCart fill
121 : * sparse matrices that encode Eqs. $(\ref{eq:C2S})$ and
122 : * $(\ref{eq:S2C})$. When working on Spherepack coefficients, the functions
123 : * FillCartToSphere and FillSphereToCart instead
124 : * fill sparse matrices that encode
125 : * \begin{align}
126 : * {\breve T}^{\tilde B}_{\ell m} &= (-1)^m\sum_{A, \ell',m'\geq 0}(-1)^{m'}
127 : * \left[
128 : * C_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}
129 : * +{\hat C}_{\ell m A}^{\ell' m'\tilde{B}}
130 : * {\breve T}^A_{\ell' m'}{}^\star
131 : * \right]\label{eq:S2CSpherepack},\\
132 : * {\breve T}^{B}_{\ell m} &= (-1)^m\sum_{\tilde{A},\ell',m'\geq 0}(-1)^{m'}
133 : * \left[
134 : * C_{\ell m\tilde{A}}^{\ell' m' B}
135 : * T^{\tilde A}_{\ell' m'}
136 : * +{\hat C}_{\ell m\tilde{A}}^{\ell' m' B}
137 : * {\breve T}^{\tilde A}_{\ell' m'}{}^\star
138 : * \right]\label{eq:C2SSpherepack}.
139 : * \end{align}
140 : *
141 : * ## Filtering
142 : *
143 : * Consider starting with coefficients $T^{\tilde A}_{\ell' m'}$,
144 : * transforming to spin-weighted harmonics, applying a filter to
145 : * the spin-weighted harmonic coefficients, and transforming back.
146 : * We can write this entire operation as
147 : * \begin{align}
148 : * F_{\ell m \tilde{D}}^{\ell'' m''\tilde{A}} &=
149 : * \sum_{\ell'=0}^{\ell_{\mathrm{cut}}^--1} C^{\ell' m' \tilde{A}}_{\ell m B}
150 : * C^{\ell'' m'' B}_{\ell' m' \tilde{D}}
151 : * + \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{cut}}^+}
152 : * C^{\ell' m' \tilde{A}}_{\ell m B}
153 : * C^{\ell'' m'' B}_{\ell' m' \tilde{D}} f(\ell'),
154 : * \end{align}
155 : * where $f(\ell')$ is a filter function,
156 : * $\ell_{\mathrm{cut}}^{-}$ is the smallest $\ell$ mode that is filtered,
157 : * and
158 : * $\ell_{\mathrm{cut}}^{+}$ is the largest $\ell$ mode that is retained.
159 : *
160 : * There are two cases we care about.
161 : * The first is simple "Heaviside filtering", where
162 : * $\ell_{\mathrm{cut}}^{-} = \ell_{\mathrm{cut}}^{+}$ and $f(\ell')=1$.
163 : *
164 : * The second case is
165 : * \begin{align}
166 : * f(\ell') &= \exp\left[-\alpha
167 : * \left(\frac{\ell'}{\ell_{\mathrm{cut}}^++1}\right)^{2 \sigma}\right],\\
168 : * \ell_{\mathrm{cut}}^- &= \mathrm{ceil}\left[(\ell_{\mathrm{cut}}^++1)
169 : * \left(\frac{\epsilon}{\alpha}\right)^{1/(2\sigma)}\right],
170 : * \end{align}
171 : * where $\alpha$ is a parameter we choose to be 36, $\sigma$ is an integer
172 : * parameter typically between 28 and 32, and
173 : * $\epsilon$ is machine roundoff.
174 : * Note that the filter remains smooth at $\ell' = \ell_{\mathrm{cut}}^-$ and
175 : * it reduces to the simple filter as $\sigma \to \infty$.
176 : *
177 : * As with the transformations, we sum over only nonnegative $m$, so we can
178 : * write the action of the filter as
179 : * \begin{align}
180 : * T^{\tilde B}_{\ell m} &= \sum_{{\tilde A}, \ell',m'\geq 0}
181 : * \left[
182 : * F_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
183 : * T^{\tilde A}_{\ell' m'}
184 : * +{\hat F}_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
185 : * T^{\tilde A}_{\ell' m'}{}^\star
186 : * \right]\label{eq:Filter},
187 : * \end{align}
188 : * where
189 : * \begin{align}
190 : * {\hat F}_{\ell m\tilde{A}}^{\ell' m' \tilde{B}} &=
191 : * (1-\delta_{0m'})F_{\ell m \tilde{A}}^{\ell' -m' \tilde{B}}
192 : * (-1)^{m'}.
193 : * \end{align}
194 : *
195 : * When filtering Spherepack coefficients, the expression is
196 : * \begin{align}
197 : * {\breve T}^{\tilde B}_{\ell m} &=
198 : * (-1)^m\sum_{{\tilde A}, \ell',m'\geq 0}(-1)^{m'}
199 : * \left[
200 : * F_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
201 : * {\breve T}^{\tilde A}_{\ell' m'}
202 : * +{\hat F}_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
203 : * {\breve T}^{\tilde A}_{\ell' m'}{}^\star
204 : * \right]\label{eq:FilterSpherepack}.
205 : * \end{align}
206 : *
207 : * The function FillFilter encapsulates Eq. $(\ref{eq:Filter})$ or
208 : * Eq. $(\ref{eq:FilterSpherepack})$,
209 : * depending what type of coefficients are being filtered.
210 : */
211 : namespace ylm::TensorYlm {
212 :
213 : /// Instances of this class are passed to FillFilter,
214 : /// FillCartToSphere, and FillSphereToCart to identify how the
215 : /// coefficients are normalized. This makes a difference because the
216 : /// normalization is not just a constant factor.
217 1 : enum class CoefficientNormalization {
218 : /// The standard normalization of the spherical-harmonic coefficients.
219 : Standard,
220 : /// The Spherepack normalization. Spherepack coefficients are standard
221 : /// coefficients multiplied by $(-1)^m \sqrt{2/\pi}$.
222 : Spherepack
223 : };
224 : } // namespace ylm::TensorYlm
|
