Interpolated Mass attenuation values don't seem to match the reference table data

[samaloney]

Comparing the graphs from the NIST source page here the edges don't seem to be being interpolated properly.

Code to create plot

e = np.linspace(1e-3, 20, 10000)*u.MeV
te = MassAttenuationCoefficient('Te')
plt.plot(e, te.func(e))
plt.loglog()
plt.xlim(1e-3, 1e2)
plt.ylim(1e-2, 1e4)
plt.xlabel('Photon Energy (MeV)')
plt.ylabel(r'$\mu / \rho$ (cm$^{2}$ g$^{-1}$)')
plt.title('Te')

[samaloney]

It seems scipy's interp1d doesn't handle multivalued/non-monotonic x-values, attached plot shows the effect, evaluating the attenuation coefficient at same values used for the interpolation.

[rschwartz70]

where is the interp1d done in your codelet. If I could see the issue in detail I could determine whether my IDL code would suffer the same issue.

[samaloney]

It's these lines, where self.energy.value and self.data.value are values from the NIST tables and the interpolation function is setup as self._f and then calling self._f(<log of engery>) returns the an interpolated value.

roentgen/roentgen/absorption/material.py

Lines 312 to 320 in 4446758

	data_energy_kev = np.log10(self.energy.value)
	data_attenuation_coeff = np.log10(self.data.value)
	self._f = interpolate.interp1d(
	data_energy_kev, data_attenuation_coeff, bounds_error=False,
	fill_value=0.0
	)
	self.func = lambda x: u.Quantity(
	10 ** self._f(np.log10(x.to("keV").value)), "cm^2/g"
	)

I altered the first of any repeating data points by a tiny fraction (1e-10) and this resolved the issue but not sure if that is a good solution or not.

[rschwartz70]

interp1d should be called using assume_sorted as true because the data are already sorted

[rschwartz70]

More extensive NIST tables exist into the pair production energy range. There are several independent processes that go into the mass atten coefficients that should be independently specifiable. It would require a trivial elaboration of this code. I am going to enlist Laura and Danny on this effort.

[ehsteve]

Hey @rschwartz70 and @samaloney! sorry it took me so long to reply, I was on vacation last week. Thanks for noticing and reporting this bug. I had thought I had looked at this but apparently not closely enough (see https://roentgen.readthedocs.io/en/latest/examples/tutorial1.html)! It looks like @rschwartz70 solution works for me. I'll admit that I don't understand why it works because the documentation (https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html) states that it will sort the data itself by default so this should not be needed.

[ehsteve]

@rschwartz70 I'd be happy to expand the data set here if you can point me to the right place. Can I ask what you guys are working on? I'm happy to see that some is beginning to use this!

[rschwartz70]

This is for STIX but I'm going to use these NIST data for Be to fix the GOES16 response problem. But it also demonstrated that using XCOM data from NIST is a little better than the Berger and Seltzer tables/fits that I've been using in XSEC (ssw/packages/xray) for 40 years (pre IDL). The XCOM tables from NIST used date back to before 1995. Those tables can be found in a FITS file in SSW $SSW/packages/xray/dbase/xcom.fits read_mat_xcom.pro in xray is an idl interface into individual tables. An updated XCOM version, with ASCII tables, is found here https://physics.nist.gov/PhysRefData/Xcom/Text/download.html The DOS file is an archive that can be unpacked by 7zip. Probably the MAC version works for the rest of you. I may use the DOS files to update my tables. A word about xsec. Those were log/log stepwise fits around edges. Put out by NIST people. Probably the biggest issue is that in the XCOM version there is much more edge detail and there are no log/log coefficients, only the data and the double energy double cross-section solution at the k/l edges which caused the scipy and not the numpy interpol1d glitch I'd be so happy if we had a unified database between ssw and StevePy for xsections. Thanks

…

On Mon, Aug 31, 2020 at 9:37 AM Steven Christe ***@***.***> wrote: @rschwartz70 <https://github.com/rschwartz70> I'd be happy to expand the data set here if you can point me to the right place. Can I ask what you guys are working on? I'm happy to see that some is beginning to use this! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#24 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADMUMXYAVSSAAK7VHXOXASTSDORTDANCNFSM4QNHD3BQ> .

[ehsteve]

@samaloney would you be willing to review my proposed fix in PR #25?

[samaloney]

Just a final note on this the assume_sorted keyword only works for floats which is fine in this case as astropysQuantity are always cast to floats.

[rschwartz70]

The xcom FITS file I've been using has a typo in the spelling of Fluorine. It has used Flourine. That's not US or British English. I'm going to correct that entry and load the file back into ssw/packages/xray/dbase/xcom.fits Richard On Mon, Aug 31, 2020 at 10:15 AM Richard Schwartz <rschwartz70@gmail.com> wrote:

…

This is for STIX but I'm going to use these NIST data for Be to fix the GOES16 response problem. But it also demonstrated that using XCOM data from NIST is a little better than the Berger and Seltzer tables/fits that I've been using in XSEC (ssw/packages/xray) for 40 years (pre IDL). The XCOM tables from NIST used date back to before 1995. Those tables can be found in a FITS file in SSW $SSW/packages/xray/dbase/xcom.fits read_mat_xcom.pro in xray is an idl interface into individual tables. An updated XCOM version, with ASCII tables, is found here https://physics.nist.gov/PhysRefData/Xcom/Text/download.html The DOS file is an archive that can be unpacked by 7zip. Probably the MAC version works for the rest of you. I may use the DOS files to update my tables. A word about xsec. Those were log/log stepwise fits around edges. Put out by NIST people. Probably the biggest issue is that in the XCOM version there is much more edge detail and there are no log/log coefficients, only the data and the double energy double cross-section solution at the k/l edges which caused the scipy and not the numpy interpol1d glitch I'd be so happy if we had a unified database between ssw and StevePy for xsections. Thanks On Mon, Aug 31, 2020 at 9:37 AM Steven Christe ***@***.***> wrote: > @rschwartz70 <https://github.com/rschwartz70> I'd be happy to expand the > data set here if you can point me to the right place. Can I ask what you > guys are working on? I'm happy to see that some is beginning to use this! > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#24 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ADMUMXYAVSSAAK7VHXOXASTSDORTDANCNFSM4QNHD3BQ> > . >

[ehsteve]

This was fixed.