epsproc.sphCalc module

ePSproc spherical function calculations.

Collection of functions for calculating Spherical Tensors: Ylm, wignerD etc.

For spherical harmonics, currently using scipy.special.sph_harm

For other functions, using Moble’s spherical_functions package https://github.com/moble/spherical_functions

See tests/Spherical function testing Aug 2019.ipynb

04/12/19 Added setPolGeoms() to define frames as Xarray.

Added setADMs() to define ADMs as Xarray

02/12/19 Added basic TKQ multipole frame rotation routine. 27/08/19 Added wDcalc for Wigner D functions. 14/08/19 v1 Implmented sphCalc

epsproc.sphCalc.TKQarrayRot(TKQ, eAngs)[source]

Frame rotation for multipoles $T_{K,Q}$.

Basic frame rotation code, see https://github.com/phockett/Quantum-Metrology-with-Photoelectrons/blob/master/Alignment/Alignment-1.ipynb for examples.

Parameters:

  • TKQ (np.array) – Values defining the initial distribution, [K,Q,TKQ]

  • eAngs (list or np.array) – List of Euler angles (theta,phi,chi) defining rotated frame.

Returns:

  • TKQRot (np.array) – Multipoles $T’_{K,Q}$ in rotated frame, as an np.array [K,Q,TKQ].

  • TODO (redo with Moble’s functions, and Xarray input & output.)

  • Formalism

  • ———-

  • For the state multipoles, frame rotations are fairly straightforward

  • (Eqn. 4.41 in Blum)

  • .. math:: – begin{equation} leftlangle T(J’,J)_{KQ}^{dagger}rightrangle =sum_{q}leftlangle T(J’,J)_{Kq}^{dagger}rightrangle D(Omega)_{qQ}^{K*} end{equation}

  • Where $D(Omega)_{qQ}^{K}$ is a Wigner rotation operator, for a*

  • rotation defined by a set of Euler angles $Omega={theta,phi,chi}$.

  • Hence the multipoles transform, as expected, as irreducible tensors,

  • i.e. components $q$ are mixed by rotation, but terms of different

  • rank $K$ are not.

epsproc.sphCalc.TKQarrayRotX(TKQin, RX, form=2)[source]

Frame rotation for multipoles $T_{K,Q}$.

Basic frame rotation code, see https://github.com/phockett/Quantum-Metrology-with-Photoelectrons/blob/master/Alignment/Alignment-1.ipynb for examples.

Parameters:

  • TKQin (Xarray) – Values defining the initial distribution, [K,Q,TKQ]. Other dimensions will be propagated.

  • RX (Xarray defining frame rotations, from epsproc.setPolGeoms()) – List of Euler angles (theta,phi,chi) and corresponding quaternions defining rotated frame.

Returns:

TKQRot – Multipoles $T’_{K,Q}$ in rotated frame, as an np.array [K,Q,TKQ].

Return type:

Xarray

Formalism

For the state multipoles, frame rotations are fairly straightforward (Eqn. 4.41 in Blum):

\[\begin{equation} \left\langle T(J',J)_{KQ}^{\dagger}\right\rangle =\sum_{q}\left\langle T(J',J)_{Kq}^{\dagger}\right\rangle D(\Omega)_{qQ}^{K*} \end{equation}\]

Where $D(Omega)_{qQ}^{K*}$ is a Wigner rotation operator, for a rotation defined by a set of Euler angles $Omega={theta,phi,chi}$. Hence the multipoles transform, as expected, as irreducible tensors, i.e. components $q$ are mixed by rotation, but terms of different rank $K$ are not.

Examples

>>> vFlag = 2
>>> RX = ep.setPolGeoms(vFlag = vFlag)  # Package version
>>> RX
>>> testADMX = ep.setADMs(ADMs=[[0,0,0,1],[2,0,0,0.5]])
>>> testADMX
>>> testADMrot, wDX, wDXre = TKQarrayRotX(testADMX, RX)
>>> testADMrot
>>> testADMrot.attrs['dataType'] = 'ADM'
>>> sph, _ = sphFromBLMPlot(testADMrot, facetDim = 'Euler', plotFlag = True)
epsproc.sphCalc.setADMs(ADMs=[0, 0, 0, 1], KQSLabels=None, t=None, addS=False, name=None, tUnits='ps')[source]

Create Xarray from ADMs, or create default case ADM(K,Q,S) = [0,0,0,1].

Parameters:

  • ADMs (list or np.array, default = [0,0,0,1]) – Set of ADMs = [K, Q, S, ADM]. If multiple ADMs are provided per (K,Q,S) index, they are set to the t axis (if provided), or indexed numerically.

  • KQSLabels (list or np.array, optional, default = None) – If passed, assume ADMs are unabelled, and use (K,Q,S) indicies provided here.

  • t (list or np.array, optional, default = None) – If passed, use for dimension defining ADM sets (usually time). Defaults to numerical label if not passed, t = np.arange(0,ADMs.shape[1])

  • addS (bool, default = False) – If set, append S = 0 to ADMs. This allows for passing of [K,Q,ADM] type values (e.g. for symmetric top case)

  • name (str, optional, default = None) – Set a name for the array. If None, will be set to ‘ADM’ (same as dataType attrib)

  • tUnits (str, optional, default = 'ps') – Units for temporal axis, if set.

Returns:

ADMX – ADMs in Xarray format, dims as per epsproc.utils.ADMdimList()

Return type:

Xarray

Examples

>>> # Default case
>>> ADMX = setADMs()
>>> ADMX

>>> # Set with ranges (as an array [K,Q,S, t(0), t(1)...]), with values per t
>>> tPoints = 10
>>> ADMX = setADMs(ADMs = [[0,0,0, *np.ones(10)], [2,0,0, *np.linspace(0,1,tPoints)], [4,0,0, *np.linspace(0,0.5,tPoints)]])

>>> # With full N2 rotational wavepacket ADM set from demo data (ePSproc\datalignment), where modPath defines root...
>>> # Load ADMs for N2
>>> from scipy.io import loadmat
>>> ADMdataFile = os.path.join(modPath, 'data', 'alignment', 'N2_ADM_VM_290816.mat')
>>> ADMs = loadmat(ADMdataFile)
>>> ADMX = setADMs(ADMs = ADMs['ADM'], KQSLabels = ADMs['ADMlist'], addS = True)
>>> ADMX
epsproc.sphCalc.setBLMs(BLMs=[0, 0, 1], LMLabels=None, t=None, name=None, tUnits='ps', conformDims=False, **kwargs)[source]

Create Xarray from BLMs, or create default case BLM = [0,0,1].

Parameters:

  • BLMs (list or np.array, default = [0,0,1]) – Set of BLMs = [L, M, BLM]. If multiple BLMs are provided per (L,M) index, they are set to the E or t axis (if provided), or indexed numerically.

  • LMLabels (list or np.array, optional, default = None) – If passed, assume BLMs are unabelled, and use (L,M) indicies provided here.

  • t (list or np.array, optional, default = None) – If passed, use for dimension defining BLM sets (usually time or energy). Defaults to numerical label if not passed, t = np.arange(0,ADMs.shape[1])

  • name (str, optional, default = None) – Set a name for the array. If None, will be set to ‘ADM’ (same as dataType attrib)

  • tUnits (str, optional, default = 'ps') – Units for t axis, if set.

  • conformDims (bool, optional, default = False) – Add any missing dims to match ep.listFuncs.BLMdimList. NOT YET IMPLEMENTED - see PEMtk.toePSproc() for method.

Returns:

  • BLMX (Xarray) – BLMs in Xarray format, dims as per epsproc.utils.BLMdimList()

  • TODO

  • - Implemnt conformDims, see PEMtk.toePSproc() for method.

  • - General handling for additional dims, currently only ‘t’ supported. Should change to dict mapping.

Examples

>>> # Default case
>>> BLMX = setBLMs()
>>> BLMX

>>> # Set with ranges (as an array [L,M, t(0), t(1)...]), with values per t
>>> tPoints = 10
>>> BLMX = setBLMs(ADMs = [[0,0, *np.ones(10)], [2,0, *np.linspace(0,1,tPoints)], [4,0, *np.linspace(0,0.5,tPoints)]])
epsproc.sphCalc.setPolGeoms(eulerAngs=None, quat=None, labels=None, vFlag=2, defaultMap='canonical')[source]

Generate Xarray containing polarization geometries as Euler angles and corresponding quaternions.

Define LF > MF polarization geometry/rotations. Provide either eulerAngs or quaternions, but not both (supplied quaternions only will be used in this case).

For default case (eulerAngs = None, quat = None), 3 geometries are calculated, corresponding to z-pol, x-pol and y-pol cases. Defined by Euler angles: (p,t,c) = [0 0 0] for z-pol, (p,t,c) = [0 pi/2 0] for x-pol, (p,t,c) = [pi/2 pi/2 0] for y-pol.

Parameters:

  • eulerAngs (list or np.array of Euler angles (p(hi), t(heta), c(hi)), optional.) – List or array [p,t,c…], shape (Nx3). List or array including set labels, [label,p,t,c…], shape (Nx4)

  • quat (list or np.array of quaternions, optional.) –

  • labels (list of labels, one per set of angles. Optional.) – If not set, states will be labelled numerically.

  • vFlag (version of routine to use, optional, default = 2) – Options: - 1, use labels as sub-dimensional coord. - 2, set labels as non-dimensional coord.

  • defaultMap (str, optional, default = 'canonical') –

    Mapping to use for default case. - ‘canonical’ as listed above.

    Use this in most cases.

    • ’exy’ for use with ep.efield.epol.EfieldPol() for conversion

      of (Ex,Ey) or (El,Er) basis to canonical orientation. This assumes x-axis is main pol axis, and maps x axis > z axis for z-pol case. Note this is appropriate for linear pol light, but less so for circ pol.

    • ’exyDiag’

      As ‘exy’, plus some diagonal cases.

    • ’circ’ same as ‘canoncial’.

      For circ. pol. canonical case with (Ex,Ey) or (El,Er) basis. Here (z,x,y) correspond to propagation axis.

Returns:

  • RX (Xarray of quaternions, with Euler angles as dimensional params.)

  • To do

  • —–

  • - Better label handling, as dictionary? With mixed-type array may get issues later. – (sf.quaternion doesn’t seem to have an issue however.)

  • - Xarray MultiIndex with mixed types? – Tested with pd - not supported: >>> eulerInd = pd.MultiIndex.from_arrays([eulerAngs[:,0].T, eulerAngs[:,1:].T.astype(‘float’)], names = [‘Label’,’P’,’T’,’C’]) # Gives error: # NotImplementedError: > 1 ndim Categorical are not supported at this time

Examples

>>> # Defaults
>>> RXdefault = setPolGeoms()
>>> print(RXdefault)

>>> # Pass Eulers, no labels
>>> pRot = [1.666, 0, np.pi/2]
>>> tRot = [0, np.pi/2, np.pi/2]
>>> cRot = [-1.5, 0, 0]

>>> eulerAngs = np.array([pRot, tRot, cRot]).T

>>> RXePass = setPolGeoms(eulerAngs = eulerAngs)
>>> print(RXePass)

>>> # Pass labels separately
>>> RXePass = setPolGeoms(eulerAngs = eulerAngs, labels = ['1','23','ff'])
>>> print(RXePass)

>>> # Pass Eulers with existing labels
>>> labels = ['A','B','C']
>>> eulerAngs = np.array([labels, pRot, tRot, cRot]).T
>>> RXePass = setPolGeoms(eulerAngs = eulerAngs)
>>> print(RXePass)

>>> # Pass Quaternions and labels
>>> RXqPass = setPolGeoms(quat = RXePass, labels = labels)
>>> print(RXqPass)

>>> # Pass both - only quaternions will be used in this case, and warning displayed.
>>> RXqeTest = setPolGeoms(eulerAngs = eulerAngs, quat = RXePass, labels = labels)
>>> print(RXqeTest)
epsproc.sphCalc.sphCalc(Lmax=2, Lmin=0, res=None, angs=None, XFlag=True, fnType='sph', convention='phys', conj=False, realCSphase=False)[source]

Calculate set of spherical harmonics Ylm(theta,phi) on a grid.

Parameters:

  • Lmax (int) – Maximum L for the set. Ylm calculated for Lmin:Lmax, all m.

  • Lmin (int, optional, default 0) – Min L for the set. Ylm calculated for Lmin:Lmax, all m.

  • res (int or list, optional, default None) – (Theta, Phi) grid resolution, outputs will be of dim [res,res] (if int) or [res[0], res[1]] (if list).

  • angs (list of 2D np.arrays, [thetea, phi], optional, default None) – If passed, use these grids for calculation. NOTE: if ‘maths’ convention set, this array will be assumed to be [phi, theta]

  • XFlag (bool, optional, default True) – Flag for output. If true, output is Xarray. If false, np.arrays

  • fnType (str, optional, default = 'sph') – Currently can set to ‘sph’ for SciPy spherical harmonics, or ‘lg’ for SciPy Legendre polynomials. 23/09/22 - hacked in ‘complex’ and ‘real’ here too, but TODO: split by kind and backend. More backends to follow.

  • convention (str, optional, default = 'phys') – Set to ‘phys’ (theta from z-axis) or ‘maths’ (phi from z-axis) spherical polar conventions. (Might need some work!)

  • conj (bool, optional, default = False) – Return complex conjugate.

  • realCSphase (bool, optional, default = False) – Apply (-1)^m term in real harmonics calc. Generally will already be included in complex spherical harmonics defn., so can be left as False.

Note that either res OR angs needs to be passed.

Outputs

  • if XFlag -

YlmX

3D Xarray, dims (lm,theta,phi)

  • else -

Ylm, lm

3D np.array of values, dims (lm,theta,phi), plus list of lm pairs

Currently set for scipy.special.sph_harm as calculation routine. Note (theta, phi) definition, and normalisation.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
.. math::

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

And for real harmonics:
.. math::

begin{equation} begin{aligned} Y_{ell m}&={begin{cases}{sqrt {2}},Im [{Y_{ell }^{|m|}}]&{text{if}}mlt0 \Y_{ell }^{0}&{text{if}}m=0 \{sqrt {2}},Re [{Y_{ell }^{m}}]&{text{if}}mgt0.end{cases}} end{aligned} end{equation}

See https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form.
Demos:
- Complex: https://epsproc.readthedocs.io/en/dev/special_topics/ePSproc_docs_working_with_spherical_harmonics_200922.html
- Real: https://epsproc.readthedocs.io/en/dev/special_topics/ePSproc_docs_working_with_real_harmonics_220922.html
- For scipy backend, realCSphase = False case matches SHtools defn, see https://epsproc.readthedocs.io/en/dev/special_topics/ePSproc_docs_working_with_real_harmonics_220922.html#Converting-coefficients-real-to-complex-form

Example

>>> YlmX = sphCalc(2, res = 50)

Notes

TODO:

  • additional backends.

  • remove hard-coded dim names for flexibility (or set from ref YLMdimList())

  • remove ‘fnType’ and/or use ‘kind’ instead (SHtool notation)

epsproc.sphCalc.wDcalc(Lrange=[0, 1], Nangs=None, eAngs=None, R=None, XFlag=True, QNs=None, dlist=['lp', 'mu', 'mu0'], eNames=['P', 'T', 'C'], conjFlag=False, sfError=True, verbose=False)[source]

Calculate set of Wigner D functions D(l,m,mp; R) on a grid.

Parameters:

  • Lrange (list, optional, default [0, 1]) – Range of L to calculate parameters for. If len(Lrange) == 2 assumed to be of form [Lmin, Lmax], otherwise list is used directly. For a given l, all (m, mp) combinations are calculated.

  • QNs (np.array, optional, default = None) – List of QNs [l,m,mp] to compute Wigner D terms for. If supplied, use this instead of Lrange setting.

  • only) (Options for setting angles (use one) –

  • Nangs (int, optional, default None) – If passed, use this to define Euler angles sampled. Ranges will be set as (theta, phi, chi) = (0:pi, 0:pi/2, 0:pi) in Nangs steps.

  • eAngs (np.array, optional, default None) – If passed, use this to define Euler angles sampled. Array of angles, [theta,phi,chi], in radians

  • R (np.array, optional, default None) – If passed, use this to define Euler angles sampled. Array of quaternions, as given by quaternion.from_euler_angles(eAngs).

XFlagbool, optional, default True

Flag for output. If true, output is Xarray. If false, np.arrays

dlistlist, optional, default [‘lp’,’mu’,’mu0’]

Labels for Xarray QN dims.

eNameslist, optional, default [‘P’,’T’,’C’]

Labels for Xarray Euler dims.

conjFlagbool, optional, default = False

If true, return complex conjuage values.

sfErrorbool, optional, default = None

If not None, set sf.error_on_bad_indices = sfError If True (default case), this will raise a value error on bad QNs. If False, set = 0. See code at https://github.com/moble/spherical_functions/blob/master/spherical_functions/WignerD/__init__.py 05/05/21 - added, but CURRENTLY NOT WORKING

Outputs

  • if XFlag -

wDX

Xarray, dims (lmmp,Euler)

  • else -

wD, R, lmmp

np.arrays of values, dims (lmmp,Euler), plus list of angles and lmmp sets.

Uses Moble's spherical_functions package for wigner D function.
https://github.com/moble/spherical_functions
Moble's quaternion package for angles and conversions.
https://github.com/moble/quaternion
For testing, see https://epsproc.readthedocs.io/en/latest/tests/Spherical_function_testing_Aug_2019.html

Examples

>>> wDX1 = wDcalc(eAngs = np.array([0,0,0]))

>>> wDX2 = wDcalc(Nangs = 10)