lynntf / nist-calculators

Python implementaion of some NIST calculators
MIT License
1 stars 3 forks source link

Attenuation interpolator transition at photoelectric edges #2

Open lynntf opened 1 year ago

lynntf commented 1 year ago

Oh, another thing I want to say: the attenuation interpolator for the photoelectric edges is broken because it fails to adequately do the handoff between the log-log linear interpolation and the cubic spline interpolation regimes. This means that specifying a better energy resolution for interpolation (like instead of the default energy grid, something like np.logspace(3,11, 20000)) causes the code to break because it's considered extrapolation and scipy throws an error. Altering interpolators.py such that


def _interpolateAbsorptionEdge(data) -> Callable[[np.ndarray], np.ndarray]:
data, h5file, group = data
    data_K = h5file.get_node(group, "K").read()
    cubicSplineThreshold = np.max(data_K[NameProcess.ENERGY]) * 1e6
    x = data[NameProcess.ENERGY]
    y = data[NameProcess.PHOTOELECTRIC]
    indx = x > cubicSplineThreshold
    cs = CubicSpline(np.log(x[indx]), np.log(y[indx]), bc_type='natural')
    linear = interp1d(np.log(x[~indx]), np.log(y[~indx]), kind='linear')
> instead has inserted the line
```python
    def _interpolateAbsorptionEdge(data) -> Callable[[np.ndarray], np.ndarray]:
        data, h5file, group = data

        data_K = h5file.get_node(group, "K").read()
        cubicSplineThreshold = np.max(data_K[NameProcess.ENERGY]) * 1e6
        x = data[NameProcess.ENERGY]
        y = data[NameProcess.PHOTOELECTRIC]
        indx = x > cubicSplineThreshold
        indx = np.append(False,indx[:-1])     # <----- this line right here!
        cs = CubicSpline(np.log(x[indx]), np.log(y[indx]), bc_type='natural')
        linear = interp1d(np.log(x[~indx]), np.log(y[~indx]), kind='linear')

alleviates this problem, and now one can add as many points on the energy regime as they like. The issue is that I'm not sure this is entirely kosher... so I'm not uploading it as a pull request right now, but I wanted somewhere to write this down. Thanks!

Originally posted by @ischreiber3 in https://github.com/lynntf/nist-calculators/issues/1#issuecomment-1745599453

lynntf commented 1 year ago

@ischreiber3 I'm not seeing this behavior. Can you provide a minimal example to demonstrate it?

One thing I do see is that the spline interpolation for coherent and incoherent scattering blow up once the energy gets too small/too large. I've added linear extrapolation to all of the interpolators to address the blow-up and added warnings when data outside of the tabulated data from XCOM is requested in commit ac82696