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 : * ## Transformations
56 : *
57 : * We can define the following transformations between the expansion
58 : * coefficients $T^A_{\ell m}$ and $T^{\tilde A}_{\ell m}$:
59 : * \begin{align}
60 : * T^{\tilde B}_{\ell' m'} &=
61 : * \sum_{\ell,m,A} C_{\ell' m' A}^{\ell m\tilde{B}} T^A_{\ell m},\\
62 : * T^{B}_{\ell' m'} &= \sum_{\ell m \tilde{A}}
63 : * C_{\ell' m'\tilde{A}}^{\ell mB}
64 : * T^{\tilde A}_{\ell m},
65 : * \end{align}
66 : * where the above transformations define the coefficients
67 : * $C_{\ell' m' A}^{\ell m\tilde{B}}$ and $C_{\ell' m'\tilde{A}}^{\ell mB}$.
68 : * Analytic expressions for these coefficients are derived in Klinger
69 : * and Scheel (in prep) and will be used here to compute the coefficients.
70 : *
71 : * For real-valued tensors, the coefficients $T^A_{\ell m}$ and
72 : * $T^{\tilde A}_{\ell m}$ for negative $m$ can be computed from the
73 : * complex conjugates of the same coefficients with positive $m$. Therefore
74 : * we store and loop over only nonnegative $m$. This means we can write
75 : * \begin{align}
76 : * T^{\tilde B}_{\ell m} &= \sum_{A, \ell',m'\geq 0}
77 : * \left[
78 : * C_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}
79 : * +{\hat C}_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}{}^\star
80 : * \right]\label{eq:S2C},\\
81 : * T^{B}_{\ell m} &= \sum_{\tilde{A},\ell',m'\geq 0}
82 : * \left[
83 : * C_{\ell m\tilde{A}}^{\ell' m' B}
84 : * T^{\tilde A}_{\ell' m'}
85 : * +{\hat C}_{\ell m\tilde{A}}^{\ell' m' B}
86 : * T^{\tilde A}_{\ell' m'}{}^\star
87 : * \right]\label{eq:C2S},
88 : * \end{align}
89 : * where
90 : * \begin{align}
91 : * {\hat C}_{\ell m A}^{\ell' m'\tilde{B}} &=
92 : * (1-\delta_{0m'})C_{\ell m A^\star}^{\ell' -m'\tilde{B}}
93 : * (-1)^{S(A)+m'},\\
94 : * {\hat C}_{\ell m\tilde{A}}^{\ell' m' B} &=
95 : * (1-\delta_{0m'})C_{\ell m \tilde{A}}^{\ell' -m'B}
96 : * (-1)^{m'}.
97 : * \end{align}
98 : *
99 : * The functions FillCartToSphere and FillSphereToCart fill
100 : * sparse matrices that encode Eqs. $(\ref{eq:C2S})$ and
101 : * $(\ref{eq:S2C})$.
102 : *
103 : * ## Filtering
104 : *
105 : * Consider starting with coefficients $T^{\tilde A}_{\ell' m'}$,
106 : * transforming to spin-weighted harmonics, applying a filter to
107 : * the spin-weighted harmonic coefficients, and transforming back.
108 : * We can write this entire operation as
109 : * \begin{align}
110 : * F_{\ell m \tilde{D}}^{\ell'' m''\tilde{A}} &=
111 : * \sum_{\ell'=0}^{\ell_{\mathrm{cut}}^--1} C^{\ell' m' \tilde{A}}_{\ell m B}
112 : * C^{\ell'' m'' B}_{\ell' m' \tilde{D}}
113 : * + \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{cut}}^+}
114 : * C^{\ell' m' \tilde{A}}_{\ell m B}
115 : * C^{\ell'' m'' B}_{\ell' m' \tilde{D}} f(\ell'),
116 : * \end{align}
117 : * where $f(\ell')$ is a filter function,
118 : * $\ell_{\mathrm{cut}}^{-}$ is the smallest $\ell$ mode that is filtered,
119 : * and
120 : * $\ell_{\mathrm{cut}}^{+}$ is the largest $\ell$ mode that is retained.
121 : *
122 : * There are two cases we care about.
123 : * The first is simple "Heaviside filtering", where
124 : * $\ell_{\mathrm{cut}}^{-} = \ell_{\mathrm{cut}}^{+}$ and $f(\ell')=1$.
125 : *
126 : * The second case is
127 : * \begin{align}
128 : * f(\ell') &= \exp\left[-\alpha
129 : * \left(\frac{\ell'}{\ell_{\mathrm{cut}}^++1}\right)^{2 \sigma}\right],\\
130 : * \ell_{\mathrm{cut}}^- &= \mathrm{ceil}\left[(\ell_{\mathrm{cut}}^++1)
131 : * \left(\frac{\epsilon}{\alpha}\right)^{1/(2\sigma)}\right],
132 : * \end{align}
133 : * where $\alpha$ is a parameter we choose to be 36, $\sigma$ is an integer
134 : * parameter typically between 28 and 32, and
135 : * $\epsilon$ is machine roundoff.
136 : * Note that the filter remains smooth at $\ell' = \ell_{\mathrm{cut}}^-$ and
137 : * it reduces to the simple filter as $\sigma \to \infty$.
138 : *
139 : * As with the transformations, we sum over only nonnegative $m$, so we can
140 : * write the action of the filter as
141 : * \begin{align}
142 : * T^{\tilde B}_{\ell m} &= \sum_{{\tilde A}, \ell',m'\geq 0}
143 : * \left[
144 : * F_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
145 : * T^{\tilde A}_{\ell' m'}
146 : * +{\hat F}_{\ell m {\tilde A}}^{\ell' m'\tilde{B}}
147 : * T^{\tilde A}_{\ell' m'}{}^\star
148 : * \right]\label{eq:Filter},
149 : * \end{align}
150 : * where
151 : * \begin{align}
152 : * {\hat F}_{\ell m\tilde{A}}^{\ell' m' \tilde{B}} &=
153 : * (1-\delta_{0m'})F_{\ell m \tilde{A}}^{\ell' -m' \tilde{B}}
154 : * (-1)^{m'}.
155 : * \end{align}
156 : *
157 : * The function FillFilter encapsulates Eq. $(\ref{eq:Filter})$.
158 : */
159 1 : namespace ylm::TensorYlm {} // namespace ylm::TensorYlm
|
