Working with spherical harmonics

20/09/22

This notebook gives a brief introduction to some of the conventions and functionality for handling spherical harmonics and related functions, using low-level ePSproc routines. Note some higher-level routines are available in the base class.

Expansions in (complex) spherical harmonics

In general, expansions in complex harmonics are used, with expansion parameters \(\beta_{L,M}\). A given function is then written as:

\begin{equation} I(\theta,\phi)=\sum_{l,m}\beta_{l,m}Y_{l,m}(\theta,\phi)\label{eq:sph-harmonics} \end{equation}

For additional control and conversion, the SHtools library can also be used - this includes some useful utility functions including converting between different forms (e.g. real and complex forms), and basic plotters.

For use of real spherical harmonics, see the ‘working with real harmonics notebook’.

Expansions in Legendre polynomials

For cylindrically-symmetric cases (\(m=0\)), general expansions in Legendre polynomials are also sufficient and often used:

\begin{equation} I(\theta)=\sum_{l}\beta_{l}P^0_l(\cos(\theta))\label{eq:lg-expansion} \end{equation}

Where the “Legendre” and “spherical harmonic” expansion parameters for the harmonics defined above can be related by:

\begin{equation} \beta^{Sph}_{l,0} = \sqrt{(2L+1)/4\pi}\beta^{Lg}_{l} \end{equation}

Expansions in symmetrized harmonics

In some case, symmetrized or generalised harmonics are useful. These are defined similarly, but with expansion coeffs additionally set by (point group) symmetry. A basic implementation can be found in the PEMtk package, defined as:

\begin{equation} X_{hl}^{\Gamma\mu*}(\theta,\phi)=\sum_{\lambda}b_{hl\lambda}^{\Gamma\mu}Y_{l,\lambda}(\theta,\phi)\label{eq:symm-harmonics} \end{equation}

Where the \(b_{hl\lambda}^{\Gamma\mu}\) are defined by symmetry. General function expansions can then be written as a set of spherical harmonics including symmetrization, or an equivalent expansion in symmetrized harmonics (here not all indices may be necessary):

\begin{equation} I(\theta,\phi)=\sum_{h,l,\Gamma,\mu}\beta_{hl}^{\Gamma\mu}X_{hl}^{\Gamma\mu*}(\theta,\phi)=\sum_{h,l,\Gamma,\mu}\sum_{\lambda}b_{hl\lambda}^{\Gamma\mu}\beta_{l,\lambda}Y_{l,\lambda}(\theta,\phi)\label{eq:symm-harmonics-2} \end{equation}

Numerical implementation

Various tools are currently implemented in ePSProc (as of Sept. 2022), and are illustrated below. In particular:

ep.sphCalc contains the base routines for generation of harmonics.
ep.sphPlot implements basic plotting routines.
ep.sphFuncs.sphConv implements additional handling and tools, including conversion routines.

The base class implements some more sophisticated plotting options.

The default routines in ePSproc make use of scipy.special.sph_harm as the default calculation routine. Note the \((\theta, \phi)\) definition, and normalisation, corresponding to the usual “physics” convention (\(\theta\) defined from the z-axis, includes the Condon-Shortley phase, and orthonormalised):

\begin{equation} Y_{l,m}(\theta,\phi) = (-1)^m\sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!}}e^{i m \phi} P^m_l(\cos(\theta)) \end{equation}

For more details, see also Wikipaedia, which has matching defintions plus further discussion. Note that in the Scipy routines the Condon-Shortley phase term, \((-1)^m\), is actually included in the associated Legendre polynomial function \(P^m_l\), but is written explicitly above.

Conjugates and \(\pm m\) sign swaps are implemented in ep.sphFuncs.sphConv.sphConj as per Blum:

\begin{equation} Y_{l,m} = (-1)^m Y_{l,-m}^* \end{equation}

Imports

[1]:

import xarray as xr
import pandas as pd
import numpy as np
import epsproc as ep

# Set compact XR repr
xr.set_options(display_expand_data = False)

OMP: Info #273: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

* sparse not found, sparse matrix forms not available.
* natsort not found, some sorting functions not available.
* Hvplot not found, some hvPlotters may not be available. See https://hvplot.holoviz.org/user_guide/Gridded_Data.html for package details.
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.

* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available.

[1]:

<xarray.core.options.set_options at 0x7fe60c370b50>

Computing spherical harmonics on a grid

A basic wrapper for the backends is provided by ep.sphCalc. This computes spherical harmonics for all orders up to Lmax, for a given angular resolution or set of angles, and returns an Xarray. Details of the nature of the harmonics is output in the Xarray, as self.attrs['harmonics']. In most cases this will be used by other functions as required, or can be overridden at the function call.

[2]:

# Compute harmonics for 50x50 (theta,phi) grid
Isph = ep.sphCalc(Lmax = 2, res = 50)
Isph

[3]:

# Compute harmonics for 50x25 (theta,phi) grid
I = ep.sphCalc(Lmax = 2, res = [50,25])
I

[4]:

# The sum of these can be plotted directly... although may not be very useful
ep.sphSumPlotX(I)

*** WARNING: plot dataset has min value < 0, min = (-0.8400465549302628-0.43869471033269675j). This may be unphysical and/or result in plotting issues.
Sph plots:
Plotting with facetDims=Eke, pType=a with backend=mpl.

[4]:

[<Figure size 720x480 with 1 Axes>]

../_images/special_topics_ePSproc_docs_working_with_spherical_harmonics_200922_8_2.png

[5]:

# To plot components subselect on l first
# Set backend to Plotly for interactive panel plot
ep.sphSumPlotX(I.sel(l=2), facetDim = 'm', backend='pl');

*** WARNING: plot dataset has min value < 0, min = (-0.38607574059609784-9.456128398989102e-17j). This may be unphysical and/or result in plotting issues.
Sph plots:
Plotting with facetDims=m, pType=a with backend=pl.
*** Plotting for [1,1,0]
*** Plotting for [1,2,1]
*** Plotting for [1,3,2]
*** Plotting for [1,4,3]
*** Plotting for [1,5,4]

[6]:

#  Cart plots can be constructed directly with Xarray methods (NOTE this needs unstack() and .real or .imag)
I.real.unstack().plot(y='Theta', x='Phi', row='l', col='m', robust=True)

# For more data type options use ep.plotTypeSelector, e.g. plot phase
ep.plotTypeSelector(I.unstack('LM'), pType = 'phase').plot(y='Theta', x='Phi', row='l', col='m', robust=True)

[6]:

<xarray.plot.facetgrid.FacetGrid at 0x7fe4b5dfbe20>

../_images/special_topics_ePSproc_docs_working_with_spherical_harmonics_200922_10_1.png

../_images/special_topics_ePSproc_docs_working_with_spherical_harmonics_200922_10_2.png

[7]:

# Change plot style to default MPL
ep.basicPlotters.setPlotters(resetMpl=True)

* Set Holoviews with bokeh.

[8]:

# For Legendre Polynomials, set fnType = 'lg'
Ilg = ep.sphCalc(Lmax = 2, res = 50, fnType = 'lg')
ep.sphSumPlotX(Ilg)

Sph plots:
Plotting with facetDims=Eke, pType=a with backend=mpl.

[8]:

[<Figure size 640x480 with 1 Axes>]

../_images/special_topics_ePSproc_docs_working_with_spherical_harmonics_200922_12_2.png

[9]:

Ilg

Setting coefficients

To manually set coefficients from lists or arrays, use the setBLMs() function.

[10]:

# Set specific LM coeffs by list with setBLMs, items are [l,m,value]
from epsproc.sphCalc import setBLMs
BLM = setBLMs([[0,0,1],[1,1,1],[2,2,1]])   # Note different index
BLM

[11]:

# Set specific LM coeffs by numpy array with setBLMs, items are [l,m,value]
from epsproc.sphCalc import setBLMs
BLM = setBLMs(np.array([[0,0,1],[1,1,1],[2,2,1]]))   # Note different index
BLM

[12]:

# Set all allowed BLM up to Lmax with genLM
LM = ep.genLM(2)
# BLMall = setBLMs(np.ones([LM.shape[0],1]), LM)
BLMall = setBLMs(np.ones(LM.shape[0]), LM)
BLMall

[13]:

# Set all allowed BLM up to Lmax with genLM, random value case
LM = ep.genLM(2)
BLMallrand = setBLMs(np.random.random(LM.shape[0]), LM)
BLMallrand

Note that the function adds a t dimension, which is used as a generic index in these cases. For 1D cases it can be dropped by squeezing the array:

[14]:

BLMallrand = BLMallrand.squeeze(drop=True)
BLMallrand

[14]:

<xarray.DataArray 'BLM' (BLM: 9)>
0.494 0.44 0.2508 0.3413 0.629 0.1287 0.5041 0.7749 0.9317
Coordinates:
  * BLM      (BLM) object MultiIndex
  * l        (BLM) int64 0 1 1 1 2 2 2 2 2
  * m        (BLM) int64 0 -1 0 1 -2 -1 0 1 2
Attributes:
    dataType:   BLM
    long_name:  Beta parameters
    units:      arb
    harmonics:  {'dtype': 'Complex harmonics', 'kind': 'complex', 'normType':...

… or, for cases \(\beta_{L,M}(t)\), can be set directly (or renamed and used for other coords). Currently only a single dimension is supported here.

[15]:

# Set all allowed BLM up to Lmax with genLM, random value case, and as a function of t
# If coordinates are not specified, integer labels are used.
LM = ep.genLM(2)
BLMallrandt = setBLMs(np.random.random([LM.shape[0],5]), LM)
BLMallrandt

For a tabulation of the expansion parameters (via Pandas), use tabulateLM.

[16]:

from epsproc.sphFuncs.sphConv import tabulateLM
tabulateLM(BLMallrand)

[16]:

m	-2	-1	0	1	2
l
0	0.000000	0.000000	0.494049	0.000000	0.000000
1	0.000000	0.439969	0.250776	0.341302	0.000000
2	0.628958	0.128726	0.504075	0.774880	0.931717

Versions

[32]:

import scooby
scooby.Report(additional=['epsproc', 'holoviews', 'hvplot', 'xarray', 'matplotlib', 'bokeh'])

[32]:

Wed Sep 28 18:55:07 2022 UTC
OS	Linux	CPU(s)	32	Machine	x86_64	Architecture	64bit
RAM	50.1 GiB	Environment	Jupyter
Python 3.9.10 \| packaged by conda-forge \| (main, Feb 1 2022, 21:24:11) [GCC 9.4.0]
epsproc	1.3.2-dev	holoviews	1.14.8	hvplot	Module not found	xarray	2022.6.0
matplotlib	3.5.1	bokeh	2.4.2	numpy	1.21.5	scipy	1.8.0
IPython	8.1.1	scooby	0.5.12

[33]:

# Check current Git commit for local ePSproc version
from pathlib import Path
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1

* dev
  master
  numba-tests
0a15eb6cf4e469446efbf0ab2518d3fb1e64dfdb

[34]:

# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc

92c661789a7d2927f2b53d7266f57de70b3834fa        refs/heads/dependabot/pip/notes/envs/envs-versioned/mistune-2.0.3
fe1e9540c7b91fe571f60562acd31d8e489d491e        refs/heads/dependabot/pip/notes/envs/envs-versioned/nbconvert-6.5.1
0a15eb6cf4e469446efbf0ab2518d3fb1e64dfdb        refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e        refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee        refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57        refs/heads/revert-9-master

Working with spherical harmonics

Expansions in (complex) spherical harmonics

Expansions in Legendre polynomials

Expansions in symmetrized harmonics

Numerical implementation

Imports

Computing spherical harmonics on a grid

Setting coefficients

Calculating & plotting expansions

Conversions

Coefficient conversion

Generating conjugate terms

SHtools conversion

Versions