:mod:`commander.sphere`: Working with spherical data
=====================================================

The :mod:`commander.sphere` package contains utilities for working
with spherical data. The tools here are independent of the rest of
Commander; the rest of Commander is heavily based on the data format
and conventions present here.

Spherical harmonic conventions
------------------------------

Ordering scheme
'''''''''''''''

Across Commander we use :math:`m`-major ordering with interleaved
positive and negative $m$; i.e., the coefficients $(\ell, m)$ are
ordered first by $|m|$, then by $\ell$, and finally by ordering
the positive $m$ before the negative $m$.

**Examples:** For all $m$ and $0 \le \ell \le 2$, the coefficient order is

.. math::
    
    \left[a_{0,0}, a_{1,0}, a_{2,0}, a_{1,1}, a_{1,-1}, a_{2,1}, a_{2,-1}, a_{2,2}, a_{2,-2} \right ],

whereas for storing only $3 \le \ell \le 4$ one gets

.. math::
    
    \left[a_{3,0}, a_{4,0}, a_{3,1}, a_{3,-1}, a_{4,1}, a_{4,-1}, a_{3,2}, a_{3,-2},
    \dots a_{4,4}, a_{4,-4} \right ].

The routines :func:`commander.sphere.mmajor.lm_to_idx` and
:func:`commander.sphere.mmajor.idx_to_lm` can translate between
array indices and $(\ell, m)$.

Spherical harmonic transforms process the coefficients one $m$ at the
time, making :math:`m`-major ordering the most efficient one for
SHTs. Memory locality aside, it is crucial that all coefficients for a
given $m$ reside on the same node.
   
In various computations $a_{\ell,m}$ and $a_{\ell,-m}$ are often
related and treated together, so it makes sense to keep them
together. When interleaving in this particular way with positive $m$
first, converting between *half-complex* and *real* spherical harmonic
formats is just a matter of in-place coefficient scaling with no
reordering (except possibly for $m=0$). In the half-complex format,
$a_{\ell,m}$ is seens as the real component and $a_{\ell,-m}$ the
imaginary.

We avoid any padding. Padding could have made for more convenient
indexing, but is impractical in linear algebra, as one is forced to
also pad matrices accordingly. Also it becomes easy to mis-count the
number of degrees of freedom for vectors and so on.

Real vs. complex spherical harmonics
''''''''''''''''''''''''''''''''''''

We deal with real fields on the sphere, and with the
traditional complex definition of spherical harmonics we have
$a_{\ell m}=(-1)^m a_{\ell -m}$ and the negative $m$ are redundant. *However*,
they must still be included in many computations (you can **not**
stack only the non-negative $m$ in a vector and do linear
algebra with the result).

Therefore, to make sure that linear algebra with spherical harmonic
vectors always work out well (and to avoid a very significant speed
penalty), we mainly use *real spherical
harmonics*, related to complex spherical harmonics by the equations

.. math::
  :nowrap:
  
  \begin{eqnarray}
     a^R_{\ell,0} &=& a^C_{\ell,0} & \\
     a^R_{\ell,m} &=& \sqrt{2} \text{Re}(a^C_{\ell,m}) & \text{for $m>0$}\\
     a^R_{\ell,m} &=& \sqrt{2} \text{Im}(a^C_{\ell,-m}) & \text{for $m<0$}\\
   \end{eqnarray}

Algebraically, the transformation from complex to real spherical harmonics
is an orthonormal transformation.

.. warning::

  While the conversion of *vectors* is trivial, the conversion of
  *matrices* between the real and complex spherical harmonic bases is
  a bit less trivial, since one needs to apply the basis change
  operator on both sides of the matrix.  In the case of sparse
  matrices it may even sometimes be better to stick with the full expanded
  complex vectors, as the number of non-zero elements are smaller.

When it comes to storage, we do not define any *half-complex* format
(i.e., where one stores only non-negative $m$ ), as one can simply
read the real coefficients and correct for the $\sqrt{2}$ on the
fly. Both *real* and *complex* storage formats are in use, both using
the :math:`m`-major ordering described above with either real or
complex coefficients.

Reference
---------

.. automodule:: commander.sphere