Method development for geometric functions
26/02/20
Aims:
Develop \(\beta_{L,M}\) formalism.
Develop corresponding numerical methods.
Speed things up (see low-level benchmarking notebook).
Analyse geometric terms.
Setup
[50]:
# Imports
import numpy as np
import pandas as pd
import xarray as xr
from functools import lru_cache # For function result caching
# Special functions
# from scipy.special import sph_harm
import spherical_functions as sf
import quaternion
# Performance & benchmarking libraries
# from joblib import Memory
# import xyzpy as xyz
import numba as nb
# Timings with ttictoc
# https://github.com/hector-sab/ttictoc
# from ttictoc import TicToc
# Package fns.
# For module testing, include path to module here
import sys
import os
modPath = r'D:\code\github\ePSproc' # Win test machine
# modPath = r'/home/femtolab/github/ePSproc/' # Linux test machine
sys.path.append(modPath)
import epsproc as ep
# TODO: tidy this up!
from epsproc.util import matEleSelector
from epsproc.geomFunc import geomCalc
Exploring Wigner 3js
In photoionization calculations, there is a lot of angular momentum coupling to deal with. Typically, 4 to 6 Wigner 3j terms appear (depending on the formalism), and/or higher-order terms in cases where couplings are included.
\begin{equation}
W = \left(\begin{array}{ccc}
l & l' & L\\
m & m' & M
\end{array}\right)
\end{equation}
Since this is, effectively, a 6D space, dimensions \((l_{max}, l_{max}, 2l_{max}, 2l_{max}+1, 2l_{max}+1, 4l_{max}+1)\) things can get large quickly. For small \(l_{max}\) it’s easy to look at some values directly…
(For more details on 3j symbols, see Wikipedia and Wolfram Mathworld; for more on the numerics, see the test notebook, and benchmarks.)
TODO: test numerically some of the symmetry properties here, by summation.
[ ]:
[109]:
# Calculate some values.
# w3jTable will output all values up to l=lp=Lmax (hence L=2Lmax)
lmax = 1
w3jlist = geomCalc.w3jTable(Lmax = lmax, form = '2d') # For form = '2d', the function will output only valid entries as a coordinate table
print(w3jlist.shape)
print(f'Max value: {w3jlist[:,-1].max()}, min value: {w3jlist[:,-1].min()}\n')
# Print the table - output format has rows (l, lp, L, m, mp, M, 3j)
print(w3jlist)
(26, 7)
Max value: 1.0, min value: -0.5773502691896257
[[ 0. 0. 0. 0. 0. 0.
1. ]
[ 0. 1. 1. 0. -1. 1.
0.57735027]
[ 0. 1. 1. 0. 0. 0.
-0.57735027]
[ 0. 1. 1. 0. 1. -1.
0.57735027]
[ 1. 0. 1. -1. 0. 1.
0.57735027]
[ 1. 0. 1. 0. 0. 0.
-0.57735027]
[ 1. 0. 1. 1. 0. -1.
0.57735027]
[ 1. 1. 2. -1. -1. 2.
0.4472136 ]
[ 1. 1. 1. -1. 0. 1.
0.40824829]
[ 1. 1. 2. -1. 0. 1.
-0.31622777]
[ 1. 1. 0. -1. 1. 0.
0.57735027]
[ 1. 1. 1. -1. 1. 0.
-0.40824829]
[ 1. 1. 2. -1. 1. 0.
0.18257419]
[ 1. 1. 1. 0. -1. 1.
-0.40824829]
[ 1. 1. 2. 0. -1. 1.
-0.31622777]
[ 1. 1. 0. 0. 0. 0.
-0.57735027]
[ 1. 1. 1. 0. 0. 0.
0. ]
[ 1. 1. 2. 0. 0. 0.
0.36514837]
[ 1. 1. 1. 0. 1. -1.
0.40824829]
[ 1. 1. 2. 0. 1. -1.
-0.31622777]
[ 1. 1. 0. 1. -1. 0.
0.57735027]
[ 1. 1. 1. 1. -1. 0.
0.40824829]
[ 1. 1. 2. 1. -1. 0.
0.18257419]
[ 1. 1. 1. 1. 0. -1.
-0.40824829]
[ 1. 1. 2. 1. 0. -1.
-0.31622777]
[ 1. 1. 2. 1. 1. -2.
0.4472136 ]]
[110]:
# Recalculate and set to Xarray output format, then plot with ep.lmPlot()
w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')
# Check number of valid entries matches basic table above
print(f'Number of valid (non-NaN) elements: {w3j.count()}')
# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
Number of valid (non-NaN) elements: <xarray.DataArray 'w3jStacked' ()>
array(26)
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
[111]:
# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[111]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221eedaca20>
[53]:
# Print out values by QNs (Pandas table)
daPlotpd
[53]:
| L | 0 | 1 | 2 | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| M | 0 | -1 | 0 | 1 | -2 | -1 | 0 | 1 | 2 | |||
| l | m | lp | mp | |||||||||
| 0 | 0 | 0 | 0 | 1.00000 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | -1 | NaN | NaN | NaN | 0.577350 | NaN | NaN | NaN | NaN | NaN | ||
| 0 | NaN | NaN | -0.577350 | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | NaN | 0.577350 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | -1 | 0 | 0 | NaN | NaN | NaN | 0.577350 | NaN | NaN | NaN | NaN | NaN |
| 1 | -1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.447214 | ||
| 0 | NaN | NaN | NaN | 0.408248 | NaN | NaN | NaN | -0.316228 | NaN | |||
| 1 | 0.57735 | NaN | -0.408248 | NaN | NaN | NaN | 0.182574 | NaN | NaN | |||
| 0 | 0 | 0 | NaN | NaN | -0.577350 | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | -1 | NaN | NaN | NaN | -0.408248 | NaN | NaN | NaN | -0.316228 | NaN | ||
| 0 | -0.57735 | NaN | 0.000000 | NaN | NaN | NaN | 0.365148 | NaN | NaN | |||
| 1 | NaN | 0.408248 | NaN | NaN | NaN | -0.316228 | NaN | NaN | NaN | |||
| 1 | 0 | 0 | NaN | 0.577350 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | -1 | 0.57735 | NaN | 0.408248 | NaN | NaN | NaN | 0.182574 | NaN | NaN | ||
| 0 | NaN | -0.408248 | NaN | NaN | NaN | -0.316228 | NaN | NaN | NaN | |||
| 1 | NaN | NaN | NaN | NaN | 0.447214 | NaN | NaN | NaN | NaN | |||
This ends up as a relatively sparse array, since many combinations are invalid (do not follow angular momentum selection rules), hence there are many NaN terms.
The results can also be output as a 6D sparse array, using the Sparse library.
[54]:
# Calculate and output in Sparse array format
w3jSparse = geomCalc.w3jTable(Lmax = lmax, form = 'ndsparse')
w3jSparse
[54]:
| Format | coo |
|---|---|
| Data Type | float64 |
| Shape | (2, 2, 3, 2, 2, 3) |
| nnz | 25 |
| Density | 0.1736111111111111 |
| Read-only | True |
| Size | 1.4K |
| Storage ratio | 1.2 |
Here nnz is the number of non-zero elements.
[55]:
# Try a larger Lmax and plot only.
lmax = 3
w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')
# Check number of valid entries matches basic table above
print(f'Number of valid (non-NaN) elements: {w3j.count()}')
# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
Number of valid (non-NaN) elements: <xarray.DataArray 'w3jStacked' ()>
array(820)
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
[56]:
# Resort axis by (l,lp)
plotDimsRed = ['l', 'lp', 'm', 'mp']
xDim = {'LM':['L','M']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
[57]:
# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[57]:
<epsproc._sns_matrixMod.ClusterGrid at 0x22183e4dfd0>
This clearly shows that the valid terms become sparser at higher \(l\), and the couplings become smaller.
Structure can also be examined using other methods, e.g. correlation functions (see, for example, Seaborn Discovering structure in heatmap data. The example here shows Panda’s standard Pearson correlation coefficient, which may (or may not) be particularly meaningful here… but does show structures.
[58]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[58]:
<epsproc._sns_matrixMod.ClusterGrid at 0x22183f526a0>
[59]:
# Test some other Seaborn methods... these likely won't scale well for large lmax!
import seaborn as sns
# Recalculate for small lmax
lmax = 1
w3j = geomCalc.w3jTable(Lmax = lmax, form = 'xdaLM')
# Set parameters to restack the Xarray into (L,M) pairs
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(w3j, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
# Try sns pairplot
# sns.pairplot(daPlotpd.fillna(0).T) # Big grids!
# sns.pairplot(daPlotpd.fillna(0)) # OK, not particularly informative
# sns.pairplot(daPlotpd.fillna(0), hue = 'l') # Doesn't work - multindex issue?
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
\(E_{P,R}\) tensor
The coupling of two 1-photon terms can be written as a tensor contraction:
\begin{equation}
E_{PR}(\hat{e})=[e\otimes e^{*}]_{R}^{P}=[P]^{\frac{1}{2}}\sum_{p}(-1)^{R}\left(\begin{array}{ccc}
1 & 1 & P\\
p & R-p & -R
\end{array}\right)e_{p}e_{R-p}^{*}\label{eq:EPR-defn-1}
\end{equation}
Where \(e_{p}\) and \(e_{R-p}\) define the field strengths for the polarizations \(p\) and \(R-p\), which are coupled into the spherical tensor \(E_{PR}\).
[60]:
# Calculate EPR terms, all QNs, with field strengths e = 1
EPRtable = geomCalc.EPR(form = '2d')
# Output values as list, [l, lp, P, p, R-p, R, EPR]
print(EPRtable)
[[ 1. 1. 0. -1. 1. 0.
0.57735027]
[ 1. 1. 1. -1. 0. 1.
-0.70710678]
[ 1. 1. 1. -1. 1. 0.
-0.70710678]
[ 1. 1. 2. -1. -1. 2.
1. ]
[ 1. 1. 2. -1. 0. 1.
0.70710678]
[ 1. 1. 2. -1. 1. 0.
0.40824829]
[ 1. 1. 0. 0. 0. 0.
-0.57735027]
[ 1. 1. 1. 0. -1. 1.
0.70710678]
[ 1. 1. 1. 0. 1. -1.
-0.70710678]
[ 1. 1. 2. 0. -1. 1.
0.70710678]
[ 1. 1. 2. 0. 0. 0.
0.81649658]
[ 1. 1. 2. 0. 1. -1.
0.70710678]
[ 1. 1. 0. 1. -1. 0.
0.57735027]
[ 1. 1. 1. 1. -1. 0.
0.70710678]
[ 1. 1. 1. 1. 0. -1.
0.70710678]
[ 1. 1. 2. 1. -1. 0.
0.40824829]
[ 1. 1. 2. 1. 0. -1.
0.70710678]
[ 1. 1. 2. 1. 1. -2.
1. ]]
As before, we can visualise these values…
[61]:
EPRX = geomCalc.EPR(form = 'xarray')
# Set parameters to restack the Xarray into (L,M) pairs
# plotDimsRed = ['l', 'p', 'lp', 'R-p']
plotDimsRed = ['p', 'R-p']
xDim = {'PR':['P','R']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(EPRX, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
# Seem to have some all-NaN cols persisting here, not sure why...
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
No handles with labels found to put in legend.
[62]:
# Tabulate
daPlotpd.dropna(axis = 1, how = 'all')
[62]:
| P | 0 | 1 | 2 | |||||||
|---|---|---|---|---|---|---|---|---|---|---|
| R | 0 | -1 | 0 | 1 | -2 | -1 | 0 | 1 | 2 | |
| p | R-p | |||||||||
| -1 | -1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 1.0 |
| 0 | NaN | NaN | NaN | -0.707107 | NaN | NaN | NaN | 0.707107 | NaN | |
| 1 | 0.57735 | NaN | -0.707107 | NaN | NaN | NaN | 0.408248 | NaN | NaN | |
| 0 | -1 | NaN | NaN | NaN | 0.707107 | NaN | NaN | NaN | 0.707107 | NaN |
| 0 | -0.57735 | NaN | NaN | NaN | NaN | NaN | 0.816497 | NaN | NaN | |
| 1 | NaN | -0.707107 | NaN | NaN | NaN | 0.707107 | NaN | NaN | NaN | |
| 1 | -1 | 0.57735 | NaN | 0.707107 | NaN | NaN | NaN | 0.408248 | NaN | NaN |
| 0 | NaN | 0.707107 | NaN | NaN | NaN | 0.707107 | NaN | NaN | NaN | |
| 1 | NaN | NaN | NaN | NaN | 1.0 | NaN | NaN | NaN | NaN | |
[63]:
# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[63]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221812b0ac8>
[64]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[64]:
<epsproc._sns_matrixMod.ClusterGrid at 0x2218055fdd8>
\(B_{L,M}\) term
The coupling of the partial wave pairs, \(|l,m\rangle\) and \(|l',m'\rangle\), into the observable set of \(\{L,M\}\) is defined by a tensor contraction with two 3j terms.
\begin{equation}
B_{L,M}=(-1)^{m}\left(\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc}
l & l' & L\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
l & l' & L\\
-m & m' & M
\end{array}\right)
\end{equation}
(See notebook ePSproc_wigner3j_dataStructures_260220.ipynb on Bemo for additonal dev details)
Note that this term is equivalent, effectively, as the triple intergral over spherical harmonics (e.g. eqn. 3.119 in Zare):
\begin{equation}
\intop_{0}^{2\pi}\intop_{0}^{\pi}Y_{J_{3}M_{3}}(\theta,\phi)Y_{J_{2}M_{2}}(\theta,\phi)Y_{J_{1}M_{1}}(\theta,\phi)\sin\theta d\theta d\phi=\left(\frac{(2J_{1}+1)(2J_{2}+1)(2J_{3}+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc}
J_{1} & J_{2} & J_{3}\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
J_{1} & J_{2} & J_{3}\\
M_{1} & M_{2} & M_{3}
\end{array}\right)
\end{equation}
… or contraction over a pair of harmonics into a resultant harmonic (e.g. eqns. C.21, C.22 in Blum), and this is how the term arises in the derivation.
\begin{equation}
Y_{J_{1}M_{1}}(\theta,\phi)Y_{J_{2}M_{2}}(\theta,\phi)=\sum_{J_{3}M_{3}}\left(\frac{(2J_{1}+1)(2J_{2}+1)(2J_{3}+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc}
J_{1} & J_{2} & J_{3}\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
J_{1} & J_{2} & J_{3}\\
M_{1} & M_{2} & M_{3}
\end{array}\right)Y^*_{J_{3}M_{3}}(\theta,\phi)
\end{equation}
(Again, for more details on see Wikipedia and Wolfram Mathworld.)
Note also some defns. use conjugate spherical harmonics, can convert as (Blum C.21):
\begin{equation}
\beta_{L,M}^{\mu_{i},\mu_{f}}Y_{LM}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})=\beta_{L,-M}^{\mu_{i},\mu_{f}}(-1)^{M}Y_{L,-M}(\theta_{\hat{k}},\phi_{\hat{k}})\label{eq:sph-conj-conv}
\end{equation}
[65]:
Lmax = 1
BLMtable = geomCalc.betaTerm(Lmax = Lmax, form = 'xdaLM') # Output as stacked Xarray
[66]:
BLMtable
[66]:
<xarray.DataArray 'w3jStacked' (lSet: 6, mSet: 9)>
array([[ nan, nan, nan, nan, 0.28209479,
nan, nan, nan, nan],
[ nan, nan, nan, -0.28209479, 0.28209479,
-0.28209479, nan, nan, nan],
[ nan, 0.28209479, nan, nan, 0.28209479,
nan, nan, 0.28209479, nan],
[ nan, nan, 0.28209479, nan, 0.28209479,
nan, 0.28209479, nan, nan],
[ nan, nan, nan, nan, nan,
nan, nan, nan, nan],
[-0.30901936, 0.21850969, -0.12615663, -0.21850969, 0.25231325,
-0.21850969, -0.12615663, 0.21850969, -0.30901936]])
Coordinates:
* mSet (mSet) MultiIndex
- m (mSet) int64 -1 -1 -1 0 0 0 1 1 1
- mp (mSet) int64 -1 0 1 -1 0 1 -1 0 1
- M (mSet) int64 2 1 0 1 0 -1 0 -1 -2
* lSet (lSet) MultiIndex
- l (lSet) int64 0 0 1 1 1 1
- lp (lSet) int64 0 1 0 1 1 1
- L (lSet) int64 0 1 1 0 1 2
Attributes:
dataType: Wigner3j[67]:
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
# daPlot = BLMtable.unstack().stack({'LM':['L','M']})
# ep.lmPlot(w3jXcombMult.unstack().stack({'LM':['l','m']}), plotDims=['lp', 'L', 'mp', 'M'], xDim='LM', SFflag = False)
# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable.unstack().stack(xDim), plotDims=plotDimsRed, xDim=xDim, SFflag = False, squeeze = False)
# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim)
# ep.lmPlot(w3jXcombMult.unstack(), plotDims=['lp', 'L', 'mp', 'M'], xDim='L', SFflag = False)
# daPlotpd = daPlot.unstack().stack(plotDim = plotDimsRed).to_pandas().dropna(axis = 1).T
[68]:
# daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable.unstack().stack(xDim), plotDims=plotDimsRed, xDim=xDim, pType = 'r')
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
[69]:
daPlotpd
[69]:
| L | 0 | 1 | 2 | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| M | 0 | -1 | 0 | 1 | -2 | -1 | 0 | 1 | 2 | |||
| l | m | lp | mp | |||||||||
| 0 | 0 | 0 | 0 | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | -1 | NaN | NaN | NaN | -0.282095 | NaN | NaN | NaN | NaN | NaN | ||
| 0 | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | NaN | -0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | -1 | 0 | 0 | NaN | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN |
| 1 | -1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | -0.309019 | ||
| 0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.21851 | NaN | |||
| 1 | 0.282095 | NaN | NaN | NaN | NaN | NaN | -0.126157 | NaN | NaN | |||
| 0 | 0 | 0 | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | -1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | -0.21851 | NaN | ||
| 0 | 0.282095 | NaN | NaN | NaN | NaN | NaN | 0.252313 | NaN | NaN | |||
| 1 | NaN | NaN | NaN | NaN | NaN | -0.21851 | NaN | NaN | NaN | |||
| 1 | 0 | 0 | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | -1 | 0.282095 | NaN | NaN | NaN | NaN | NaN | -0.126157 | NaN | NaN | ||
| 0 | NaN | NaN | NaN | NaN | NaN | 0.21851 | NaN | NaN | NaN | |||
| 1 | NaN | NaN | NaN | NaN | -0.309019 | NaN | NaN | NaN | NaN | |||
[70]:
# A complementary visulization is to call directly the sns.clustermap plot, use clustering and plot by category labels - see https://seaborn.pydata.org/index.html
# (ep.lpPlot uses a modified version of this routine.)
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T, center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[70]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221eed13a20>
[71]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[71]:
<epsproc._sns_matrixMod.ClusterGrid at 0x22183ac17b8>
[72]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr(), center=0, cmap="vlag", row_cluster=True, col_cluster=True)
[72]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221839db400>
What do we learn here…? This is a nice way to visualize the selection rules into the observable: for instance, only terms \(l=l'\) and \(m=-m'\) contribute to the overall photoinoization cross-section term (\(L=0, M=0\)). However, since these terms are fairly simply followed algebraically in this case, via the rules inherent in the 3j product, this is not particularly insightful. These visualizations will become more useful when dealing with real sets of matrix elements, and specific polarization geometries, which will further modulate the \(B_{L,M}\) terms.
[73]:
# PLOT AS A FN. of dl, dm... should just need to recalc dims or add a separate label for this?
# dl = BLMtable.l - BLMtable.lp
# dm = BLMtable.m - BLMtable.mp
BLMdeltas = BLMtable.unstack().stack({'LM':['L','M'], 'dl':['l','lp'], 'dm':['m','mp']})
dl = BLMdeltas.l - BLMdeltas.lp
dm = BLMdeltas.m - BLMdeltas.mp
# Set new dims & plot
BLMdeltas = BLMdeltas.assign_coords(dl=dl, dm=dm) # This works, but not plotting correctly for some reason - maybe dim names?
plotDimsRed = ['dl','dm']
# BLMdeltas = BLMdeltas.assign_coords(dl=dl, dm=dm).rename({'dl':'l', 'dm':'m'}) # Hmm, same result - values present in daPlotpd not showing up here, or losing correlations...?
# plotDimsRed = ['l','m']
xDim = {'LM':['L','M']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMdeltas, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
No handles with labels found to put in legend.
[74]:
daPlotpd
[74]:
| L | 0 | 1 | 2 | |||||||
|---|---|---|---|---|---|---|---|---|---|---|
| M | 0 | -1 | 0 | 1 | -2 | -1 | 0 | 1 | 2 | |
| dl | dm | |||||||||
| 0 | 0 | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| -1 | 1 | NaN | NaN | NaN | -0.282095 | NaN | NaN | NaN | NaN | NaN |
| 0 | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | |
| -1 | NaN | -0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | -1 | NaN | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN |
| 0 | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | |
| 1 | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| 0 | 0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | -0.309019 |
| -1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.21851 | NaN | |
| -2 | 0.282095 | NaN | NaN | NaN | NaN | NaN | -0.126157 | NaN | NaN | |
| 1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | -0.21851 | NaN | |
| 0 | 0.282095 | NaN | NaN | NaN | NaN | NaN | 0.252313 | NaN | NaN | |
| -1 | NaN | NaN | NaN | NaN | NaN | -0.21851 | NaN | NaN | NaN | |
| 2 | 0.282095 | NaN | NaN | NaN | NaN | NaN | -0.126157 | NaN | NaN | |
| 1 | NaN | NaN | NaN | NaN | NaN | 0.21851 | NaN | NaN | NaN | |
| 0 | NaN | NaN | NaN | NaN | -0.309019 | NaN | NaN | NaN | NaN | |
[75]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0))
[75]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221840ae048>
Note correlation of terms, with final \(L\) and \(M\) odd or even terms grouped according to odd or even \(\Delta l\) and \(\Delta m\) respectively.
[76]:
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr())
[76]:
<epsproc._sns_matrixMod.ClusterGrid at 0x22181071d68>
[79]:
# Test correlation fns.
# This fails with NaNs present it seems
# ep.snsMatMod.clustermap(daPlot.dropna(dim='plotDim', how='all').dropna(dim='LM',how='all').fillna(0).to_pandas().T.corr())
# ep.snsMatMod.clustermap(daPlot.dropna(dim='plotDim', how='all').dropna(dim='LM',how='all').fillna(0).to_pandas().corr())
ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr())
ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr())
[79]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221817b2828>
[80]:
# Switch plotting dims
plotDimsRed = ['L','M']
xDim = {'llpmmp':['l', 'm', 'lp', 'mp']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, mMax = 2, pType = 'r')
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
No handles with labels found to put in legend.
\(\Lambda\) Term
Define MF projection term, \(\Lambda_{R',R}(R_{\hat{n}})\):
\begin{equation}
\Lambda_{R',R}(R_{\hat{n}})=(-1)^{(R')}\left(\begin{array}{ccc}
1 & 1 & P\\
\mu & -\mu' & R'
\end{array}\right)D_{-R',-R}^{P}(R_{\hat{n}})
\end{equation}
This is similar to the \(E_{PR}\) term, and essentially rotates it into the MF.
[81]:
# Calculate MF projection terms, all QNs, for Euler angles corresponding to (z,x,y) cases
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = '2d')
# Output values as list, ['l', 'lp', 'P', 'mu', 'mup', 'Rp', 'R']
# print(lambdaTerm.shape)
# print(np.c_[QNs[0:10,:], lambdaTerm[0:10,:]])
[82]:
# Xarray form
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = 'xarray')
[83]:
# Check values Xarray vs. 2d case for consistency
# TO DO
[84]:
# In this case the QNs are left unstacked to facilitate multiplication across different terms - this fails for stacked (multindex) cases here.
lambdaTerm.dims
[84]:
('Rp', 'l', 'lp', 'P', 'mu', 'mup', 'R', 'Euler')
[85]:
# Plot all QNs vs. Euler angles (by label)
plotDimsRed = ['l', 'lp', 'P', 'mu', 'mup', 'R', 'Rp']
xDim = 'Labels'
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'a') # Plot abs values
daPlotPhase, daPlotpdPhase, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'phase') # Plot phase
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
Plotting data (No filename), pType=phase, thres=0.01, with Seaborn
[86]:
daPlotpd
[86]:
| Labels | z | x | y | ||||||
|---|---|---|---|---|---|---|---|---|---|
| l | lp | P | mu | mup | R | Rp | |||
| 1 | 1 | 0 | -1 | 1 | 0 | 0 | 0.577350 | 0.577350 | 0.577350 |
| 0 | 0 | 0 | 0 | 0.577350 | 0.577350 | 0.577350 | |||
| 1 | -1 | 0 | 0 | 0.577350 | 0.577350 | 0.577350 | |||
| 1 | -1 | 0 | -1 | 1 | 0.000000 | 0.204124 | 0.204124 | ||
| 0 | 1 | 0.000000 | 0.288675 | 0.288675 | |||||
| 1 | 1 | 0.408248 | 0.204124 | 0.204124 | |||||
| 1 | -1 | 0 | 0.000000 | 0.288675 | 0.288675 | ||||
| 0 | 0 | 0.408248 | 0.000000 | 0.000000 | |||||
| 1 | 0 | 0.000000 | 0.288675 | 0.288675 | |||||
| 0 | -1 | -1 | 1 | 0.000000 | 0.204124 | 0.204124 | |||
| 0 | 1 | 0.000000 | 0.288675 | 0.288675 | |||||
| 1 | 1 | 0.408248 | 0.204124 | 0.204124 | |||||
| 1 | -1 | -1 | 0.408248 | 0.204124 | 0.204124 | ||||
| 0 | -1 | 0.000000 | 0.288675 | 0.288675 | |||||
| 1 | -1 | 0.000000 | 0.204124 | 0.204124 | |||||
| 1 | -1 | -1 | 0 | 0.000000 | 0.288675 | 0.288675 | |||
| 0 | 0 | 0.408248 | 0.000000 | 0.000000 | |||||
| 1 | 0 | 0.000000 | 0.288675 | 0.288675 | |||||
| 0 | -1 | -1 | 0.408248 | 0.204124 | 0.204124 | ||||
| 0 | -1 | 0.000000 | 0.288675 | 0.288675 | |||||
| 1 | -1 | 0.000000 | 0.204124 | 0.204124 | |||||
| 2 | -1 | -1 | -2 | 2 | 0.000000 | 0.111803 | 0.111803 | ||
| -1 | 2 | 0.000000 | 0.223607 | 0.223607 | |||||
| 0 | 2 | 0.000000 | 0.273861 | 0.273861 | |||||
| 1 | 2 | 0.000000 | 0.223607 | 0.223607 | |||||
| 2 | 2 | 0.447214 | 0.111803 | 0.111803 | |||||
| 0 | -2 | 1 | 0.000000 | 0.158114 | 0.158114 | ||||
| -1 | 1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 1 | 1 | 0.316228 | 0.158114 | 0.158114 | |||||
| 2 | 1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 1 | -2 | 0 | 0.000000 | 0.111803 | 0.111803 | ||||
| 0 | 0 | 0.182574 | 0.091287 | 0.091287 | |||||
| 2 | 0 | 0.000000 | 0.111803 | 0.111803 | |||||
| 0 | -1 | -2 | 1 | 0.000000 | 0.158114 | 0.158114 | |||
| -1 | 1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 1 | 1 | 0.316228 | 0.158114 | 0.158114 | |||||
| 2 | 1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 0 | -2 | 0 | 0.000000 | 0.223607 | 0.223607 | ||||
| 0 | 0 | 0.365148 | 0.182574 | 0.182574 | |||||
| 2 | 0 | 0.000000 | 0.223607 | 0.223607 | |||||
| 1 | -2 | -1 | 0.000000 | 0.158114 | 0.158114 | ||||
| -1 | -1 | 0.316228 | 0.158114 | 0.158114 | |||||
| 1 | -1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 2 | -1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 1 | -1 | -2 | 0 | 0.000000 | 0.111803 | 0.111803 | |||
| 0 | 0 | 0.182574 | 0.091287 | 0.091287 | |||||
| 2 | 0 | 0.000000 | 0.111803 | 0.111803 | |||||
| 0 | -2 | -1 | 0.000000 | 0.158114 | 0.158114 | ||||
| -1 | -1 | 0.316228 | 0.158114 | 0.158114 | |||||
| 1 | -1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 2 | -1 | 0.000000 | 0.158114 | 0.158114 | |||||
| 1 | -2 | -2 | 0.447214 | 0.111803 | 0.111803 | ||||
| -1 | -2 | 0.000000 | 0.223607 | 0.223607 | |||||
| 0 | -2 | 0.000000 | 0.273861 | 0.273861 | |||||
| 1 | -2 | 0.000000 | 0.223607 | 0.223607 | |||||
| 2 | -2 | 0.000000 | 0.111803 | 0.111803 |
[87]:
# Test correlation fns.
# ep.snsMatMod.clustermap(daPlotpd.fillna(0).T.corr()) # This currently failing with NaNs or inf?
ep.snsMatMod.clustermap(daPlotpd.fillna(0).corr())
[87]:
<epsproc._sns_matrixMod.ClusterGrid at 0x221eeaef1d0>
\(\beta_{L,M}^{MF}\) rewrite
The various terms defined above can be used to redefine the full MF observables, expressed as a set of \(\beta_{L,M}\) coefficients.
The original (full) form for the MF equation, as implemented in ePSproc (note this is defined for an expansion in non-conjugate harmonics, \(Y_{LM}(\theta_{\hat{k}},\phi_{\hat{k}})\)):
\begin{eqnarray}
\beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{M}(-1)^{m}(-1)^{(\mu'-\mu_{0})}\left(\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}\right)^{1/2}\left(\begin{array}{ccc}
l & l' & L\\
0 & 0 & 0
\end{array}\right)\left(\begin{array}{ccc}
l & l' & L\\
-m & m' & -M
\end{array}\right)\nonumber \\
& \times & \sum_{P,R',R}(2P+1)(-1)^{(R'-R)}\left(\begin{array}{ccc}
1 & 1 & P\\
\mu & -\mu' & R'
\end{array}\right)\left(\begin{array}{ccc}
1 & 1 & P\\
\mu_{0} & -\mu_{0} & R
\end{array}\right)D_{-R',-R}^{P}(R_{\hat{n}})I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E)
\end{eqnarray}
Where \(I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)\) are the energy-dependent dipole matrix elements.
Substituting the terms above (plus some tidying up with phases etc. - see other notes for full derivation):
\begin{eqnarray}
\beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}B_{L,-M}}\nonumber \\
& \times & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e})\Lambda_{R',R}(R_{\hat{n}})}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E)
\end{eqnarray}
And can also reorder…
\begin{eqnarray*}
\beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e})}\sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}B_{L,-M}\Lambda_{R',R}(R_{\hat{n}})}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E)
\end{eqnarray*}
Add additional (function) labels to track QNs:
\begin{eqnarray*}
\beta_{L,-M}^{\mu_{i},\mu_{f}} & = & \sum_{P,R',R}{\color{red}E_{P-R}(\hat{e};\mu_{0})}\sum_{l,m,\mu}\sum_{l',m',\mu'}(-1)^{(\mu'-\mu_{0})}{\color{red}\Lambda_{R',R}(R_{\hat{n}};\mu,P,R,R')B_{L,-M}(l,l',m,m')}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)I_{l',m',\mu'}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(E)
\end{eqnarray*}
Cf. R&U rotational wavepacket case, see also Chpt. 12 in Quantum Met book for similar formalism separating geometric channels:
\begin{equation}
\varUpsilon_{L,M}^{u}=\sum_{ll'}\sum_{\lambda\lambda'}\sum_{mm'}(-1)^{m}\sqrt{\frac{(2l+1)(2l'+1)(2L+1)}{4\pi}}\left(\begin{array}{ccc}
l & l' & L\\
m & -m' & M
\end{array}\right)\left(\begin{array}{ccc}
l & l' & L\\
0 & 0 & 0
\end{array}\right)\gamma_{\alpha_{+}l\lambda ml'\lambda'm'}^{u}\label{eq:geom-channels}
\end{equation}
\begin{equation}
\Lambda^{\zeta\zeta'}=\sum_{\Gamma,\Gamma'}\sum_{\mu,\mu'}\sum_{h,h'}b_{hl\lambda}^{\Gamma\mu*}b_{h'l'\lambda'}^{\Gamma'\mu'}\boldsymbol{D}_{hl}^{\Gamma\mu*}(q,\,k)\boldsymbol{D}_{h'l'}^{\Gamma'\mu'}(q',\,k)
\end{equation}
\begin{equation}
\beta_{L,M}^{u}=\sum_{\zeta,\zeta'}\varUpsilon_{L,M}^{u,\zeta\zeta'}\Lambda^{\zeta\zeta'}
\end{equation}
Testing multiplications…
[88]:
lmax = 2
p = [0] # Set single pol state only - NOTE EPR() currently sets all ep=1, need to fix this!
# Calculate various tensors...
# *** EPR
EPRX = geomCalc.EPR(form = 'xarray', p = p)
# Set parameters to restack the Xarray into (L,M) pairs
# plotDimsRed = ['l', 'p', 'lp', 'R-p']
plotDimsRed = ['p', 'R-p']
xDim = {'PR':['P','R']}
# Plot with ep.lmPlot(), real values
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(EPRX, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
# *** Blm term to lmax
BLMtable = geomCalc.betaTerm(Lmax = lmax, form = 'xdaLM')
plotDimsRed = ['l', 'm', 'lp', 'mp']
xDim = {'LM':['L','M']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(BLMtable, plotDims=plotDimsRed, xDim=xDim, pType = 'r')
# *** Lambda term
lambdaTerm, lambdaTable, lambdaD, QNs = geomCalc.MFproj(form = 'xarray')
plotDimsRed = ['l', 'lp', 'P', 'mu', 'mup', 'R', 'Rp']
xDim = 'Labels'
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(lambdaTerm, plotDims=plotDimsRed, xDim=xDim, pType = 'a') # Plot abs values
#%% Basic multiplication (no matrix elements)
# Fix/remove QNs for multiplication
# EPRXresort = EPRX.unstack().squeeze().drop('R-p') # Not quite right - need to flatten on R-p?
# EPRXresort = EPRX.unstack().squeeze() # This removes singleton dims from data, but leaves labels.
EPRXresort = EPRX.unstack().squeeze().drop('l').drop('lp') # This removes (l,lp) dims fully.
lambdaTermResort = lambdaTerm.squeeze().drop('l').drop('lp')
BLMtableResort = BLMtable.unstack()
# Apply additional phase conventions
Rphase = True
Mphase = True
mupPhase = True
# Set phase conventions +ve <> -ve terms by in-place coord multiplications
# This should allow for simple selection of values by coords.
if Rphase:
EPRXresort['R'] *= -1
if Mphase:
BLMtableResort['M'] *= -1
# Test big mult...
mTerm = EPRXresort * lambdaTermResort * BLMtableResort # Multiplication works OK... but might be an ugly result...
# Set additional phase term, (-1)^(mup-p)
if mupPhase:
mupPhaseTerm = np.power(-1, np.abs(mTerm.mup - mTerm.p))
mTerm = mTerm * mupPhaseTerm
mTerm.attrs['dataType'] = 'MulTest'
mTerm.attrs['file'] = 'MulTest' # Temporarily adding this, not sure why this is an issue here however (not an issue for other cases...)
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
No handles with labels found to put in legend.
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev\lib\site-packages\xarray\core\nputils.py:215: RuntimeWarning:
All-NaN slice encountered
Plotting data (No filename), pType=r, thres=0.01, with Seaborn
Plotting data (No filename), pType=a, thres=0.01, with Seaborn
[89]:
mTerm.dims
[89]:
('P', 'R-p', 'R', 'Rp', 'mu', 'mup', 'Euler', 'l', 'lp', 'L', 'm', 'mp', 'M')
[90]:
mTerm.shape
[90]:
(3, 3, 3, 5, 3, 3, 3, 3, 3, 5, 5, 5, 9)
[91]:
# Try making sense of this...
# NOTE throwing values to console is very slow here due to high number of dims... calcs and plotting are OK however.
# sumDims = ['l', 'lp', 'P', 'mu', 'mup', 'Rp', 'm', 'mp'] # Define dims to sum over
sumDims = ['P', 'mu', 'mup', 'Rp', ] # Define dims to sum over
# selDims = {'p':0, 'R':0, 'R-p':0} # Use selDims to select pol state
selDims = {'R':0, 'R-p':0} # Use selDims to select pol state
# plotDimsRed = ['L','M']
plotDimsRed = ['L','M', 'l', 'lp', 'm', 'mp']
xDim = 'Labels'
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(mTerm, plotDims=plotDimsRed, xDim=xDim, sumDims=sumDims, selDims=selDims, pType = 'a') # Plot abs values
Plotting data MulTest, pType=a, thres=0.01, with Seaborn
[92]:
# Single pol state and geometry....
sumDims = ['P', 'mu', 'mup', 'Rp', ] # Define dims to sum over
# selDims = {'p':0, 'R':0, 'R-p':0,'Labels':'z'} # Use selDims to select pol state
selDims = {'R':0, 'R-p':0,'Labels':'z'} # Use selDims to select pol state
# plotDimsRed = ['L','M']
plotDimsRed = ['l', 'lp', 'm', 'mp']
xDim = {'LM':['L','M']}
daPlot, daPlotpd, legendList, gFig = ep.lmPlot(mTerm, plotDims=plotDimsRed, xDim=xDim, sumDims=sumDims, selDims=selDims, pType = 'a') # Plot abs values
# This currently includes many non-valid (L,M) - need to further reindex before plotting? They are not present in full plot above.
Plotting data MulTest, pType=a, thres=0.01, with Seaborn
[93]:
# daPlotpd # .dropna(how='any',axis=1) - many invalid null terms here??? (Hence can't drop with dropna)
daPlotpd[daPlotpd.abs() >= 0.01].dropna(how='all', axis=1) # Slightly better!
[93]:
| L | 0 | 1 | 2 | 3 | 4 | |||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| M | 0 | 1 | 0 | -1 | 2 | 1 | 0 | -1 | -2 | 3 | ... | -3 | 4 | 3 | 2 | 1 | 0 | -1 | -2 | -3 | -4 | |||
| l | lp | m | mp | |||||||||||||||||||||
| 0 | 0 | 0 | 0 | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 0 | -1 | NaN | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| 0 | NaN | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | NaN | 0.282095 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 2 | 0 | -2 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.282095 | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2 | 2 | 2 | -2 | 0.282095 | NaN | NaN | NaN | NaN | NaN | 0.180224 | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | 0.040299 | NaN | NaN | NaN | NaN |
| -1 | NaN | NaN | NaN | NaN | NaN | 0.220728 | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | 0.090112 | NaN | NaN | NaN | NaN | NaN | |||
| 0 | NaN | NaN | NaN | NaN | 0.180224 | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | 0.156078 | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | NaN | 0.238414 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
| 2 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | NaN | 0.337168 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |||
81 rows × 25 columns
[94]:
# With cluster sorting
ep.snsMatMod.clustermap(daPlotpd[daPlotpd.abs() >= 0.01].dropna(how='all', axis=1).fillna(0))
[94]:
<epsproc._sns_matrixMod.ClusterGrid at 0x2218048de48>
Calculations with matrix elements to be explored in the next installment…