Merge pull request #36 from aidancrilly/pyendf

Swapped to using ENDF files for cross section data

README.md

@@ -40,25 +40,18 @@ pip install -e .
 ## Current model specifications:
 - Primary spectrum models for DT, DD and TT
 - Elastic and inelastic (n,2n) processes for D and T
+- ENDF file interface using [ENDF-python](https://github.com/paulromano/endf-python)
 - Relativistic corrections to elastic scattering kernels
 - Scattering of all primary neutron sources
 - Inclusion of areal density asymmetry effects and variable fuel fractions
 - Backscatter edge shape effects from scattering ion kinematics
-- Preliminary scattering in H, 12C, 9Be support
 - Synthetic neutron time-of-flight tools
 
 ## Future model developments:
 - Fitting models with ion kinematic approximations
 - Pre-computed table support for backscatter edge matrix
 - Knock-on ion spectra
 
-## DT cross section and reactivity data:
-ENDF: nT elastic, nD elastic, nH elastic, T(n,2n), TT primary
-
-CENDL: D(n,2n), n12C elastic and 4.4 MeV inelastic
-
-Bosch-Hale: DT primary, DD primary
-
 ## Publications:
 The models used in this code are described in the following publications:
 

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "NeSST"
-version = "1.0.0"
+version = "1.1.0"
 description = "Neutron Scattered Spectra Tool, ICF primary and scattered neutron spectroscopy analysis code"
 readme = "README.md"
 classifiers = [
@@ -21,7 +21,8 @@ keywords = ["python", "reactor", "fusion", "power", "sankey"]
 dependencies = [
     "numpy",
     "scipy",
-    "matplotlib"
+    "matplotlib",
+    "endf"
 ]
 
 [project.urls]
@@ -35,7 +36,8 @@ test = [
 ]
 
 [tool.setuptools.package-data]
-"NeSST.data" = ["*.dat", "*.txt","*.csv"]
+"NeSST.data" = ["*.dat", "*.txt","*.csv","*.json"]
+"NeSST.data.ENDF" = ["*.endf", "*.cendl","*.brond"]
 
 [tool.setuptools]
 package-dir = {"" = "src"}

src/NeSST/constants.py

@@ -1,6 +1,8 @@
 ''' Contains some physical constants used throughout the analysis '''
 import scipy.constants as sc
 import numpy as np
+import numpy.typing as npt
+
 # Scipy constants uses CODATA2018 database
 
 # Need to swap from MeV to eV
@@ -27,4 +29,5 @@
 # Directories
 import os
 package_directory = os.path.dirname(os.path.abspath(__file__))
-data_dir = os.path.join(package_directory,"./data/")
+data_dir = os.path.join(package_directory,"./data/")
+ENDF_dir = os.path.join(data_dir,"./ENDF/")

src/NeSST/core.py

@@ -443,10 +443,10 @@ def mat_scatter_spec(mat : typing.Type[sm.material_data],
     
     """
     mat.calc_dNdEs(I_E,rhoL_func)
-    total = mat.elastic_dNdE
+    total = mat.elastic_dNdE.copy()
     if(mat.l_n2n):
         total += mat.n2n_dNdE
-    if(mat.inelastic):
+    if(mat.l_inelastic):
         total += mat.inelastic_dNdE
     return total
 

src/NeSST/cross_sections.py

@@ -1,25 +1,54 @@
 # Backend of spectral model
 import numpy as np
 from numpy.polynomial.legendre import legval
+from dataclasses import dataclass
 from NeSST.constants import *
 from NeSST.utils import *
+from scipy.interpolate import griddata
 import NeSST.collisions as col
 
 ###############################
 # Differential cross sections #
 ###############################
 
+@dataclass
+class NeSST_SDX:
+    Ein    : npt.NDArray
+    points : npt.NDArray
+    values : npt.NDArray
+
+@dataclass
+class NeSST_DDX:
+    NEin  : int
+    Ein   : list
+    Ncos  : list
+    cos   : list
+    NEout : dict
+    Eout  : dict
+    f     : dict
+    Emax  : dict
+
 # Elastic single differential cross sections
 
+def diffxsec_table_eval(sig,mu,E,table):
+
+    xi = np.column_stack((E.flatten(),mu.flatten()))
+    # Rescale is very important, gives poor results otherwise
+    interp = griddata(table.points,table.values,xi,rescale=True).reshape(mu.shape)
+
+    ans = sig*interp
+    ans[np.abs(mu) > 1.0] = 0.0
+    return ans
+
 # Interpolate the legendre coefficients (a_l) of the differential cross section
 # See https://t2.lanl.gov/nis/endf/intro20.html
-def interp_Tlcoeff(dx_spline,E_vec):
+def interp_Tlcoeff(legendre_dx_spline,E_vec):
     size = [E_vec.shape[0]]
-    NTl = len(dx_spline)
+    NTl = len(legendre_dx_spline)
     size.append(NTl)
     Tlcoeff = np.zeros(size)
     for i in range(NTl):
-        Tlcoeff[:,i] = dx_spline[i](E_vec)
+        Tlcoeff[:,i] = legendre_dx_spline[i](E_vec)
     return Tlcoeff,NTl
 
 # Evaluate the differential cross section by combining legendre and cross section
@@ -39,12 +68,12 @@ def diffxsec_legendre_eval(sig,mu,coeff):
 # CoM frame differential cross section wrapper fucntion
 def f_dsdO(Ein_vec,mu,material):
 
-    dx_spline = material.dx_spline
-    Tlcoeff_interp,Nl = interp_Tlcoeff(dx_spline,Ein_vec)
-    Tlcoeff_interp    = 0.5*(2*np.arange(0,Nl)+1)*Tlcoeff_interp
-
     sig = material.sigma(Ein_vec)
 
+    legendre_dx_spline = material.legendre_dx_spline
+    Tlcoeff_interp,Nl = interp_Tlcoeff(legendre_dx_spline,Ein_vec)
+    Tlcoeff_interp    = 0.5*(2*np.arange(0,Nl)+1)*Tlcoeff_interp
+
     dsdO = diffxsec_legendre_eval(sig,mu,Tlcoeff_interp)
     return dsdO
 
@@ -54,151 +83,26 @@ def dsigdOmega(A,Ein,Eout,Ein_vec,muin,muout,vf,material):
     return f_dsdO(Ein_vec,mu_CoM,material)
 
 # Inelastic double differential cross sections
-
-def ENDF_format(x):
-    x_after_decimal = x.split('.')[1]
-    if('+' in x_after_decimal):
-        sign = '+'
-    elif('-' in x_after_decimal):
-        sign = '-'
-    exp  = x.split(sign)[-1]
-    num  = x[:-(len(exp)+1)]
-    if(sign == '-'):
-        return float(num)*10**(-int(exp))
-    elif(sign == '+'):
-        return float(num)*10**(+int(exp))
-    else:
-        print("Strange ENDF float detected....")
-        print(x)
-        return 0.0
-
 # Reads and interpolated data saved in the ENDF interpreted data format
 class doubledifferentialcrosssection_data:
 
-    def __init__(self,filexsec,fileddx,ENDF):
-        self.filexsec = filexsec
-        self.fileddx = fileddx
-        self.ENDF = ENDF
-        if(ENDF):
-            self.read_xsec_file()
-            self.read_ddx_file()
-        else:
-            self.read_ddx_file_csv()
-
-    def read_xsec_file(self):
-        with open(self.filexsec,"r") as f:
-            file = f.read()
-            # Read number of points
-            self.NEin_xsec = int(file.split()[0])
-            self.Ein_xsec  = np.zeros(self.NEin_xsec)
-            self.xsec      = np.zeros(self.NEin_xsec)
-            counter        = 0
-            data = "".join(file.split("\n")[5:]).split()
-            E = data[::2]
-            x = data[1::2]
-            for i in range(self.NEin_xsec):
-                self.Ein_xsec[i] = ENDF_format(E[i])
-                self.xsec[i]     = ENDF_format(x[i])
-        self.xsec_interp = interpolate_1d(self.Ein_xsec,self.xsec,method='linear',bounds_error=False,fill_value=0.0)
-
-    def read_ddx_file(self):
-        with open(self.fileddx,"r") as f:
-            file = f.read().split("\n")
-            # Read number of points
-            self.NEin_ddx = int(file[1].split()[0])
-            self.Ein_ddx  = np.zeros(self.NEin_ddx)
-            # Read data
-            Ecounter = 0
-            Ccounter = 0
-            read_cosine = False
-            read_energy = False
-            read_data   = False
-            self.Ncos_ddx  = []
-            self.cos_ddx   = []
-            self.NEout_ddx = {}
-            self.Eout_ddx  = {}
-            self.f_ddx     = {}
-            self.f_ddx_interp = {}
-            self.Emax_ddx     = {}
-            # Read in all the array axes
-            for Lcounter,line in enumerate(file):
-                line_split = line.split()
-                if(line_split != []):
-                    # Read number of cosines for given incoming energy
-                    if(read_cosine):
-                        self.Ncos_ddx.append(int(line_split[0]))
-                        self.cos_ddx.append(np.zeros(int(line_split[0])))
-                        read_cosine = False
-                        Ccounter = 0
-                    # Read number of energies for given incoming energy and cosine
-                    if(read_energy):
-                        NEout = int(line_split[0])
-                        self.NEout_ddx[(Ecounter-1,Ccounter-1)] = NEout
-                        self.Eout_ddx[(Ecounter-1,Ccounter-1)]  = np.zeros(NEout)
-                        self.f_ddx[(Ecounter-1,Ccounter-1)]     = np.zeros(NEout)
-                        read_energy = False
-                        idx1 = Lcounter + 4
-                        idx2 = idx1 + int(np.ceil(NEout/3)) + 1
-                        read_data   = True
-                    # Read in the data
-                    if(read_data):
-                        data = "".join(file[idx1:idx2]).split()
-                        E = data[::2]
-                        x = data[1::2]
-                        for i in range(NEout):
-                            self.Eout_ddx[(Ecounter-1,Ccounter-1)][i] = ENDF_format(E[i])
-                            self.f_ddx[(Ecounter-1,Ccounter-1)][i]    = ENDF_format(x[i])
-                        self.f_ddx_interp[(Ecounter-1,Ccounter-1)] = interpolate_1d(self.Eout_ddx[(Ecounter-1,Ccounter-1)],self.f_ddx[(Ecounter-1,Ccounter-1)],method='linear',bounds_error=False,fill_value=0.0)
-                        self.Emax_ddx[(Ecounter-1,Ccounter-1)] = np.max(self.Eout_ddx[(Ecounter-1,Ccounter-1)])
-                        read_data = False
-                    # Read incoming energy
-                    if(line_split[0] == 'Energy:'):
-                        self.Ein_ddx[Ecounter] = ENDF_format(line_split[1])
-                        Ecounter += 1
-                    # Prep for cosine number read in
-                    elif(" ".join(line_split) == 'Cosine Interpolation:'):
-                        read_cosine = True
-                    # Read number of secondary energies
-                    elif(" ".join(line_split) == 'Secondary-Energy Interpolation:'):
-                        read_energy = True
-                    elif('Cosine:' in line):
-                        line_split_c = line.split(":")[1]
-                        self.cos_ddx[Ecounter-1][Ccounter] = ENDF_format(line_split_c)
-                        Ccounter += 1
-
-    def read_ddx_file_csv(self):
-        self.NEin_ddx = 2
-        self.Ein_ddx  = np.array([0.0,14.0e6])
-        data = np.loadtxt(self.fileddx,delimiter=',',skiprows=1)
-        angles,counts = np.unique(data[:,-1],return_counts=True)
-        cos = np.cos(angles[::-1]*np.pi/180.)
-        NC = cos.shape[0]
-        self.Ncos_ddx  = [NC,NC]
-        self.cos_ddx   = [cos,cos]
-        E_prev = 0.0
-        self.NEout_ddx = {}
-        self.Eout_ddx  = {}
-        self.f_ddx     = {}
+    def __init__(self,ENDF_LAW6_xsec_data,ENDF_LAW6_dxsec_data):
+        self.xsec_interp = interpolate_1d(ENDF_LAW6_xsec_data['E'],ENDF_LAW6_xsec_data['sig'],method='linear',bounds_error=False,fill_value=0.0)
+
+        self.NEin_ddx  = ENDF_LAW6_dxsec_data['DDX'].NEin
+        self.Ein_ddx   = ENDF_LAW6_dxsec_data['DDX'].Ein
+        self.Ncos_ddx  = ENDF_LAW6_dxsec_data['DDX'].Ncos
+        self.cos_ddx   = ENDF_LAW6_dxsec_data['DDX'].cos
+        self.NEout_ddx = ENDF_LAW6_dxsec_data['DDX'].NEout
+        self.Eout_ddx  = ENDF_LAW6_dxsec_data['DDX'].Eout
+        self.f_ddx     = ENDF_LAW6_dxsec_data['DDX'].f
+        self.Emax_ddx  = ENDF_LAW6_dxsec_data['DDX'].Emax
+
+        # Build interpolator dict
         self.f_ddx_interp = {}
-        self.Emax_ddx     = {}
-        idx = data[:,0].shape[0]
-        i = 0
-        for ic in range(NC-1,-1,-1):
-            NEout = counts[ic]
-            self.NEout_ddx[(0,i)] = NEout
-            self.Eout_ddx[(0,i)]  = data[idx-NEout:idx,1]
-            self.f_ddx[(0,i)]     = np.zeros(NEout)
-            self.NEout_ddx[(1,i)] = NEout
-            self.Eout_ddx[(1,i)]  = data[idx-NEout:idx,1]
-            # From barns to mbarns, from sr to per cosine, from number of neutrons to cross section
-            self.f_ddx[(1,i)]     = 0.5*(2*np.pi)*data[idx-NEout:idx,0]/1e3
-            self.f_ddx_interp[(0,i)] = interpolate_1d(self.Eout_ddx[(0,i)],self.f_ddx[(0,i)],method='linear',bounds_error=False,fill_value=0.0)
-            self.Emax_ddx[(0,i)] = np.max(self.Eout_ddx[(0,i)])
-            self.f_ddx_interp[(1,i)] = interpolate_1d(self.Eout_ddx[(1,i)],self.f_ddx[(1,i)],method='linear',bounds_error=False,fill_value=0.0)
-            self.Emax_ddx[(1,i)] = np.max(self.Eout_ddx[(1,i)])
-            idx -= NEout
-            i   += 1
-        self.xsec_interp = lambda x : 1.
+        for Ecounter in range(self.NEin_ddx):
+            for Ccounter in range(self.Ncos_ddx[Ecounter]):
+                self.f_ddx_interp[(Ecounter,Ccounter)] = interpolate_1d(self.Eout_ddx[(Ecounter,Ccounter)],self.f_ddx[(Ecounter,Ccounter)],method='linear',bounds_error=False,fill_value=0.0)
 
     # Interpolate using Unit Base Transform
     def interpolate(self,Ein,mu,Eout):
@@ -285,31 +189,14 @@ def regular_grid(self,Ein,mu,Eout):
 
 class doubledifferentialcrosssection_LAW6:
 
-    def __init__(self,filexsec,A_i,A_e,A_t,A_p,A_tot,Q_react):
-        self.A_i     = A_i
-        self.A_e     = A_e
-        self.A_t     = A_t
-        self.A_p     = A_p
-        self.A_tot   = A_tot
-        self.Q_react = Q_react
-        self.filexsec = filexsec
-        self.read_xsec_file()
-
-    def read_xsec_file(self):
-        with open(self.filexsec,"r") as f:
-            file = f.read()
-            # Read number of points
-            self.NEin_xsec = int(file.split()[0])
-            self.Ein_xsec  = np.zeros(self.NEin_xsec)
-            self.xsec      = np.zeros(self.NEin_xsec)
-            counter        = 0
-            data = "".join(file.split("\n")[5:]).split()
-            E = data[::2]
-            x = data[1::2]
-            for i in range(self.NEin_xsec):
-                self.Ein_xsec[i] = ENDF_format(E[i])
-                self.xsec[i]     = ENDF_format(x[i])
-        self.xsec_interp = interpolate_1d(self.Ein_xsec,self.xsec,method='linear',bounds_error=False,fill_value=0.0)
+    def __init__(self,ENDF_LAW6_xsec_data,ENDF_LAW6_dxsec_data):
+        self.A_i     = ENDF_LAW6_dxsec_data['A_i']
+        self.A_e     = ENDF_LAW6_dxsec_data['A_e']
+        self.A_t     = ENDF_LAW6_dxsec_data['A_t']
+        self.A_p     = ENDF_LAW6_dxsec_data['A_p']
+        self.A_tot   = ENDF_LAW6_dxsec_data['A_tot']
+        self.Q_react = ENDF_LAW6_dxsec_data['Q_react']
+        self.xsec_interp = interpolate_1d(ENDF_LAW6_xsec_data['E'],ENDF_LAW6_xsec_data['sig'],method='linear',bounds_error=False,fill_value=0.0)
 
     def ddx(self,Ein,mu,Eout):
         E_star = Ein*self.A_i*self.A_e/(self.A_t+self.A_i)**2

src/NeSST/data/Be9.json

@@ -0,0 +1,9 @@
+{
+    "n-004_Be_009.endf" :
+    {
+        "total" : true,
+        "elastic" : true,
+        "inelastic" : true,
+        "n2n" : true
+    }
+}

src/NeSST/data/C12.json

@@ -0,0 +1,9 @@
+{
+    "n-006_C_012.endf" :
+    {
+        "total" : true,
+        "elastic" : true,
+        "inelastic" : true,
+        "n2n" : false
+    }
+}