In [1]:
import sys
import os
import time
import numpy as np
# For module testing, include path to module here
modPath = r'D:\code\github\ePSproc'
sys.path.append(modPath)
import epsproc as ep
In [2]:
# Load data from modPath\data
dataPath = os.path.join(modPath, 'data')
# Scan data dir
dataSet = ep.readMatEle(fileBase = dataPath)
The $\beta_{L,M}$ parameters are defined as:
$ \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} $
Calculations use ep.mfblm(), which will calculate all values at each energy point for the supplied dataset. This may take a while in some cases due to multiple nested sums - this code will be parallelised in future.
Calculate $\beta_{LM}$ as function of energy.
In [3]:
daIn = dataSet[0].copy()
# BLMXeN2 = ep.mfblm(daIn[:, 1:4], selDims = {'Type':'L'}, thres = 1e-4) # Subselected on Eke
start = time.time()
BLMXeN2 = ep.mfblm(daIn, selDims = {'Type':'L'}, thres = 1e-4) # Run for all Eke
end = time.time()
print('Elapsed time = {0} seconds, for {1} energy points.'.format((end-start), BLMXeN2.Eke.size))
In [4]:
BLMXeN2
Out[4]:
In [5]:
# Plot using Xarray functionality with thresholding
BLMXeN2.where(np.abs(BLMXeN2) > 1e-4, drop = True).real.squeeze().plot.line(x='Eke')
Out[5]:
In [6]:
# Plot values normalised by B00
# This seems to work... probably a more elegant solution here, since this assumes dimension order.
normBLM = BLMXeN2/BLMXeN2[:,:,0]
normBLM.where(np.abs(normBLM) > 1e-4, drop = True).real.squeeze().plot.line(x='Eke')
Out[6]:
In [7]:
# Calculate & plot MFPADs from BLMs
# def MFPAD_BLM(BLMXin):
# # Calculate YLMs
# YLMX = ep.sphCalc(BLMXin.l.max(), res=50)
# YLMX = YLMX.rename({'LM':'BLM'}) # Switch naming for multiplication & plotting
# MFPAD = BLMXin*YLMX
# MFPAD = MFPAD.rename({'BLM':'LM'})
# return MFPAD
# MFPAD = MFPAD_BLM(BLMXeN2)
MFPAD_N2, _ = ep.sphFromBLMPlot(BLMXeN2)
In [8]:
# Plot single E with matplotlib
ep.sphSumPlotX(MFPAD_N2.sel({'Eke':1.1}).squeeze(), pType = 'r', backend = 'mpl')
Out[8]:
In [9]:
# Plot MFPAD surfaces vs E
print('N2 test data, MFPADs vs E')
MFPAD_N2.sum('LM').squeeze().isel(Eke=slice(0,10,2)).real.plot(x='Theta',y='Phi', col='Eke')
Out[9]:
In [10]:
# Plot multiple E with matplotlib
ep.sphSumPlotX(MFPAD_N2.isel(Eke=slice(0,50,10)), pType = 'r', backend = 'mpl')
Out[10]:
In [11]:
# Plot multiple E with plotly (work in progress!)
# Broken since adding Euler angs...?
_ = ep.sphSumPlotX(MFPAD_N2.squeeze().isel(Eke=slice(0,50,10)), pType = 'r', backend = 'pl')
Benchmark results, single energy, multiple polarization geometries.
In [12]:
daIn = dataSet[1].copy()
# Set threshold for matrix elements & BLMs - this can speed up calculations significantly, but will affect accuracy.
thres = 1e-4
# For all pol geoms
pRot = [0, 0, np.pi/2]
tRot = [0, np.pi/2, np.pi/2]
cRot = [0, 0, 0]
eAngs = np.array([pRot, tRot, cRot]).T # List form to use later, rows per set of angles
ts = []
BLMXeNO2list = []
for eIn in range(0,3):
start = time.time()
BLMXeNO2list.append(ep.mfblm(daIn, selDims = {'Type':'L'}, eAngs = eAngs[eIn,:], thres=thres))
ts.append(time.time()-start)
print('Elapsed time = {0} seconds'.format(ts[-1]))
In [13]:
# Stack results - should improve on this and add eAngs stacking to mfblm code, but OK for testing
import xarray as xr
BLMXeNO2 = xr.combine_nested(BLMXeNO2list,'Euler')
# xr.combine_nested(Tlm, concat_dim=['Euler'])
# Plot values normalised by B00
# This seems to work... probably a more elegant solution here, since this assumes dimensions.
normBLM = BLMXeNO2/BLMXeNO2[:,:,0]
# With mag checking to avoid spurious division by small B00 terms...
# normBLM = BLMXeNO2.where(np.abs(BLMXeNO2[:,:,0]) > 1e-4, drop = True)
# normBLM = normBLM/normBLM[:,:,0]
# Replace multi-index with linear index for plotting (otherwise get coord errors)
normBLM['Euler'] = np.arange(0,normBLM.Euler.size)
normBLM.where(np.abs(normBLM) > 1e-2, drop = True).squeeze().real.plot.line('-x',x='Euler')
Out[13]:
In [14]:
MFPAD_NO2, _ = ep.sphFromBLMPlot(BLMXeNO2)
# Plot multiple E with matplotlib
ep.sphSumPlotX(MFPAD_NO2, pType = 'r', backend = 'mpl', facetDim = 'Euler')
Out[14]:
Compare against results from ep.mfpad, see main "ePSproc demo notebook" for calculation details.
In [15]:
print('MFPADs for test NO2 dataset (single energy, (z,x,y) pol states)')
TX, TlmX = ep.mfpad(dataSet[1])
# Plot for each pol geom (symmetry)
for n in range(0,3):
ep.sphSumPlotX(TX[n].sum('Sym').squeeze()**2, pType = 'a')
(Visual comparison looks OK, still some benchmarks to to for numerical re/im comparisons, which currently show some differences.)