Maximum beta is not Emax/kT when exporting S(a,b) from Info object

[marquezj]

For future reference after discussion with @tkittel.

(This one qualifies as "this is a feature, not a bug").

When exporting S(a,b) to use in other applications, for instance exporting to ENDF-6 files, the beta grid does not extend to completely cover the range of energies from 0 to Emax. The rule of thumb to calculate the range of alpha and beta, given by Mattes and Keinert is:

$\beta_{max} = E_{max} /kT$
$\alpha_{max} = 4 \beta_{max}$

(here we use the definition of $\alpha$ used in NCrystal)

Now, this is an upper bound for the beta range. But, depending on the mass of the scatterer, that beta range can be outside of the energy exchanges that can be possible from the kinematics. The whole range will only be accessible for $^{1}$H.

In NCrystal alpha and beta grid generation is adaptive, as it can be seen here, using a fix number of points given by vdoslux. The methodology is complex, but by checking the shape of the contribution from the different phonons, the algorithm ensures that the points are used where there are needed.

import NCrystal as NC
import numpy as np
import matplotlib.pyplot as plt

cfg = 'Al_sg225.ncmat;temp=293.6K'

mat_info = NC.createInfo(cfg)
mat_di = mat_info.dyninfos[0]
mat_knl = mat_di.loadKernel(vdoslux=3)

alpha = mat_knl['alpha']
beta = mat_knl['beta']
sab = mat_knl['sab']
Emax = mat_knl['suggestedEmax']
kT = mat_di.temperature*NC.constant_boltzmann

XX,YY = np.meshgrid(alpha, beta)
sab.shape = np.shape(XX)
plt.pcolormesh(XX,YY,sab, cmap = plt.colormaps['jet'])
plt.xlabel('alpha')
plt.ylabel('beta')
plt.colorbar()

mattes_betamax = mat_knl['suggestedEmax']/kT
mattes_alphamax = 4*mattes_betamax # NCrystal does not divide by the atomic weight ratio

print(f'Emax: {Emax:.2f}')
print(f'Alpha range: {min(alpha):.2f}, {max(alpha):.2f}')
print(f'Beta range: {min(beta):.2f}, {max(beta):.2f}')

print(f'Mattes alpha max: 4Emax/kT = {mattes_alphamax:.2f}')
print(f'Mattes beta max: Emax/kT = {mattes_betamax:.2f}')

b = np.linspace(-mattes_betamax, max(beta), 1000)

alpha_max = (2*Emax+kT*b+2*np.sqrt(Emax*(Emax+kT*b)))/kT
alpha_min = (2*Emax+kT*b-2*np.sqrt(Emax*(Emax+kT*b)))/kT

plt.plot(alpha_max, b, 'r-')
plt.plot(alpha_min, b, 'r-')

> Emax: 5.00 eV
> Alpha range: 0.00, 839.70
> Beta range: -94.74, 20.00
> Mattes alpha max: 4Emax/kT = 790.50
> Mattes beta max: Emax/kT = 197.62