From 0a4b4ec77093019a14d6e8ff8a24ff8cb21c7234 Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Thu, 24 Jul 2025 12:23:03 +0100 Subject: [PATCH 1/5] Performed ruff formatting --- .github/workflows/python-publish.yml | 4 +- .gitignore | 5 +- .pre-commit-config.yaml | 50 ++++ README.md | 2 +- pyproject.toml | 2 +- ruff.toml | 49 ++++ src/NeSST/__init__.py | 2 - src/NeSST/collisions.py | 313 +++++++++++--------- src/NeSST/constants.py | 50 ++-- src/NeSST/core.py | 419 +++++++++++++++------------ src/NeSST/cross_sections.py | 235 ++++++++------- src/NeSST/data/Be9.json | 2 +- src/NeSST/data/C12.json | 2 +- src/NeSST/data/H1.json | 2 +- src/NeSST/data/H2.json | 4 +- src/NeSST/data/H3.json | 4 +- src/NeSST/data/TT_reac_Hale.dat | 2 +- src/NeSST/endf_interface.py | 310 +++++++++++--------- src/NeSST/fitting.py | 356 +++++++++++++---------- src/NeSST/spectral_model.py | 367 +++++++++++++---------- src/NeSST/time_of_flight.py | 343 +++++++++++++--------- src/NeSST/utils.py | 42 +-- tests/test_core.py | 16 +- tests/test_endf.py | 16 +- tests/test_time_of_flight.py | 34 ++- tests/test_utils.py | 68 +++-- tools/table_build.py | 14 +- 27 files changed, 1543 insertions(+), 1170 deletions(-) create mode 100644 .pre-commit-config.yaml create mode 100644 ruff.toml diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index e96f09b..b4795a0 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -3,10 +3,10 @@ name: Upload Python Package on: release: types: [created] - + # Allows you to run this workflow manually from the Actions tab workflow_dispatch: - + jobs: deploy: runs-on: ubuntu-latest diff --git a/.gitignore b/.gitignore index 071b8f3..67c128b 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,7 @@ __pycache__ NeSST/__pycache__/ example/.ipynb_checkpoints/ *.egg-info -dist/ \ No newline at end of file +dist/ +.vscode/ + +notebook/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..8251763 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,50 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v5.0.0 + hooks: + - id: check-added-large-files + additional_dependencies: [--isolated] + args: ["--maxkb=2000"] + # Add exceptions here, as a regex + exclude: "" + + - id: check-json + additional_dependencies: [--isolated] + + - id: check-toml + additional_dependencies: [--isolated] + + - id: check-yaml + additional_dependencies: [--isolated] + + - id: detect-private-key + additional_dependencies: [--isolated] + + - id: end-of-file-fixer + additional_dependencies: [--isolated] + + - id: trailing-whitespace + additional_dependencies: [--isolated] + + +- repo: https://github.com/henryiii/validate-pyproject-schema-store + rev: 2025.05.12 + hooks: + - id: validate-pyproject + +- repo: https://github.com/astral-sh/ruff-pre-commit + # Ruff version. + rev: v0.11.12 + hooks: + # Run the linter. + - id: ruff + args: [--fix] + types_or: [pyi, python, jupyter] + # Ignore global python configuration for private registry and install hooks from public index + # Add for each hook + # Reference: https://github.com/pre-commit/pre-commit/issues/1454#issuecomment-1816328894 + additional_dependencies: [--isolated] + # Run the formatter. + - id: ruff-format + types_or: [pyi, python, jupyter] + additional_dependencies: [--isolated] diff --git a/README.md b/README.md index a50af57..bddc06b 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ E-mail: ac116@ic.ac.uk ## Installation -- Easier method: Install from PyPI +- Easier method: Install from PyPI ``` pip install NeSST diff --git a/pyproject.toml b/pyproject.toml index a422305..8c31281 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,7 +32,7 @@ dependencies = [ [project.optional-dependencies] test = [ - "pytest", "jupyter", "matplotlib" + "pre-commit","pytest", "jupyter", "matplotlib" ] [tool.setuptools.package-data] diff --git a/ruff.toml b/ruff.toml new file mode 100644 index 0000000..28c9b45 --- /dev/null +++ b/ruff.toml @@ -0,0 +1,49 @@ +# Set to the lowest supported Python version. +target-version = "py311" + +# Set the target line length for formatting. +line-length = 120 + +# Exclude a variety of commonly ignored directories. +extend-exclude = [ + "example/*.ipynb" +] + +src = ["."] + +[lint] +# Select and/or ignore rules for linting. +# Full list of available rules: https://docs.astral.sh/ruff/rules/ +extend-select = [ + "B", # Flake8 bugbear + "E", # Pycodestyle errors + "F", # Pyflakes + "I", # Isort + "NPY", # Numpy + "PT", # Pytest + "RUF", # Ruff-specific rules + "UP", # Pyupgrade + "W", # Pycodestyle warnings +] +ignore = [ + "NPY002", # Replace legacy `np.random.rand` call with `np.random.Generator` + "E731", # Do not assign a lambda expression, use a def + "F403", # * import used; unable to detect undefined names + "F405", # x may be undefined, or defined from star imports + "B007", # Loop control variable `i` not used within loop body + "UP006", # tuple, dict instead of Tuple, Dict + "UP015", # Unnecessary mode argument in open + "UP035", # typing.Tuple/Dict deprecated + "PT018", # Assertion should be broken down into multiple parts + "RUF002", # En dash + "B028", # No explicit `stacklevel` keyword argument found +] + + +[lint.pycodestyle] +max-line-length = 150 # Allow some flexibility in line lengths +max-doc-length = 150 + +[format] +# Enable reformatting of code snippets in docstrings. +docstring-code-format = true diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index 2bc8749..f1002f4 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,5 +1,3 @@ ################## from .constants import * from .core import * -import NeSST.fitting -import NeSST.time_of_flight \ No newline at end of file diff --git a/src/NeSST/collisions.py b/src/NeSST/collisions.py index b2e6400..b57b9f6 100644 --- a/src/NeSST/collisions.py +++ b/src/NeSST/collisions.py @@ -1,4 +1,5 @@ import numpy as np + from NeSST.constants import * classical_collisions = None @@ -7,174 +8,212 @@ # Relativistic Collisions # ########################### + def gamma(beta): - return 1.0/np.sqrt(1-beta**2) - -def p(m,beta): - g = gamma(beta) - return g*m*beta - -def E(m,beta): - mom = p(m,beta) - return np.sqrt(m**2+mom**2) - -def mom_invariant(m1,m2,beta1,beta2,cos): - E1 = E(m1,beta1) - E2 = E(m2,beta2) - p1 = p(m1,beta1) - p2 = p(m2,beta2) - return np.sqrt(m1**2+m2**2+2*E1*E2-2*p1*p2*cos) - -def rel_lab_scattering_cosine(m1,m2,beta1,beta2,beta3,cos12,cos23): - # m3 == m1 - E1 = E(m1,beta1); E2 = E(m2,beta2); E3 = E(m1,beta3) - p1 = p(m1,beta1); p2 = p(m2,beta2); p3 = p(m1,beta3) - W = mom_invariant(m1,m2,beta1,beta2,cos12) - mu0 = (E3*(E1+E2)-0.5*W**2-p2*p3*cos23+0.5*(m2**2-m1**2))/(p1*p3) - return mu0 + return 1.0 / np.sqrt(1 - beta**2) + + +def p(m, beta): + g = gamma(beta) + return g * m * beta + + +def E(m, beta): + mom = p(m, beta) + return np.sqrt(m**2 + mom**2) + + +def mom_invariant(m1, m2, beta1, beta2, cos): + E1 = E(m1, beta1) + E2 = E(m2, beta2) + p1 = p(m1, beta1) + p2 = p(m2, beta2) + return np.sqrt(m1**2 + m2**2 + 2 * E1 * E2 - 2 * p1 * p2 * cos) + + +def rel_lab_scattering_cosine(m1, m2, beta1, beta2, beta3, cos12, cos23): + # m3 == m1 + E1 = E(m1, beta1) + E2 = E(m2, beta2) + E3 = E(m1, beta3) + p1 = p(m1, beta1) + p2 = p(m2, beta2) + p3 = p(m1, beta3) + W = mom_invariant(m1, m2, beta1, beta2, cos12) + mu0 = (E3 * (E1 + E2) - 0.5 * W**2 - p2 * p3 * cos23 + 0.5 * (m2**2 - m1**2)) / (p1 * p3) + return mu0 + # With mu_in == +1, mu_out == mu_0 -def rel_mu_out(m1,m2,beta1,beta2,beta3): - # m3 == m1 - E1 = E(m1,beta1); E2 = E(m2,beta2); E3 = E(m1,beta3) - p1 = p(m1,beta1); p2 = p(m2,beta2); p3 = p(m1,beta3) - W = mom_invariant(m1,m2,beta1,beta2,1.0) - mu = (E3*(E1+E2)-0.5*W**2+0.5*(m2**2-m1**2))/(p1*p3+p2*p3) - return mu - -def rel_CoM_scattering_cosine(m1,m2,beta1,beta2,beta3,cos12,cos23): - # m3 == m1 - # Lab frame quantities - E1 = E(m1,beta1); E2 = E(m2,beta2); E3 = E(m1,beta3) - p1 = p(m1,beta1); p2 = p(m2,beta2); p3 = p(m1,beta3) - W = mom_invariant(m1,m2,beta1,beta2,cos12) - mu0 = (E3*(E1+E2)-0.5*W**2-p2*p3*cos23+0.5*(m2**2-m1**2))/(p1*p3) - # CoM quantities - betac = np.sqrt(p1**2+p2**2+2*p1*p2*cos12)/(E1+E2) - gammac = gamma(betac) - beta_p1 = (p1**2+p2*p1*cos12)/(E1+E2) - beta_p3 = (p1*p3*mu0+p2*p3*cos23)/(E1+E2) - # Evaluate - numerator = p1*p3*mu0+gammac**2*(beta_p1*beta_p3+betac**2*E1*E3-(E3*beta_p1+E1*beta_p3)) - denominator = 0.25*(W+(m1**2-m2**2)/W)**2-m1**2 - muc = numerator/denominator - return muc - -def rel_mucE_jacobian(m1,m2,beta1,beta2,beta3,cos12,cos23): - # m3 == m1 - # Lab frame quantities - E1 = E(m1,beta1); E2 = E(m2,beta2); E3 = E(m1,beta3) - p1 = p(m1,beta1); p2 = p(m2,beta2); p3 = p(m1,beta3) - W = mom_invariant(m1,m2,beta1,beta2,cos12) - # Evaluate - numerator = E2-p2*E3*cos23/p3 - denominator = 0.25*(W+(m1**2-m2**2)/W)**2-m1**2 - g = numerator/denominator - return g +def rel_mu_out(m1, m2, beta1, beta2, beta3): + # m3 == m1 + E1 = E(m1, beta1) + E2 = E(m2, beta2) + E3 = E(m1, beta3) + p1 = p(m1, beta1) + p2 = p(m2, beta2) + p3 = p(m1, beta3) + W = mom_invariant(m1, m2, beta1, beta2, 1.0) + mu = (E3 * (E1 + E2) - 0.5 * W**2 + 0.5 * (m2**2 - m1**2)) / (p1 * p3 + p2 * p3) + return mu + + +def rel_CoM_scattering_cosine(m1, m2, beta1, beta2, beta3, cos12, cos23): + # m3 == m1 + # Lab frame quantities + E1 = E(m1, beta1) + E2 = E(m2, beta2) + E3 = E(m1, beta3) + p1 = p(m1, beta1) + p2 = p(m2, beta2) + p3 = p(m1, beta3) + W = mom_invariant(m1, m2, beta1, beta2, cos12) + mu0 = (E3 * (E1 + E2) - 0.5 * W**2 - p2 * p3 * cos23 + 0.5 * (m2**2 - m1**2)) / (p1 * p3) + # CoM quantities + betac = np.sqrt(p1**2 + p2**2 + 2 * p1 * p2 * cos12) / (E1 + E2) + gammac = gamma(betac) + beta_p1 = (p1**2 + p2 * p1 * cos12) / (E1 + E2) + beta_p3 = (p1 * p3 * mu0 + p2 * p3 * cos23) / (E1 + E2) + # Evaluate + numerator = p1 * p3 * mu0 + gammac**2 * (beta_p1 * beta_p3 + betac**2 * E1 * E3 - (E3 * beta_p1 + E1 * beta_p3)) + denominator = 0.25 * (W + (m1**2 - m2**2) / W) ** 2 - m1**2 + muc = numerator / denominator + return muc + + +def rel_mucE_jacobian(m1, m2, beta1, beta2, beta3, cos12, cos23): + # m3 == m1 + # Lab frame quantities + E2 = E(m2, beta2) + E3 = E(m1, beta3) + p2 = p(m2, beta2) + p3 = p(m1, beta3) + W = mom_invariant(m1, m2, beta1, beta2, cos12) + # Evaluate + numerator = E2 - p2 * E3 * cos23 / p3 + denominator = 0.25 * (W + (m1**2 - m2**2) / W) ** 2 - m1**2 + g = numerator / denominator + return g + # Conversions -def Ekin_2_beta(Ek,m): - x = Ek/m+1 - beta = gamma_2_beta(x) - return beta +def Ekin_2_beta(Ek, m): + x = Ek / m + 1 + beta = gamma_2_beta(x) + return beta + def gamma_2_beta(g): - return np.sqrt(1.0-1.0/g**2) + return np.sqrt(1.0 - 1.0 / g**2) + def v_2_beta(v): - return v/c + return v / c + def beta_2_normtime(beta): - return 1.0/beta + return 1.0 / beta + -def beta_2_Ekin(beta,m): - Etot = E(m,beta) - return Etot-m +def beta_2_Ekin(beta, m): + Etot = E(m, beta) + return Etot - m -def Jacobian_dEdnorm_t(E,m): - beta = Ekin_2_beta(E,m) - gam = gamma(beta) - return m*(gam*beta)**3 -def velocity_addition_to_Ekin(Ek,m,u): - beta_frame = u/c - beta = Ekin_2_beta(Ek,m) - beta = (beta+beta_frame)/(1+beta*beta_frame) - Ek = beta_2_Ekin(beta,m) - return Ek +def Jacobian_dEdnorm_t(E, m): + beta = Ekin_2_beta(E, m) + gam = gamma(beta) + return m * (gam * beta) ** 3 + + +def velocity_addition_to_Ekin(Ek, m, u): + beta_frame = u / c + beta = Ekin_2_beta(Ek, m) + beta = (beta + beta_frame) / (1 + beta * beta_frame) + Ek = beta_2_Ekin(beta, m) + return Ek + ######################## # Classical Collisions # ######################## -def cla_lab_scattering_cosine(A,Ein,Eout,muin,muout,vf): - vout = sqrtE_2_v*np.sqrt(Eout) - vin = sqrtE_2_v*np.sqrt(Ein) - mu0_star = 0.5*((A+1)*np.sqrt(Eout/Ein)-(A-1)*np.sqrt(Ein/Eout)) - return mu0_star+A*vf/vout*muin-A*vf/vin*muout + +def cla_lab_scattering_cosine(A, Ein, Eout, muin, muout, vf): + vout = sqrtE_2_v * np.sqrt(Eout) + vin = sqrtE_2_v * np.sqrt(Ein) + mu0_star = 0.5 * ((A + 1) * np.sqrt(Eout / Ein) - (A - 1) * np.sqrt(Ein / Eout)) + return mu0_star + A * vf / vout * muin - A * vf / vin * muout + # With mu_in == +1, mu_out == mu_0 -def cla_mu_out(A,Ein,Eout,vf): - vout = sqrtE_2_v*np.sqrt(Eout) - vin = sqrtE_2_v*np.sqrt(Ein) - mu0_star = 0.5*((A+1)*np.sqrt(Eout/Ein)-(A-1)*np.sqrt(Ein/Eout)) - return (mu0_star+A*vf/vout)/(1+A*vf/vin) - -def cla_CoM_scattering_cosine(A,Ein,Eout,muin,muout,vf): - vout = sqrtE_2_v*np.sqrt(Eout) - vin = sqrtE_2_v*np.sqrt(Ein) - v_ratio = (vout**2-2*vf*vout*muout+vf**2)/(vin**2-2*vf*vin*muin+vf**2) - return (A+1)**2*v_ratio/(2*A)-(A**2+1)/(2*A) - -def cla_mucE_jacobian(A,Ein,Eout,muin,muout,vf): - vout = sqrtE_2_v*np.sqrt(Eout) - vin = sqrtE_2_v*np.sqrt(Ein) - alpha = ((A-1)/(A+1))**2 - g0 = 2.0/((1-alpha)*Ein) - vcorr = ((1-vf*muout/vout)/(1-2*vf*muin/vin+(vf/vin)**2)) - return g0*vcorr +def cla_mu_out(A, Ein, Eout, vf): + vout = sqrtE_2_v * np.sqrt(Eout) + vin = sqrtE_2_v * np.sqrt(Ein) + mu0_star = 0.5 * ((A + 1) * np.sqrt(Eout / Ein) - (A - 1) * np.sqrt(Ein / Eout)) + return (mu0_star + A * vf / vout) / (1 + A * vf / vin) + + +def cla_CoM_scattering_cosine(A, Ein, Eout, muin, muout, vf): + vout = sqrtE_2_v * np.sqrt(Eout) + vin = sqrtE_2_v * np.sqrt(Ein) + v_ratio = (vout**2 - 2 * vf * vout * muout + vf**2) / (vin**2 - 2 * vf * vin * muin + vf**2) + return (A + 1) ** 2 * v_ratio / (2 * A) - (A**2 + 1) / (2 * A) + + +def cla_mucE_jacobian(A, Ein, Eout, muin, muout, vf): + vout = sqrtE_2_v * np.sqrt(Eout) + vin = sqrtE_2_v * np.sqrt(Ein) + alpha = ((A - 1) / (A + 1)) ** 2 + g0 = 2.0 / ((1 - alpha) * Ein) + vcorr = (1 - vf * muout / vout) / (1 - 2 * vf * muin / vin + (vf / vin) ** 2) + return g0 * vcorr + ####################### # Interface functions # ####################### + # Change in flux due to different relative velocity between target and scatterer # |vn-vf|/vn -def flux_change(Ein,muin,vf): - vin = sqrtE_2_v*np.sqrt(Ein) - del_f = np.sqrt(1-2*vf*muin/vin+(vf/vin)**2) - return del_f +def flux_change(Ein, muin, vf): + vin = sqrtE_2_v * np.sqrt(Ein) + del_f = np.sqrt(1 - 2 * vf * muin / vin + (vf / vin) ** 2) + return del_f + # Slowing down kernel -def g(A,Ein,Eout,muin,muout,vf): - if (classical_collisions == True): - ans = cla_mucE_jacobian(A,Ein,Eout,muin,muout,vf) - else: - beta1 = Ekin_2_beta(Ein,Mn) - beta2 = v_2_beta(vf) - beta3 = Ekin_2_beta(Eout,Mn) - ans = rel_mucE_jacobian(Mn,A*Mn,beta1,beta2,beta3,muin,muout) - return ans +def g(A, Ein, Eout, muin, muout, vf): + if classical_collisions: + ans = cla_mucE_jacobian(A, Ein, Eout, muin, muout, vf) + else: + beta1 = Ekin_2_beta(Ein, Mn) + beta2 = v_2_beta(vf) + beta3 = Ekin_2_beta(Eout, Mn) + ans = rel_mucE_jacobian(Mn, A * Mn, beta1, beta2, beta3, muin, muout) + return ans + # Centre of mass cosine -def muc(A,Ein,Eout,muin,muout,vf): - if (classical_collisions == True): - ans = cla_CoM_scattering_cosine(A,Ein,Eout,muin,muout,vf) - else: - beta1 = Ekin_2_beta(Ein,Mn) - beta2 = v_2_beta(vf) - beta3 = Ekin_2_beta(Eout,Mn) - ans = rel_CoM_scattering_cosine(Mn,A*Mn,beta1,beta2,beta3,muin,muout) - return ans +def muc(A, Ein, Eout, muin, muout, vf): + if classical_collisions: + ans = cla_CoM_scattering_cosine(A, Ein, Eout, muin, muout, vf) + else: + beta1 = Ekin_2_beta(Ein, Mn) + beta2 = v_2_beta(vf) + beta3 = Ekin_2_beta(Eout, Mn) + ans = rel_CoM_scattering_cosine(Mn, A * Mn, beta1, beta2, beta3, muin, muout) + return ans -# Outgoing neutron cosine -def mu_out(A,Ein,Eout,vf): - if (classical_collisions == True): - ans = cla_mu_out(A,Ein,Eout,vf) - else: - beta1 = Ekin_2_beta(Ein,Mn) - beta2 = v_2_beta(vf) - beta3 = Ekin_2_beta(Eout,Mn) - ans = rel_mu_out(Mn,A*Mn,beta1,beta2,beta3) - return ans +# Outgoing neutron cosine +def mu_out(A, Ein, Eout, vf): + if classical_collisions: + ans = cla_mu_out(A, Ein, Eout, vf) + else: + beta1 = Ekin_2_beta(Ein, Mn) + beta2 = v_2_beta(vf) + beta3 = Ekin_2_beta(Eout, Mn) + ans = rel_mu_out(Mn, A * Mn, beta1, beta2, beta3) + return ans diff --git a/src/NeSST/constants.py b/src/NeSST/constants.py index 6703b63..cb7fb68 100644 --- a/src/NeSST/constants.py +++ b/src/NeSST/constants.py @@ -1,37 +1,37 @@ -import scipy.constants as sc +import os + import numpy as np -import numpy.typing as npt +import scipy.constants as sc -''' Contains some physical constants used throughout the analysis ''' +""" Contains some physical constants used throughout the analysis """ # Scipy constants uses CODATA2018 database # Need to swap from MeV to eV -amu = sc.value('atomic mass constant energy equivalent in MeV')*1e6 -c = sc.c -Me = sc.value('electron mass energy equivalent in MeV')*1e6 -Mn = sc.value('neutron mass energy equivalent in MeV')*1e6 -Mp = sc.value('proton mass energy equivalent in MeV')*1e6 -Md = sc.value('deuteron mass energy equivalent in MeV')*1e6 -Mt = sc.value('triton mass energy equivalent in MeV')*1e6 -MHe3 = sc.value('helion mass energy equivalent in MeV')*1e6 -MHe4 = sc.value('alpha particle mass energy equivalent in MeV')*1e6 -MC = 12.011*amu -MBe = 9.012182*amu +amu = sc.value("atomic mass constant energy equivalent in MeV") * 1e6 +c = sc.c +Me = sc.value("electron mass energy equivalent in MeV") * 1e6 +Mn = sc.value("neutron mass energy equivalent in MeV") * 1e6 +Mp = sc.value("proton mass energy equivalent in MeV") * 1e6 +Md = sc.value("deuteron mass energy equivalent in MeV") * 1e6 +Mt = sc.value("triton mass energy equivalent in MeV") * 1e6 +MHe3 = sc.value("helion mass energy equivalent in MeV") * 1e6 +MHe4 = sc.value("alpha particle mass energy equivalent in MeV") * 1e6 +MC = 12.011 * amu +MBe = 9.012182 * amu hbar = sc.hbar -qe = sc.e -fine_structure = sc.fine_structure -sqrtE_2_v = np.sqrt(2*sc.e/sc.m_n) +qe = sc.e +fine_structure = sc.fine_structure +sqrtE_2_v = np.sqrt(2 * sc.e / sc.m_n) Mn_kg = sc.m_n -sigmabarn = 1e-28 # barns to m^2 -E0_DT = ((Md+Mt)**2+Mn**2-MHe4**2)/(2*(Md+Mt))-Mn -E0_DD = ((Md+Md)**2+Mn**2-MHe3**2)/(2*(Md+Md))-Mn +sigmabarn = 1e-28 # barns to m^2 +E0_DT = ((Md + Mt) ** 2 + Mn**2 - MHe4**2) / (2 * (Md + Mt)) - Mn +E0_DD = ((Md + Md) ** 2 + Mn**2 - MHe3**2) / (2 * (Md + Md)) - Mn # Directories -import os package_directory = os.path.dirname(os.path.abspath(__file__)) -data_dir = os.path.join(package_directory,"./data/") -ENDF_dir = os.path.join(data_dir,"./ENDF/") +data_dir = os.path.join(package_directory, "./data/") +ENDF_dir = os.path.join(data_dir, "./ENDF/") # Materials -default_mat_list = ['H','D','T','C12','Be9'] +default_mat_list = ["H", "D", "T", "C12", "Be9"] available_materials = [] -mat_dict = {} \ No newline at end of file +mat_dict = {} diff --git a/src/NeSST/core.py b/src/NeSST/core.py index 1e24a7b..4c94559 100644 --- a/src/NeSST/core.py +++ b/src/NeSST/core.py @@ -4,16 +4,19 @@ # Standard libraries import typing -import numpy.typing as npt import warnings -import numpy as np from glob import glob from os.path import basename -# NeSST libraries -from NeSST.constants import * + +import numpy as np +import numpy.typing as npt + import NeSST.collisions as col import NeSST.spectral_model as sm +# NeSST libraries +from NeSST.constants import * + # Global variable defaults col.classical_collisions = False # Atomic fraction of D and T in scattering medium and source @@ -24,30 +27,32 @@ # Energies, temperatures eV # Velocities m/s + ############################ # Default material load in # ############################ def initialise_material_data(label): # Aliases - if(label == 'H'): - json = 'H1.json' - elif(label == 'D'): - json = 'H2.json' - elif(label == 'T'): - json = 'H3.json' + if label == "H": + json = "H1.json" + elif label == "D": + json = "H2.json" + elif label == "T": + json = "H3.json" # Parse name else: - mat_jsons = glob(data_dir+'*.json') - mats = [basename(f).split('.')[0] for f in mat_jsons] - if(label in mats): - json = label+'.json' + mat_jsons = glob(data_dir + "*.json") + mats = [basename(f).split(".")[0] for f in mat_jsons] + if label in mats: + json = label + ".json" else: - print("Material label '"+label+"' not recognised") + print("Material label '" + label + "' not recognised") return - mat_data = sm.material_data(label,json) + mat_data = sm.material_data(label, json) mat_dict[label] = mat_data available_materials.append(label) + for mat in default_mat_list: initialise_material_data(mat) @@ -55,8 +60,9 @@ def initialise_material_data(label): # Primary spectral shapes & reactivities # ########################################## + # Gaussian "Brysk" -def QBrysk(Ein : npt.NDArray, mean : float, variance : float) -> npt.NDArray: +def QBrysk(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray: """Calculates the primary spectrum with a Brysk shape i.e. Gaussian Args: Ein (numpy.array) : array of energy values on which to compute spectrum @@ -65,13 +71,14 @@ def QBrysk(Ein : npt.NDArray, mean : float, variance : float) -> npt.NDArray: Returns: numpy.array : array with Gaussian spectrum on array Ein - + """ - spec = np.exp(-(Ein-mean)**2/2.0/variance)/np.sqrt(2*np.pi*variance) + spec = np.exp(-((Ein - mean) ** 2) / 2.0 / variance) / np.sqrt(2 * np.pi * variance) return spec + # Ballabio -def QBallabio(Ein : npt.NDArray, mean : float, variance : float) -> npt.NDArray: +def QBallabio(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray: """Calculates the primary spectrum with a Ballabio shape i.e. modified Gaussian See equations 44 - 46 of Ballabio et al. @@ -82,18 +89,19 @@ def QBallabio(Ein : npt.NDArray, mean : float, variance : float) -> npt.NDArray: Returns: numpy.array : array with modified Gaussian spectrum on array Ein - + """ - common_factor = 1-1.5*variance/mean**2 - Ebar = mean*np.sqrt(common_factor) - sig2 = 4.0/3.0*mean**2*(np.sqrt(common_factor)-common_factor) - norm = np.sqrt(2*np.pi*variance) - spec = np.exp(-2.0*Ebar*(np.sqrt(Ein)-np.sqrt(Ebar))**2/sig2)/norm + common_factor = 1 - 1.5 * variance / mean**2 + Ebar = mean * np.sqrt(common_factor) + sig2 = 4.0 / 3.0 * mean**2 * (np.sqrt(common_factor) - common_factor) + norm = np.sqrt(2 * np.pi * variance) + spec = np.exp(-2.0 * Ebar * (np.sqrt(Ein) - np.sqrt(Ebar)) ** 2 / sig2) / norm return spec + # TT spectral shape -def dNdE_TT(E : npt.NDArray, Tion : float) -> npt.NDArray: - """Calculates the TT primary spectrum with Doppler broadening effect as +def dNdE_TT(E: npt.NDArray, Tion: float) -> npt.NDArray: + """Calculates the TT primary spectrum with Doppler broadening effect as calculated in Appelbe et al. HEDP 2016 Args: @@ -102,13 +110,15 @@ def dNdE_TT(E : npt.NDArray, Tion : float) -> npt.NDArray: Returns: numpy.array : array with normalised TT spectral shape at energies E - + """ - return sm.TT_2dinterp(E,Tion) + return sm.TT_2dinterp(E, Tion) + -def yield_from_dt_yield_ratio(reaction : str, dt_yield : float, Tion : float, - frac_D: float = frac_D_default, frac_T: float = frac_T_default) -> float: - """ Reactivity ratio to predict yield from the DT yield assuming same volume and burn time +def yield_from_dt_yield_ratio( + reaction: str, dt_yield: float, Tion: float, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> float: + """Reactivity ratio to predict yield from the DT yield assuming same volume and burn time rate_ij = (f_{i}*f_{j}*sigmav_{i,j}(T))/(1+delta_{i,j}) # dN/dVdt yield_ij = (rate_ij/rate_dt)*yield_dt @@ -132,25 +142,26 @@ def yield_from_dt_yield_ratio(reaction : str, dt_yield : float, Tion : float, if Tion < 0: raise ValueError("Tion (temperature of the ions) can not be below 0") - - if sum([frac_D, frac_T]) != 1.: - msg = (f'The frac_D ({frac_D_default}) and frac_T ({frac_T_default}) ' - 'arguments on the yield_from_dt_yield_ration method do not sum to 1.') + + if sum([frac_D, frac_T]) != 1.0: + msg = f"The frac_D ({frac_D_default}) and frac_T ({frac_T_default}) arguments on the yield_from_dt_yield_ration method do not sum to 1." warnings.warn(msg) - if reaction == 'tt': - ratio = (0.5*frac_T*sm.reac_TT(Tion))/(frac_D*sm.reac_DT(Tion)) - ratio = 2.* ratio # Two neutrons are generated for each reaction - elif reaction == 'dd': - ratio = (0.5*frac_D*sm.reac_DD(Tion))/(frac_T*sm.reac_DT(Tion)) + if reaction == "tt": + ratio = (0.5 * frac_T * sm.reac_TT(Tion)) / (frac_D * sm.reac_DT(Tion)) + ratio = 2.0 * ratio # Two neutrons are generated for each reaction + elif reaction == "dd": + ratio = (0.5 * frac_D * sm.reac_DD(Tion)) / (frac_T * sm.reac_DT(Tion)) else: raise ValueError(f'reaction should be either "dd" or "tt" not {reaction}') - return ratio*dt_yield + return ratio * dt_yield + -def yields_normalised(Tion : float, frac_D: float = frac_D_default, frac_T: float = frac_T_default - ) -> typing.Tuple[float, float, float]: - """ Assuming same volume and burn time, find fractional yields of DT, DD and TT respectively +def yields_normalised( + Tion: float, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> typing.Tuple[float, float, float]: + """Assuming same volume and burn time, find fractional yields of DT, DD and TT respectively Uses default models for reactivities Note that the TT reaction produces two neutrons. @@ -170,24 +181,25 @@ def yields_normalised(Tion : float, frac_D: float = frac_D_default, frac_T: floa if Tion < 0: raise ValueError("Tion (temperature of the ions) can not be below 0") - - if sum([frac_D, frac_T]) != 1.: - msg = (f'The frac_D ({frac_D_default}) and frac_T ({frac_T_default}) ' - 'arguments on the yield_from_dt_yield_ration method do not sum to 1.') + + if sum([frac_D, frac_T]) != 1.0: + msg = f"The frac_D ({frac_D_default}) and frac_T ({frac_T_default}) arguments on the yield_from_dt_yield_ration method do not sum to 1." warnings.warn(msg) - unnormed_dt_yield = frac_D*frac_T*sm.reac_DT(Tion) - unnormed_dd_yield = 0.5*frac_D*frac_D*sm.reac_DD(Tion) - unnormed_tt_yield = 2.0*0.5*frac_T*frac_T*sm.reac_TT(Tion) + unnormed_dt_yield = frac_D * frac_T * sm.reac_DT(Tion) + unnormed_dd_yield = 0.5 * frac_D * frac_D * sm.reac_DD(Tion) + unnormed_tt_yield = 2.0 * 0.5 * frac_T * frac_T * sm.reac_TT(Tion) - tot_yield = unnormed_dt_yield+unnormed_dd_yield+unnormed_tt_yield + tot_yield = unnormed_dt_yield + unnormed_dd_yield + unnormed_tt_yield + + return unnormed_dt_yield / tot_yield, unnormed_dd_yield / tot_yield, unnormed_tt_yield / tot_yield - return unnormed_dt_yield/tot_yield,unnormed_dd_yield/tot_yield,unnormed_tt_yield/tot_yield ############################################################################### # Ballabio fits, see Table III of L. Ballabio et al 1998 Nucl. Fusion 38 1723 # ############################################################################### + # Returns the mean and variance based on Ballabio def DTprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: """Calculates the mean energy and the variance of the neutron energy @@ -215,9 +227,7 @@ def DTprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: a4 = 1.3818 Tion_kev = Tion / 1e3 # Ballabio equation accepts KeV units - mean_shift = ( - a1 * Tion_kev ** (0.6666666666) / (1.0 + a2 * Tion_kev**a3) + a4 * Tion_kev - ) + mean_shift = a1 * Tion_kev ** (0.6666666666) / (1.0 + a2 * Tion_kev**a3) + a4 * Tion_kev mean_shift *= 1e3 # converting back to eV mean = E0_DT + mean_shift @@ -239,6 +249,7 @@ def DTprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: return mean, stddev, variance + # Returns the mean and variance based on Ballabio def DDprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: """Calculates the mean energy and the variance of the neutron energy @@ -266,9 +277,7 @@ def DDprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: a4 = 0.81844 Tion_kev = Tion / 1e3 # Ballabio equation accepts KeV units - mean_shift = ( - a1 * Tion_kev ** (0.6666666666) / (1.0 + a2 * Tion_kev**a3) + a4 * Tion_kev - ) + mean_shift = a1 * Tion_kev ** (0.6666666666) / (1.0 + a2 * Tion_kev**a3) + a4 * Tion_kev mean_shift *= 1e3 # converting back to eV mean = E0_DD + mean_shift @@ -289,28 +298,32 @@ def DDprimspecmoments(Tion: float) -> typing.Tuple[float, float, float]: return mean, stddev, variance -def neutron_velocity_addition(Ek,u): - return col.velocity_addition_to_Ekin(Ek,Mn,u) + +def neutron_velocity_addition(Ek, u): + return col.velocity_addition_to_Ekin(Ek, Mn, u) + ####################################### # DT scattered spectra initialisation # ####################################### -def init_DT_scatter(Eout : npt.NDArray, Ein : npt.NDArray): + +def init_DT_scatter(Eout: npt.NDArray, Ein: npt.NDArray): """Initialise the scattering matrices for D and T materials Args: Ein (numpy.array): the array on incoming neutron energies Eout (numpy.array): the array on outgoing neutron energies - + """ - mat_dict['D'].init_energy_grids(Eout,Ein) - mat_dict['T'].init_energy_grids(Eout,Ein) - mat_dict['D'].init_station_scatter_matrices() - mat_dict['T'].init_station_scatter_matrices() + mat_dict["D"].init_energy_grids(Eout, Ein) + mat_dict["T"].init_energy_grids(Eout, Ein) + mat_dict["D"].init_station_scatter_matrices() + mat_dict["T"].init_station_scatter_matrices() + -def init_DT_ionkin_scatter(varr : npt.NDArray, nT: bool = False, nD: bool = False): - """Initialise the scattering matrices including the effect of ion +def init_DT_ionkin_scatter(varr: npt.NDArray, nT: bool = False, nD: bool = False): + """Initialise the scattering matrices including the effect of ion velocities in the kinematics N.B. the static ion scattering matrices must already be calculated @@ -319,40 +332,41 @@ def init_DT_ionkin_scatter(varr : npt.NDArray, nT: bool = False, nD: bool = Fals Args: Ein (numpy.array): the array on incoming neutron energies Eout (numpy.array): the array on outgoing neutron energies - + """ - if(nT): - if(mat_dict['T'].Ein is None): + if nT: + if mat_dict["T"].Ein is None: print("nT - Needed to initialise energy grids - see init_DT_scatter") else: - mat_dict['T'].full_scattering_matrix_create(varr) - if(nD): - if(mat_dict['D'].Ein is None): + mat_dict["T"].full_scattering_matrix_create(varr) + if nD: + if mat_dict["D"].Ein is None: print("nD - Needed to initialise energy grids - see init_DT_scatter") else: - mat_dict['D'].full_scattering_matrix_create(varr) + mat_dict["D"].full_scattering_matrix_create(varr) + -def calc_DT_ionkin_primspec_rhoL_integral(I_E : npt.NDArray, rhoL_func=None, nT: bool = False, nD: bool = False): - if(nT): - if(mat_dict['T'].vvec is None): +def calc_DT_ionkin_primspec_rhoL_integral(I_E: npt.NDArray, rhoL_func=None, nT: bool = False, nD: bool = False): + if nT: + if mat_dict["T"].vvec is None: print("nT - Needed to initialise velocity grid - see init_DT_ionkin_scatter") else: - if(rhoL_func is not None): - mat_dict['T'].scattering_matrix_apply_rhoLfunc(rhoL_func) - mat_dict['T'].matrix_primspec_int(I_E) - if(nD): - if(mat_dict['D'].vvec is None): + if rhoL_func is not None: + mat_dict["T"].scattering_matrix_apply_rhoLfunc(rhoL_func) + mat_dict["T"].matrix_primspec_int(I_E) + if nD: + if mat_dict["D"].vvec is None: print("nD - Needed to initialise velocity grid - see init_DT_ionkin_scatter") else: - if(rhoL_func is not None): - mat_dict['D'].scattering_matrix_apply_rhoLfunc(rhoL_func) - mat_dict['D'].matrix_primspec_int(I_E) + if rhoL_func is not None: + mat_dict["D"].scattering_matrix_apply_rhoLfunc(rhoL_func) + mat_dict["D"].matrix_primspec_int(I_E) ################################### # General material initialisation # ################################### -def init_mat_scatter(Eout : npt.NDArray, Ein : npt.NDArray, mat_label : str): +def init_mat_scatter(Eout: npt.NDArray, Ein: npt.NDArray, mat_label: str): """General material version of init_DT_scatter as specified by material label N.B. the mat_lable must match those in available_materials_dict @@ -361,219 +375,232 @@ def init_mat_scatter(Eout : npt.NDArray, Ein : npt.NDArray, mat_label : str): Ein (numpy.array): the array on incoming neutron energies Eout (numpy.array): the array on outgoing neutron energies mat_label (str) : material label - - + + """ mat = mat_dict[mat_label] - mat.init_energy_grids(Eout,Ein) + mat.init_energy_grids(Eout, Ein) mat.init_station_scatter_matrices() return mat + ########################################### # Single Evalutation Scattered Spectra # ########################################### -def DT_sym_scatter_spec(I_E : npt.NDArray - ,frac_D: float = frac_D_default, frac_T: float = frac_T_default - ) -> typing.Tuple[npt.NDArray, - typing.Tuple[npt.NDArray,npt.NDArray,npt.NDArray,npt.NDArray]]: - """Calculates the single scattered neutron spectrum for DT given a + +def DT_sym_scatter_spec( + I_E: npt.NDArray, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> typing.Tuple[npt.NDArray, typing.Tuple[npt.NDArray, npt.NDArray, npt.NDArray, npt.NDArray]]: + """Calculates the single scattered neutron spectrum for DT given a primary neutron spectrum of I_E from isotropic areal density This requires the scattering matrices to have been pre-calculated - The primary neutron spectrum, I_E, is assumed to be on the same energy grid as + The primary neutron spectrum, I_E, is assumed to be on the same energy grid as the incoming energy grid used to calculate the scattering matrices Args: I_E (numpy.array): the neutron spectrum at Ein energies Returns: - Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components + Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components (nD,nT,Dn2n,Tn2n) - + """ - rhoL_func = lambda x : np.ones_like(x) - mat_dict['D'].calc_dNdEs(I_E,rhoL_func) - mat_dict['T'].calc_dNdEs(I_E,rhoL_func) - nD = frac_D*mat_dict['D'].elastic_dNdE - nT = frac_T*mat_dict['T'].elastic_dNdE - Dn2n = frac_D*mat_dict['D'].n2n_dNdE - Tn2n = frac_T*mat_dict['T'].n2n_dNdE - total = nD+nT+Dn2n+Tn2n - return total,(nD,nT,Dn2n,Tn2n) - -def DT_asym_scatter_spec(I_E : npt.NDArray, rhoL_func : callable, - frac_D: float = frac_D_default, frac_T: float = frac_T_default - ) -> typing.Tuple[npt.NDArray, - typing.Tuple[npt.NDArray,npt.NDArray,npt.NDArray,npt.NDArray]]: - """Calculates the single scattered neutron spectrum for DT given a + rhoL_func = lambda x: np.ones_like(x) + mat_dict["D"].calc_dNdEs(I_E, rhoL_func) + mat_dict["T"].calc_dNdEs(I_E, rhoL_func) + nD = frac_D * mat_dict["D"].elastic_dNdE + nT = frac_T * mat_dict["T"].elastic_dNdE + Dn2n = frac_D * mat_dict["D"].n2n_dNdE + Tn2n = frac_T * mat_dict["T"].n2n_dNdE + total = nD + nT + Dn2n + Tn2n + return total, (nD, nT, Dn2n, Tn2n) + + +def DT_asym_scatter_spec( + I_E: npt.NDArray, rhoL_func: callable, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> typing.Tuple[npt.NDArray, typing.Tuple[npt.NDArray, npt.NDArray, npt.NDArray, npt.NDArray]]: + """Calculates the single scattered neutron spectrum for DT given a primary neutron spectrum of I_E from anisotropic areal density This requires the scattering matrices to have been pre-calculated - The primary neutron spectrum, I_E, is assumed to be on the same energy grid as + The primary neutron spectrum, I_E, is assumed to be on the same energy grid as the incoming energy grid used to calculate the scattering matrices - The areal density function rhoL_func needs to be a callable function with a + The areal density function rhoL_func needs to be a callable function with a single argument (cosine[theta]) Args: I_E (numpy.array): the neutron spectrum at Ein energies - rhoL_func (callable): must be a single argument function f(x), + rhoL_func (callable): must be a single argument function f(x), where x e [-1,1] and f(x) e [0,inf] and int f(x) dx = 1 frac_D (float) : fraction of D in fuel frac_T (float) : fraction of T in fuel Returns: - Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components + Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components (nD,nT,Dn2n,Tn2n) - + """ - mat_dict['D'].calc_dNdEs(I_E,rhoL_func) - mat_dict['T'].calc_dNdEs(I_E,rhoL_func) - nD = frac_D*mat_dict['D'].elastic_dNdE - nT = frac_T*mat_dict['T'].elastic_dNdE - Dn2n = frac_D*mat_dict['D'].n2n_dNdE - Tn2n = frac_T*mat_dict['T'].n2n_dNdE - total = nD+nT+Dn2n+Tn2n - return total,(nD,nT,Dn2n,Tn2n) - -def DT_scatter_spec_w_ionkin(I_E : npt.NDArray, vbar : float, dv : float, rhoL_func : callable, - frac_D: float = frac_D_default, frac_T: float = frac_T_default - ) -> typing.Tuple[npt.NDArray, - typing.Tuple[npt.NDArray,npt.NDArray,npt.NDArray,npt.NDArray]]: - """Calculates the single scattered neutron spectrum for DT given a - primary neutron spectrum of I_E from anisotropic areal density and including + mat_dict["D"].calc_dNdEs(I_E, rhoL_func) + mat_dict["T"].calc_dNdEs(I_E, rhoL_func) + nD = frac_D * mat_dict["D"].elastic_dNdE + nT = frac_T * mat_dict["T"].elastic_dNdE + Dn2n = frac_D * mat_dict["D"].n2n_dNdE + Tn2n = frac_T * mat_dict["T"].n2n_dNdE + total = nD + nT + Dn2n + Tn2n + return total, (nD, nT, Dn2n, Tn2n) + + +def DT_scatter_spec_w_ionkin( + I_E: npt.NDArray, + vbar: float, + dv: float, + rhoL_func: callable, + frac_D: float = frac_D_default, + frac_T: float = frac_T_default, +) -> typing.Tuple[npt.NDArray, typing.Tuple[npt.NDArray, npt.NDArray, npt.NDArray, npt.NDArray]]: + """Calculates the single scattered neutron spectrum for DT given a + primary neutron spectrum of I_E from anisotropic areal density and including ion velocities kinematics This requires the scattering matrices with ion kinematics to have been pre-calculated - The primary neutron spectrum, I_E, is assumed to be on the same energy grid as + The primary neutron spectrum, I_E, is assumed to be on the same energy grid as the incoming energy grid used to calculate the scattering matrices - The areal density function rhoL_func needs to be a callable function with a + The areal density function rhoL_func needs to be a callable function with a single argument (cosine[theta]) Args: I_E (numpy.array): the neutron spectrum at Ein energies vbar (float) : mean velocity of the scattering ions dv (float) : standard deviation velocity of the scattering ions - rhoL_func (callable): must be a single argument function f(x), + rhoL_func (callable): must be a single argument function f(x), where x e [-1,1] and f(x) e [0,inf] and int f(x) dx = 1 frac_D (float) : fraction of D in fuel frac_T (float) : fraction of T in fuel Returns: - Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components + Tuple of numpy.arrays: the total scattered spectrum and a tuple of the components (nD,nT,Dn2n,Tn2n) - + """ - rhoL_func = lambda x : np.ones_like(x) - mat_dict['D'].calc_dNdEs(I_E,rhoL_func) - mat_dict['T'].calc_dNdEs(I_E,rhoL_func) - if(mat_dict['D'].vvec is None): - dNdE_nD = mat_dict['D'].elastic_dNdE + rhoL_func = lambda x: np.ones_like(x) + mat_dict["D"].calc_dNdEs(I_E, rhoL_func) + mat_dict["T"].calc_dNdEs(I_E, rhoL_func) + if mat_dict["D"].vvec is None: + dNdE_nD = mat_dict["D"].elastic_dNdE else: - dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(mat_dict['D'].Eout,vbar,dv) - if(mat_dict['T'].vvec is None): - dNdE_nT = mat_dict['T'].elastic_dNdE + dNdE_nD = mat_dict["D"].matrix_interpolate_gaussian(mat_dict["D"].Eout, vbar, dv) + if mat_dict["T"].vvec is None: + dNdE_nT = mat_dict["T"].elastic_dNdE else: - dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(mat_dict['T'].Eout,vbar,dv) - nD = frac_D*dNdE_nD - nT = frac_T*dNdE_nT - Dn2n = frac_D*mat_dict['D'].n2n_dNdE - Tn2n = frac_T*mat_dict['T'].n2n_dNdE - total = nD+nT+Dn2n+Tn2n - return total,(nD,nT,Dn2n,Tn2n) - -def DT_transmission(rhoL : float, E_in : npt.NDArray, rhoL_func : callable, - frac_D: float = frac_D_default, frac_T: float = frac_T_default) -> npt.NDArray: + dNdE_nT = mat_dict["T"].matrix_interpolate_gaussian(mat_dict["T"].Eout, vbar, dv) + nD = frac_D * dNdE_nD + nT = frac_T * dNdE_nT + Dn2n = frac_D * mat_dict["D"].n2n_dNdE + Tn2n = frac_T * mat_dict["T"].n2n_dNdE + total = nD + nT + Dn2n + Tn2n + return total, (nD, nT, Dn2n, Tn2n) + + +def DT_transmission( + rhoL: float, E_in: npt.NDArray, rhoL_func: callable, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> npt.NDArray: """ Calculates the straight line transmission of primary fusion sources through the DT areal density - The areal density function rhoL_func needs to be a callable function with a + The areal density function rhoL_func needs to be a callable function with a single argument (cosine[theta]) Args: rhoL (float): the (4-pi averaged) areal density of the DT E_in (numpy.array): the energy array of the primary neutron spectrum - rhoL_func (callable): must be a single argument function f(x), + rhoL_func (callable): must be a single argument function f(x), where x e [-1,1] and f(x) e [0,inf] and int f(x) dx = 1 frac_D (float) : fraction of D in fuel frac_T (float) : fraction of T in fuel Returns: numpy.array: the total transmission coefficient as calculated - + exp[ - A_1S rhoL_func(1) (f_D sigma_{D,tot}(Ein) + f_T sigma_{T,tot}(Ein)) ] """ - A_1S = rhoR_2_A1s(rhoL,frac_D=frac_D,frac_T=frac_T) - tot_xsec = frac_D*mat_dict['D'].sigma_tot(E_in) + frac_T*mat_dict['T'].sigma_tot(E_in) - return np.exp(-A_1S*rhoL_func(1.0)*tot_xsec) + A_1S = rhoR_2_A1s(rhoL, frac_D=frac_D, frac_T=frac_T) + tot_xsec = frac_D * mat_dict["D"].sigma_tot(E_in) + frac_T * mat_dict["T"].sigma_tot(E_in) + return np.exp(-A_1S * rhoL_func(1.0) * tot_xsec) -def mat_scatter_spec(mat : typing.Type[sm.material_data], - I_E : npt.NDArray, rhoL_func : callable) -> npt.NDArray: - """Calculates a material's single scattered neutron spectrum given a + +def mat_scatter_spec(mat: typing.Type[sm.material_data], I_E: npt.NDArray, rhoL_func: callable) -> npt.NDArray: + """Calculates a material's single scattered neutron spectrum given a primary neutron spectrum of I_E from anisotropic areal density This requires the scattering matrices to have been pre-calculated - The primary neutron spectrum, I_E, is assumed to be on the same energy grid as + The primary neutron spectrum, I_E, is assumed to be on the same energy grid as the incoming energy grid used to calculate the scattering matrices - The areal density function rhoL_func needs to be a callable function with a + The areal density function rhoL_func needs to be a callable function with a single argument (cosine[theta]) Args: I_E (numpy.array): the neutron spectrum at Ein energies - rhoL_func (callable): must be a single argument function f(x), + rhoL_func (callable): must be a single argument function f(x), where x e [-1,1] and f(x) e [0,inf] and int f(x) dx = 1 Returns: numpy.array: the total scattered spectrum - + """ - mat.calc_dNdEs(I_E,rhoL_func) + mat.calc_dNdEs(I_E, rhoL_func) total = mat.elastic_dNdE.copy() - if(mat.l_n2n): + if mat.l_n2n: total += mat.n2n_dNdE - if(mat.l_inelastic): + if mat.l_inelastic: total += mat.inelastic_dNdE return total -def mat_transmission(mat : typing.Type[sm.material_data], - rhoL : float, E_in : npt.NDArray, rhoL_func : callable) -> npt.NDArray: + +def mat_transmission( + mat: typing.Type[sm.material_data], rhoL: float, E_in: npt.NDArray, rhoL_func: callable +) -> npt.NDArray: """ Calculates the straight line transmission of primary fusion sources through a material's areal density - The areal density function rhoL_func needs to be a callable function with a + The areal density function rhoL_func needs to be a callable function with a single argument (cosine[theta]) Args: mat (material_data) : Material data class for attenuating medium rhoL (float): the (4-pi averaged) areal density of the DT E_in (numpy.array): the energy array of the primary neutron spectrum - rhoL_func (callable): must be a single argument function f(x), + rhoL_func (callable): must be a single argument function f(x), where x e [-1,1] and f(x) e [0,inf] and int f(x) dx = 1 Returns: numpy.array: the total transmission coefficient as calculated - + exp[ - A_1S rhoL_func(1) (f_D sigma_{D,tot}(Ein) + f_T sigma_{T,tot}(Ein)) ] """ tot_xsec = mat.sigma_tot(E_in) A_1S = mat.rhoR_2_A1s(rhoL) - return np.exp(-A_1S*rhoL_func(1.0)*tot_xsec) + return np.exp(-A_1S * rhoL_func(1.0) * tot_xsec) + ############################### # Full model fitting function # ############################### -def calc_DT_sigmabar(Ein : npt.NDArray, I_E : npt.NDArray, - frac_D: float = frac_D_default, frac_T: float = frac_T_default) -> float: + +def calc_DT_sigmabar( + Ein: npt.NDArray, I_E: npt.NDArray, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> float: """Calculates the spectral-averaged cross section for DT Args: @@ -585,11 +612,15 @@ def calc_DT_sigmabar(Ein : npt.NDArray, I_E : npt.NDArray, Returns: float : the spectrally averaged total DT cross section """ - sigmabar = frac_D*np.trapz(mat_dict['D'].sigma_tot(Ein)*I_E,Ein)+frac_T*np.trapz(mat_dict['T'].sigma_tot(Ein)*I_E,Ein) + sigmabar = frac_D * np.trapezoid(mat_dict["D"].sigma_tot(Ein) * I_E, Ein) + frac_T * np.trapezoid( + mat_dict["T"].sigma_tot(Ein) * I_E, Ein + ) return sigmabar -def rhoR_2_A1s(rhoR : typing.Union[float,npt.NDArray], - frac_D: float = frac_D_default,frac_T: float = frac_T_default) -> typing.Union[float,npt.NDArray]: + +def rhoR_2_A1s( + rhoR: float | npt.NDArray, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> float | npt.NDArray: """Calculates the scattering amplitude given a DT areal density in kg/m^2 Args: @@ -600,13 +631,15 @@ def rhoR_2_A1s(rhoR : typing.Union[float,npt.NDArray], Returns: float : the scattering amplitude for single scattering """ - mbar = (frac_D*sm.A_D+frac_T*sm.A_T)*Mn_kg - A_1S = rhoR*(sigmabarn/mbar) + mbar = (frac_D * sm.A_D + frac_T * sm.A_T) * Mn_kg + A_1S = rhoR * (sigmabarn / mbar) return A_1S -def A1s_2_rhoR(A_1S : typing.Union[float,npt.NDArray], - frac_D: float = frac_D_default,frac_T: float = frac_T_default) -> typing.Union[float,npt.NDArray]: - """Calculates the DT areal density in kg/m^2 given a scattering amplitude + +def A1s_2_rhoR( + A_1S: float | npt.NDArray, frac_D: float = frac_D_default, frac_T: float = frac_T_default +) -> float | npt.NDArray: + """Calculates the DT areal density in kg/m^2 given a scattering amplitude Args: A_1S (float): the scattering amplitude for single scattering @@ -616,6 +649,6 @@ def A1s_2_rhoR(A_1S : typing.Union[float,npt.NDArray], Returns: float : the DT areal density in kg/m^2 """ - mbar = (frac_D*sm.A_D+frac_T*sm.A_T)*Mn_kg - rhoR = A_1S/(sigmabarn/mbar) + mbar = (frac_D * sm.A_D + frac_T * sm.A_T) * Mn_kg + rhoR = A_1S / (sigmabarn / mbar) return rhoR diff --git a/src/NeSST/cross_sections.py b/src/NeSST/cross_sections.py index 3371884..d0938b9 100644 --- a/src/NeSST/cross_sections.py +++ b/src/NeSST/cross_sections.py @@ -1,129 +1,144 @@ # Backend of spectral model +from dataclasses import dataclass + 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 +from NeSST.constants import * +from NeSST.utils import * ############################### # Differential cross sections # ############################### + @dataclass class NeSST_SDX: - Ein : npt.NDArray - points : npt.NDArray - values : npt.NDArray + 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 + 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())) +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) + interp = griddata(table.points, table.values, xi, rescale=True).reshape(mu.shape) - ans = sig*interp + 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(legendre_dx_spline,E_vec): +def interp_Tlcoeff(legendre_dx_spline, E_vec): size = [E_vec.shape[0]] NTl = len(legendre_dx_spline) size.append(NTl) Tlcoeff = np.zeros(size) for i in range(NTl): - Tlcoeff[:,i] = legendre_dx_spline[i](E_vec) - return Tlcoeff,NTl + Tlcoeff[:, i] = legendre_dx_spline[i](E_vec) + return Tlcoeff, NTl + # Evaluate the differential cross section by combining legendre and cross section # See https://t2.lanl.gov/nis/endf/intro20.html -def diffxsec_legendre_eval(sig,mu,coeff): +def diffxsec_legendre_eval(sig, mu, coeff): c = coeff.T ans = np.zeros_like(mu) - if(len(mu.shape) == 1): - ans = sig*legval(mu,c,tensor=False) - elif(len(mu.shape) == 2): - ans = sig*legval(mu,c[:,None,:],tensor=False) - elif(len(mu.shape) == 3): - ans = sig*legval(mu,c[:,None,None,:],tensor=False) + if len(mu.shape) == 1: + ans = sig * legval(mu, c, tensor=False) + elif len(mu.shape) == 2: + ans = sig * legval(mu, c[:, None, :], tensor=False) + elif len(mu.shape) == 3: + ans = sig * legval(mu, c[:, None, None, :], tensor=False) ans[np.abs(mu) > 1.0] = 0.0 return ans -# CoM frame differential cross section wrapper fucntion -def f_dsdO(Ein_vec,mu,material): +# CoM frame differential cross section wrapper fucntion +def f_dsdO(Ein_vec, mu, material): 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 + 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) + dsdO = diffxsec_legendre_eval(sig, mu, Tlcoeff_interp) return dsdO + # Differential cross section even larger wrapper function -def dsigdOmega(A,Ein,Eout,Ein_vec,muin,muout,vf,material): - mu_CoM = col.muc(A,Ein,Eout,muin,muout,vf) - return f_dsdO(Ein_vec,mu_CoM,material) +def dsigdOmega(A, Ein, Eout, Ein_vec, muin, muout, vf, material): + mu_CoM = col.muc(A, Ein, Eout, muin, muout, vf) + return f_dsdO(Ein_vec, mu_CoM, material) + # Inelastic double differential cross sections # Reads and interpolated data saved in the ENDF interpreted data format class doubledifferentialcrosssection_data: + 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 - 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 = {} 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) + 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): - + def interpolate(self, Ein, mu, Eout): # Find indices # Energies - if(Ein == np.amax(self.Ein_ddx)): + if Ein == np.amax(self.Ein_ddx): Eidx2 = self.NEin_ddx - 1 Eidx1 = Eidx2 - 1 - elif(Ein < np.amin(self.Ein_ddx)): + elif Ein < np.amin(self.Ein_ddx): return 0.0 else: Eidx2 = np.argmax(self.Ein_ddx > Ein) Eidx1 = Eidx2 - 1 # Angles - if(mu == +1.0): - Cidx12 = self.Ncos_ddx[Eidx1]-1 + if mu == +1.0: + Cidx12 = self.Ncos_ddx[Eidx1] - 1 Cidx11 = Cidx12 - 1 - Cidx22 = self.Ncos_ddx[Eidx2]-1 + Cidx22 = self.Ncos_ddx[Eidx2] - 1 Cidx21 = Cidx22 - 1 - elif(mu == -1.0): + elif mu == -1.0: Cidx12 = 1 Cidx11 = 0 Cidx22 = 1 @@ -135,79 +150,81 @@ def interpolate(self,Ein,mu,Eout): Cidx21 = Cidx22 - 1 # Find interpolation factors - mu_x1 = (mu-self.cos_ddx[Eidx1][Cidx11])/(self.cos_ddx[Eidx1][Cidx12]-self.cos_ddx[Eidx1][Cidx11]) - mu_x2 = (mu-self.cos_ddx[Eidx2][Cidx21])/(self.cos_ddx[Eidx2][Cidx22]-self.cos_ddx[Eidx2][Cidx21]) - Ein_x = (Ein-self.Ein_ddx[Eidx1])/(self.Ein_ddx[Eidx2]-self.Ein_ddx[Eidx1]) + mu_x1 = (mu - self.cos_ddx[Eidx1][Cidx11]) / (self.cos_ddx[Eidx1][Cidx12] - self.cos_ddx[Eidx1][Cidx11]) + mu_x2 = (mu - self.cos_ddx[Eidx2][Cidx21]) / (self.cos_ddx[Eidx2][Cidx22] - self.cos_ddx[Eidx2][Cidx21]) + Ein_x = (Ein - self.Ein_ddx[Eidx1]) / (self.Ein_ddx[Eidx2] - self.Ein_ddx[Eidx1]) x_112 = mu_x1 - x_111 = (1-x_112) + x_111 = 1 - x_112 x_222 = mu_x2 - x_221 = (1-x_222) + x_221 = 1 - x_222 x_2 = Ein_x - x_1 = (1-x_2) + x_1 = 1 - x_2 # Unit base transform - E_h11 = self.Emax_ddx[(Eidx1,Cidx11)] - E_h12 = self.Emax_ddx[(Eidx1,Cidx12)] - E_h21 = self.Emax_ddx[(Eidx2,Cidx21)] - E_h22 = self.Emax_ddx[(Eidx2,Cidx22)] - E_h1 = E_h11 + mu_x1*(E_h12-E_h11) - E_h2 = E_h21 + mu_x2*(E_h22-E_h21) - E_high = E_h1 + Ein_x*(E_h2-E_h1) - if(E_high == 0.0): + E_h11 = self.Emax_ddx[(Eidx1, Cidx11)] + E_h12 = self.Emax_ddx[(Eidx1, Cidx12)] + E_h21 = self.Emax_ddx[(Eidx2, Cidx21)] + E_h22 = self.Emax_ddx[(Eidx2, Cidx22)] + E_h1 = E_h11 + mu_x1 * (E_h12 - E_h11) + E_h2 = E_h21 + mu_x2 * (E_h22 - E_h21) + E_high = E_h1 + Ein_x * (E_h2 - E_h1) + if E_high == 0.0: return 0.0 - J_111 = self.Emax_ddx[(Eidx1,Cidx11)]/E_high - J_112 = self.Emax_ddx[(Eidx1,Cidx12)]/E_high - J_221 = self.Emax_ddx[(Eidx2,Cidx21)]/E_high - J_222 = self.Emax_ddx[(Eidx2,Cidx22)]/E_high + J_111 = self.Emax_ddx[(Eidx1, Cidx11)] / E_high + J_112 = self.Emax_ddx[(Eidx1, Cidx12)] / E_high + J_221 = self.Emax_ddx[(Eidx2, Cidx21)] / E_high + J_222 = self.Emax_ddx[(Eidx2, Cidx22)] / E_high # Find unit base transformed energy - Eout_111 = Eout*J_111 - Eout_112 = Eout*J_112 - Eout_221 = Eout*J_221 - Eout_222 = Eout*J_222 + Eout_111 = Eout * J_111 + Eout_112 = Eout * J_112 + Eout_221 = Eout * J_221 + Eout_222 = Eout * J_222 - f_111 = self.f_ddx_interp[(Eidx1,Cidx11)](Eout_111)*J_111 - f_112 = self.f_ddx_interp[(Eidx1,Cidx12)](Eout_112)*J_112 - f_221 = self.f_ddx_interp[(Eidx2,Cidx21)](Eout_221)*J_221 - f_222 = self.f_ddx_interp[(Eidx2,Cidx22)](Eout_222)*J_222 + f_111 = self.f_ddx_interp[(Eidx1, Cidx11)](Eout_111) * J_111 + f_112 = self.f_ddx_interp[(Eidx1, Cidx12)](Eout_112) * J_112 + f_221 = self.f_ddx_interp[(Eidx2, Cidx21)](Eout_221) * J_221 + f_222 = self.f_ddx_interp[(Eidx2, Cidx22)](Eout_222) * J_222 - f_1 = x_111*f_111+x_112*f_112 - f_2 = x_221*f_221+x_222*f_222 + f_1 = x_111 * f_111 + x_112 * f_112 + f_2 = x_221 * f_221 + x_222 * f_222 - f_ddx = x_1*f_1+x_2*f_2 + f_ddx = x_1 * f_1 + x_2 * f_2 return f_ddx - def regular_grid(self,Ein,mu,Eout): - self.rgrid_shape = (Ein.shape[0],mu.shape[0],Eout.shape[0]) + def regular_grid(self, Ein, mu, Eout): + self.rgrid_shape = (Ein.shape[0], mu.shape[0], Eout.shape[0]) self.rgrid = np.zeros(self.rgrid_shape) for i in range(Ein.shape[0]): for j in range(mu.shape[0]): - self.rgrid[i,j,:] = 2.*self.xsec_interp(Ein[i])*self.interpolate(Ein[i],mu[j],Eout) + self.rgrid[i, j, :] = 2.0 * self.xsec_interp(Ein[i]) * self.interpolate(Ein[i], mu[j], Eout) -class doubledifferentialcrosssection_LAW6: - 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 - E_a = self.A_t*Ein/(self.A_p+self.A_t)+self.Q_react - E_max = (self.A_tot-1.0)*E_a/self.A_tot - C3 = 4.0/(np.pi*E_max*E_max) - square_bracket_term = E_max-(E_star+Eout-2*mu*np.sqrt(E_star*Eout)) +class doubledifferentialcrosssection_LAW6: + 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 + E_a = self.A_t * Ein / (self.A_p + self.A_t) + self.Q_react + E_max = (self.A_tot - 1.0) * E_a / self.A_tot + C3 = 4.0 / (np.pi * E_max * E_max) + square_bracket_term = E_max - (E_star + Eout - 2 * mu * np.sqrt(E_star * Eout)) square_bracket_term[square_bracket_term < 0.0] = 0.0 - f_ddx = C3*np.sqrt(Eout*square_bracket_term) + f_ddx = C3 * np.sqrt(Eout * square_bracket_term) return f_ddx - def regular_grid(self,Ein,mu,Eout): - Ei,Mm,Eo = np.meshgrid(Ein,mu,Eout,indexing='ij') - self.rgrid = 2.*self.xsec_interp(Ei)*self.ddx(Ei,Mm,Eo) \ No newline at end of file + def regular_grid(self, Ein, mu, Eout): + Ei, Mm, Eo = np.meshgrid(Ein, mu, Eout, indexing="ij") + self.rgrid = 2.0 * self.xsec_interp(Ei) * self.ddx(Ei, Mm, Eo) diff --git a/src/NeSST/data/Be9.json b/src/NeSST/data/Be9.json index b7e0d5a..ea7f237 100644 --- a/src/NeSST/data/Be9.json +++ b/src/NeSST/data/Be9.json @@ -6,4 +6,4 @@ "inelastic" : true, "n2n" : true } -} \ No newline at end of file +} diff --git a/src/NeSST/data/C12.json b/src/NeSST/data/C12.json index 2bd848f..c2a607d 100644 --- a/src/NeSST/data/C12.json +++ b/src/NeSST/data/C12.json @@ -6,4 +6,4 @@ "inelastic" : true, "n2n" : false } -} \ No newline at end of file +} diff --git a/src/NeSST/data/H1.json b/src/NeSST/data/H1.json index 9ed8d56..b900986 100644 --- a/src/NeSST/data/H1.json +++ b/src/NeSST/data/H1.json @@ -6,4 +6,4 @@ "inelastic" : false, "n2n" : false } -} \ No newline at end of file +} diff --git a/src/NeSST/data/H2.json b/src/NeSST/data/H2.json index c44602c..9abb1ec 100644 --- a/src/NeSST/data/H2.json +++ b/src/NeSST/data/H2.json @@ -7,11 +7,11 @@ "n2n" : false } , - "n-001_H_002.cendl" : + "n-001_H_002.cendl" : { "total" : false, "elastic" : false, "inelastic" : false, "n2n" : true } -} \ No newline at end of file +} diff --git a/src/NeSST/data/H3.json b/src/NeSST/data/H3.json index 8aa6bbd..c50d4ea 100644 --- a/src/NeSST/data/H3.json +++ b/src/NeSST/data/H3.json @@ -6,11 +6,11 @@ "inelastic" : false, "n2n" : false }, - "n-001_H_003.brond" : + "n-001_H_003.brond" : { "total" : false, "elastic" : false, "inelastic" : false, "n2n" : true } -} \ No newline at end of file +} diff --git a/src/NeSST/data/TT_reac_Hale.dat b/src/NeSST/data/TT_reac_Hale.dat index 913f94d..03e1cb3 100644 --- a/src/NeSST/data/TT_reac_Hale.dat +++ b/src/NeSST/data/TT_reac_Hale.dat @@ -118,4 +118,4 @@ 7.9433e-01 9.3787e-17 8.5770e-01 9.4076e-17 9.2612e-01 9.4019e-17 -1.0000e+00 9.3646e-17 \ No newline at end of file +1.0000e+00 9.3646e-17 diff --git a/src/NeSST/endf_interface.py b/src/NeSST/endf_interface.py index 2cfa81e..bb058eb 100644 --- a/src/NeSST/endf_interface.py +++ b/src/NeSST/endf_interface.py @@ -1,10 +1,12 @@ -import endf -from NeSST.constants import * -from NeSST.cross_sections import NeSST_SDX, NeSST_DDX -import numpy as np import json from dataclasses import dataclass +import endf +import numpy as np + +from NeSST.constants import * +from NeSST.cross_sections import NeSST_DDX, NeSST_SDX + # Some defs # See https://t2.lanl.gov/nis/endf/mts.html @@ -17,189 +19,222 @@ MT_N2NXSECTION = 16 MT_INELXSECTION_START = 51 + @dataclass class ENDFManifest: - total : bool - elastic : bool - inelastic : bool - n2n : bool + total: bool + elastic: bool + inelastic: bool + n2n: bool + def retrieve_ENDF_data(json_file): - mainfest = convert_json_to_manifest(data_dir+json_file) + mainfest = convert_json_to_manifest(data_dir + json_file) ENDF_data = retrieve_ENDF_data_from_manifest(mainfest) return ENDF_data + def convert_json_to_manifest(json_file): - with open(json_file, 'r') as js: + with open(json_file, "r") as js: json_dict = json.load(js) manifest_dict = {} for filename, file_manifest in json_dict.items(): - manifest_dict[filename] = ENDFManifest(total=file_manifest['total'], - elastic=file_manifest['elastic'], - inelastic=file_manifest['inelastic'], - n2n=file_manifest['n2n']) + manifest_dict[filename] = ENDFManifest( + total=file_manifest["total"], + elastic=file_manifest["elastic"], + inelastic=file_manifest["inelastic"], + n2n=file_manifest["n2n"], + ) return manifest_dict + def retrieve_ENDF_data_from_manifest(manifest): - return_dict = {'interactions' : ENDFManifest(total=False, - elastic=False, - inelastic=False, - n2n=False)} + return_dict = {"interactions": ENDFManifest(total=False, elastic=False, inelastic=False, n2n=False)} for filename, file_manifest in manifest.items(): + mat = endf.Material(ENDF_dir + filename) - mat = endf.Material(ENDF_dir+filename) - - if(file_manifest.total == True): - if(return_dict['interactions'].total): - print(f'{filename}: Warning, loading total cross section data from multiple files...') + if file_manifest.total: + if return_dict["interactions"].total: + print(f"{filename}: Warning, loading total cross section data from multiple files...") else: - return_dict['interactions'].total = True + return_dict["interactions"].total = True - total_dataset = mat.section_data[MF_XSECTION,MT_TOTXSECTION] - return_dict['A'] = total_dataset['AWR'] - total_table = total_dataset['sigma'] + total_dataset = mat.section_data[MF_XSECTION, MT_TOTXSECTION] + return_dict["A"] = total_dataset["AWR"] + total_table = total_dataset["sigma"] total_E, total_sigma = total_table.x, total_table.y - return_dict['total_xsec'] = {'E' : total_E, 'sig' : total_sigma} + return_dict["total_xsec"] = {"E": total_E, "sig": total_sigma} - if(file_manifest.elastic == True): - if(return_dict['interactions'].elastic): - print(f'{filename}: Warning, loading elastic cross section data from multiple files...') + if file_manifest.elastic: + if return_dict["interactions"].elastic: + print(f"{filename}: Warning, loading elastic cross section data from multiple files...") else: - return_dict['interactions'].elastic = True + return_dict["interactions"].elastic = True - elastic_table = mat.section_data[MF_XSECTION,MT_ELXSECTION]['sigma'] + elastic_table = mat.section_data[MF_XSECTION, MT_ELXSECTION]["sigma"] elastic_E, elastic_sigma = elastic_table.x, elastic_table.y - return_dict['elastic_xsec'] = {'E' : elastic_E, 'sig' : elastic_sigma} - - LTT = mat.section_data[MF_ENERGYANGLEDISTS,MT_ELXSECTION]['LTT'] - - if(LTT == 0): # Isotropic - return_dict['elastic_dxsec'] = {'legendre' : True, 'E' : np.array([elastic_E[0],elastic_E[-1]]), 'a_l' : np.zeros((2,1)), 'N_l' : 1} - - elif(LTT == 1 or LTT == 3): # Legendre - if(LTT == 3): - print(f'{filename}: Warning, LTT ({LTT}) for elastic scattering, using Legendre only') - elastic_table = mat.section_data[MF_ENERGYANGLEDISTS,MT_ELXSECTION]['legendre'] - elastic_E, elastic_al = elastic_table['E'], convert_to_NeSST_legendre_format(elastic_table['a_l']) - return_dict['elastic_dxsec'] = {'legendre' : True, 'E' : elastic_E, 'a_l' : elastic_al, 'N_l' : elastic_al.shape[1]} - - elif(LTT == 2): # Tabulated - elastic_table = mat.section_data[MF_ENERGYANGLEDISTS,MT_ELXSECTION]['tabulated'] - elastic_E = elastic_table['E'] + return_dict["elastic_xsec"] = {"E": elastic_E, "sig": elastic_sigma} + + LTT = mat.section_data[MF_ENERGYANGLEDISTS, MT_ELXSECTION]["LTT"] + + if LTT == 0: # Isotropic + return_dict["elastic_dxsec"] = { + "legendre": True, + "E": np.array([elastic_E[0], elastic_E[-1]]), + "a_l": np.zeros((2, 1)), + "N_l": 1, + } + + elif LTT == 1 or LTT == 3: # Legendre + if LTT == 3: + print(f"{filename}: Warning, LTT ({LTT}) for elastic scattering, using Legendre only") + elastic_table = mat.section_data[MF_ENERGYANGLEDISTS, MT_ELXSECTION]["legendre"] + elastic_E, elastic_al = elastic_table["E"], convert_to_NeSST_legendre_format(elastic_table["a_l"]) + return_dict["elastic_dxsec"] = { + "legendre": True, + "E": elastic_E, + "a_l": elastic_al, + "N_l": elastic_al.shape[1], + } + + elif LTT == 2: # Tabulated + elastic_table = mat.section_data[MF_ENERGYANGLEDISTS, MT_ELXSECTION]["tabulated"] + elastic_E = elastic_table["E"] NE = elastic_E.shape[0] - elastic_Ein, elastic_mu, elastic_f = [] , [] , [] + elastic_Ein, elastic_mu, elastic_f = [], [], [] for iE in range(NE): - elastic_mu.append(elastic_table['mu'][iE].x) - elastic_f.append(elastic_table['mu'][iE].y) - Nmu = elastic_table['mu'][iE].x.shape[0] - elastic_Ein.append(np.repeat(elastic_E[iE],Nmu)) + elastic_mu.append(elastic_table["mu"][iE].x) + elastic_f.append(elastic_table["mu"][iE].y) + Nmu = elastic_table["mu"][iE].x.shape[0] + elastic_Ein.append(np.repeat(elastic_E[iE], Nmu)) - points = np.column_stack((np.concatenate(elastic_Ein),np.concatenate(elastic_mu))) + points = np.column_stack((np.concatenate(elastic_Ein), np.concatenate(elastic_mu))) values = np.concatenate(elastic_f) - SDX = NeSST_SDX(Ein=elastic_E,points=points,values=values) - return_dict['elastic_dxsec'] = {'legendre' : False, 'SDX' : SDX} + SDX = NeSST_SDX(Ein=elastic_E, points=points, values=values) + return_dict["elastic_dxsec"] = {"legendre": False, "SDX": SDX} else: - print(f'{filename}: Warning, LTT ({LTT}) for elastic scattering not recognised...') + print(f"{filename}: Warning, LTT ({LTT}) for elastic scattering not recognised...") - if(file_manifest.inelastic == True): - if(return_dict['interactions'].inelastic): - print(f'{filename}: Warning, loading inelastic cross section data from multiple files...') + if file_manifest.inelastic: + if return_dict["interactions"].inelastic: + print(f"{filename}: Warning, loading inelastic cross section data from multiple files...") else: - return_dict['interactions'].inelastic = True + return_dict["interactions"].inelastic = True i_inelastic = 1 while True: - MT = MT_INELXSECTION_START+(i_inelastic-1) + MT = MT_INELXSECTION_START + (i_inelastic - 1) try: - inelastic_table = mat.section_data[MF_XSECTION,MT]['sigma'] + inelastic_table = mat.section_data[MF_XSECTION, MT]["sigma"] inelastic_E, inelastic_sigma = inelastic_table.x, inelastic_table.y - return_dict[f'inelastic_xsec_n{i_inelastic}'] = {'E' : inelastic_E, 'sig' : inelastic_sigma} - inelastic_Q = mat.section_data[MF_XSECTION,MT]['QI'] - except: + return_dict[f"inelastic_xsec_n{i_inelastic}"] = {"E": inelastic_E, "sig": inelastic_sigma} + inelastic_Q = mat.section_data[MF_XSECTION, MT]["QI"] + except: # noqa: E722 break - LTT = mat.section_data[MF_ENERGYANGLEDISTS,MT]['LTT'] - - if(LTT == 0): # Isotropic - return_dict[f'inelastic_dxsec_n{i_inelastic}'] = {'legendre' : True, 'E' : np.array([inelastic_E[0],inelastic_E[-1]]), 'a_l' : np.zeros((2,1)), 'N_l' : 1, 'Q' : inelastic_Q} - - elif(LTT == 1 or LTT == 3): # Legendre - if(LTT == 3): - print(f'{filename}: Warning, LTT ({LTT}) for inelastic scattering, level {i_inelastic}, using Legendre only') - inelastic_table = mat.section_data[MF_ENERGYANGLEDISTS,MT]['legendre'] - inelastic_E, inelastic_al = inelastic_table['E'], convert_to_NeSST_legendre_format(inelastic_table['a_l']) - return_dict[f'inelastic_dxsec_n{i_inelastic}'] = {'legendre' : True, 'E' : inelastic_E, 'a_l' : inelastic_al, 'N_l' : inelastic_al.shape[1], 'Q' : inelastic_Q} - - elif(LTT == 2): # Tabulated - inelastic_table = mat.section_data[MF_ENERGYANGLEDISTS,MT]['tabulated'] - inelastic_E = inelastic_table['E'] + LTT = mat.section_data[MF_ENERGYANGLEDISTS, MT]["LTT"] + + if LTT == 0: # Isotropic + return_dict[f"inelastic_dxsec_n{i_inelastic}"] = { + "legendre": True, + "E": np.array([inelastic_E[0], inelastic_E[-1]]), + "a_l": np.zeros((2, 1)), + "N_l": 1, + "Q": inelastic_Q, + } + + elif LTT == 1 or LTT == 3: # Legendre + if LTT == 3: + print( + f"{filename}: Warning, LTT ({LTT}) for inelastic scattering, level {i_inelastic}, using Legendre only" + ) + inelastic_table = mat.section_data[MF_ENERGYANGLEDISTS, MT]["legendre"] + inelastic_E, inelastic_al = ( + inelastic_table["E"], + convert_to_NeSST_legendre_format(inelastic_table["a_l"]), + ) + return_dict[f"inelastic_dxsec_n{i_inelastic}"] = { + "legendre": True, + "E": inelastic_E, + "a_l": inelastic_al, + "N_l": inelastic_al.shape[1], + "Q": inelastic_Q, + } + + elif LTT == 2: # Tabulated + inelastic_table = mat.section_data[MF_ENERGYANGLEDISTS, MT]["tabulated"] + inelastic_E = inelastic_table["E"] NE = inelastic_E.shape[0] - inelastic_mu, Nmu, inelastic_f = [] , [] , [] + inelastic_mu, Nmu, inelastic_f = [], [], [] for iE in range(NE): - inelastic_mu.append(inelastic_table['mu'][iE].x) - inelastic_f.append(inelastic_table['mu'][iE].y) - Nmu.append(inelastic_table['mu'][iE].x.shape[0]) - SDX = NeSST_SDX(NEin=NE,Ein=inelastic_E,Ncos=Nmu,cos=inelastic_mu,f=inelastic_f) - return_dict[f'inelastic_dxsec_n{i_inelastic}'] = {'legendre' : False, 'SDX' : SDX, 'Q' : inelastic_Q} + inelastic_mu.append(inelastic_table["mu"][iE].x) + inelastic_f.append(inelastic_table["mu"][iE].y) + Nmu.append(inelastic_table["mu"][iE].x.shape[0]) + SDX = NeSST_SDX(NEin=NE, Ein=inelastic_E, Ncos=Nmu, cos=inelastic_mu, f=inelastic_f) + return_dict[f"inelastic_dxsec_n{i_inelastic}"] = {"legendre": False, "SDX": SDX, "Q": inelastic_Q} else: - print(f'{filename}: Warning, LTT ({LTT}) for inelastic scattering, level {i_inelastic}, not recognised...') + print( + f"{filename}: Warning, LTT ({LTT}) for inelastic scattering, level {i_inelastic}, not recognised..." + ) i_inelastic += 1 - return_dict['n_inelastic'] = i_inelastic-1 + return_dict["n_inelastic"] = i_inelastic - 1 - if(file_manifest.n2n == True): - if(return_dict['interactions'].n2n): - print(f'{filename}: Warning, loading n,2n cross section data from multiple files...') + if file_manifest.n2n: + if return_dict["interactions"].n2n: + print(f"{filename}: Warning, loading n,2n cross section data from multiple files...") else: - return_dict['interactions'].n2n = True + return_dict["interactions"].n2n = True - n2n_table = mat.section_data[MF_XSECTION,MT_N2NXSECTION]['sigma'] + n2n_table = mat.section_data[MF_XSECTION, MT_N2NXSECTION]["sigma"] n2n_E, n2n_sigma = n2n_table.x, n2n_table.y - return_dict['n2n_xsec'] = {'E' : n2n_E, 'sig' : n2n_sigma} + return_dict["n2n_xsec"] = {"E": n2n_E, "sig": n2n_sigma} - n2n_Q = mat.section_data[MF_XSECTION,MT_N2NXSECTION]['QI'] - n2n_table = mat.section_data[MF_PRODENERGYANGLEDISTS,MT_N2NXSECTION] - n2n_products = n2n_table['products'] + n2n_Q = mat.section_data[MF_XSECTION, MT_N2NXSECTION]["QI"] + n2n_table = mat.section_data[MF_PRODENERGYANGLEDISTS, MT_N2NXSECTION] + n2n_products = n2n_table["products"] for n2n_product in n2n_products: - if(n2n_product['AWP'] == 1.0): # is neutron? - if(n2n_product['LAW'] == 6): - return_dict['n2n_dxsec'] = {'LAW' : 6, - 'A_i' : n2n_product['AWP'], - 'A_e' : n2n_product['AWP'], - 'A_t' : n2n_table['AWR'], - 'A_p' : n2n_product['AWP'], - 'A_tot' : n2n_product['distribution']['APSX'], - 'Q_react' : n2n_Q} - elif(n2n_product['LAW'] == 7): - DDX = convert_to_NeSST_LAW7_format(n2n_product['distribution']) - return_dict['n2n_dxsec'] = {'LAW' : 7, - 'DDX' : DDX, - 'Q_react' : n2n_Q} + if n2n_product["AWP"] == 1.0: # is neutron? + if n2n_product["LAW"] == 6: + return_dict["n2n_dxsec"] = { + "LAW": 6, + "A_i": n2n_product["AWP"], + "A_e": n2n_product["AWP"], + "A_t": n2n_table["AWR"], + "A_p": n2n_product["AWP"], + "A_tot": n2n_product["distribution"]["APSX"], + "Q_react": n2n_Q, + } + elif n2n_product["LAW"] == 7: + DDX = convert_to_NeSST_LAW7_format(n2n_product["distribution"]) + return_dict["n2n_dxsec"] = {"LAW": 7, "DDX": DDX, "Q_react": n2n_Q} else: print(f"{filename}: n2n data LAW = {n2n_product['LAW']}, NeSST cannot use...") return return_dict -def retrieve_total_cross_section_from_ENDF_file(filename): - mat = endf.Material(ENDF_dir+filename) +def retrieve_total_cross_section_from_ENDF_file(filename): + mat = endf.Material(ENDF_dir + filename) - total_dataset = mat.section_data[MF_XSECTION,MT_TOTXSECTION] - total_table = total_dataset['sigma'] + total_dataset = mat.section_data[MF_XSECTION, MT_TOTXSECTION] + total_table = total_dataset["sigma"] total_E, total_sigma = total_table.x, total_table.y - return total_E, total_sigma + return total_E, total_sigma + def convert_to_NeSST_legendre_format(ENDF_al): max_Nl = max([len(a) for a in ENDF_al]) - uniform_len_al = [np.pad(a,(0,max_Nl-len(a))) for a in ENDF_al] + uniform_len_al = [np.pad(a, (0, max_Nl - len(a))) for a in ENDF_al] return np.array(uniform_len_al) + def convert_to_NeSST_LAW7_format(ENDF_dist): - NE = ENDF_dist['NE'] + NE = ENDF_dist["NE"] Ein = [] Ncos = [] cos = [] @@ -208,29 +243,20 @@ def convert_to_NeSST_LAW7_format(ENDF_dist): f = {} Emax = {} for iE in range(NE): - dist_Ein = ENDF_dist['distribution'][iE] - Ein.append(dist_Ein['E']) - NMU = dist_Ein['NMU'] + dist_Ein = ENDF_dist["distribution"][iE] + Ein.append(dist_Ein["E"]) + NMU = dist_Ein["NMU"] Ncos.append(NMU) cos_arr = np.zeros(NMU) for iM in range(NMU): - dist_mu = dist_Ein['mu'][iM] - cos_arr[iM] = dist_mu['mu'] - NEout[(iE,iM)] = dist_mu['f'].x.shape[0] - Eout[(iE,iM)] = dist_mu['f'].x - f[(iE,iM)] = dist_mu['f'].y - Emax[(iE,iM)] = np.amax(dist_mu['f'].x) + dist_mu = dist_Ein["mu"][iM] + cos_arr[iM] = dist_mu["mu"] + NEout[(iE, iM)] = dist_mu["f"].x.shape[0] + Eout[(iE, iM)] = dist_mu["f"].x + f[(iE, iM)] = dist_mu["f"].y + Emax[(iE, iM)] = np.amax(dist_mu["f"].x) cos.append(cos_arr) - DDX = NeSST_DDX( - NEin = NE, - Ein = Ein, - Ncos = Ncos, - cos = cos, - NEout = NEout, - Eout = Eout, - f = f, - Emax = Emax - ) - - return DDX \ No newline at end of file + DDX = NeSST_DDX(NEin=NE, Ein=Ein, Ncos=Ncos, cos=cos, NEout=NEout, Eout=Eout, f=f, Emax=Emax) + + return DDX diff --git a/src/NeSST/fitting.py b/src/NeSST/fitting.py index 9b05218..1f08967 100644 --- a/src/NeSST/fitting.py +++ b/src/NeSST/fitting.py @@ -1,164 +1,200 @@ -from NeSST.core import * from NeSST.constants import * +from NeSST.core import * from NeSST.utils import * + class DT_fit_function: - """ - A class which constructs various simple models for the full spectrum in DT for fitting data - - User must provide energy grids (and velocity grids if including ion kinematics) - - One can create a model function from the following list of approximations: - - -- Symmetric areal density - -- Asymmetric Mode 1 areal density - - The primary spectra are assumed isotropic and with moments defined by a single temperature - - This class doesn't represent the full set of spectrum models which can be produced by NeSST! Just a common few... - """ - def __init__(self,E_DTspec,E_sspec,vion_arr = None): - self.E_DTspec = E_DTspec - self.E_sspec = E_sspec - print("### Initialising data on energy grids... ###") - init_DT_scatter(E_sspec,E_DTspec) - if(vion_arr is not None): - self.ion_kinematics = True - self.vion_arr = vion_arr - print("### Initialising scattering matrices on ion velocity grid... ###") - init_DT_ionkin_scatter(vion_arr,nT=True,nD=True) - else: - self.ion_kinematics = False - print("### Init Done. ###") - - def set_primary_Tion(self,Tion): - if(Tion < 0.1): - print("~~ WARNING Low Tion (< 100 eV) ~~") - self.Tion = Tion - - self.DTmean,_,self.DTvar = DTprimspecmoments(Tion) - self.DDmean,_,self.DDvar = DDprimspecmoments(Tion) - - Y_DT = 1.0 - Y_DD = yield_from_dt_yield_ratio('dd',Y_DT,Tion) - Y_TT = yield_from_dt_yield_ratio('tt',Y_DT,Tion) - - self.dNdE_DT = Y_DT*QBrysk(self.E_DTspec,self.DTmean,self.DTvar) # Brysk shape i.e. Gaussian - self.dNdE_DD = Y_DD*QBrysk(self.E_sspec ,self.DDmean,self.DDvar) # Brysk shape i.e. Gaussian - self.dNdE_TT = Y_TT*dNdE_TT(self.E_sspec,Tion) - - self.I_DT = interpolate_1d(self.E_DTspec,self.dNdE_DT,fill_value=0.0,bounds_error=False) - self.I_DD = interpolate_1d(self.E_sspec,self.dNdE_DD) - self.I_TT = interpolate_1d(self.E_sspec,self.dNdE_TT) - - def init_symmetric_model(self): - """ - Creates a callable model function for a symmetric areal density distribution - """ - - rhoL_func = lambda x : np.ones_like(x) - - if(self.ion_kinematics): - calc_DT_ionkin_primspec_rhoL_integral(self.dNdE_DT,nT=True,nD=True) - mat_dict['D'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func) - mat_dict['T'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func) - else: - mat_dict['D'].calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func) - mat_dict['T'].calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func) - mat_dict['D'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func) - mat_dict['T'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func) - - dNdE_Dn2n = interpolate_1d(self.E_sspec,mat_dict['D'].n2n_dNdE,fill_value=0.0,bounds_error=False) - dNdE_Tn2n = interpolate_1d(self.E_sspec,mat_dict['T'].n2n_dNdE,fill_value=0.0,bounds_error=False) - - if(self.ion_kinematics): - def model(E,rhoL,vbar,dv,fT,fD,Yn): - """ - Symmetric areal density model with scattering ion velocity distribution with mean and std dev, vbar and dv in m/s - """ - A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT) - dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(E,vbar,dv) - dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(E,vbar,dv) - dNdE_tot = A_1S*(fT*dNdE_nT+fD*dNdE_nD+fD*dNdE_Dn2n(E)+fT*dNdE_Tn2n(E)) - return Yn*(dNdE_tot+(fD/fT)*(frac_T_default/frac_D_default)*self.I_DD(E)+(fT/fD)*(frac_D_default/frac_T_default)*self.I_TT(E)) - else: - """ Incomplete """ - def model(E,rhoL,Ts,fT,fD,Yn): - """ - Symmetric areal density model with scattering temperature Ts, in keV - """ - A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT) - dNdE_nT = mat_dict['T'].elastic_dNdE.copy() - dNdE_nD = mat_dict['D'].elastic_dNdE.copy() - if(Ts > 0.1): - T_MeV = Ts/1e3 - E_nT0 = ((sm.A_T-1.0)/(sm.A_T+1.0))**2*self.DTmean - dE_nT = np.sqrt(8.0*sm.A_T*E_nT0/(sm.A_T+1.0)**2*T_MeV) - E_nD0 = ((sm.A_D-1.0)/(sm.A_D+1.0))**2*self.DTmean - dE_nD = np.sqrt(8.0*sm.A_D*E_nD0/(sm.A_D+1.0)**2*T_MeV) - - - dNdE_nT = interpolate_1d(self.E_sspec,dNdE_nT,fill_value=0.0,bounds_error=False) - dNdE_nD = interpolate_1d(self.E_sspec,dNdE_nD,fill_value=0.0,bounds_error=False) - - dNdE_tot = A_1S*(fT*dNdE_nT(E)+fD*dNdE_nD(E)+fD*dNdE_Dn2n(E)+fT*dNdE_Tn2n(E)) - return Yn*(dNdE_tot+(fD/fT)*(frac_T_default/frac_D_default)*self.I_DD(E)+(fT/fD)*(frac_D_default/frac_T_default)*self.I_TT(E)) - - self.model = model - - def init_modeone_model(self,P1_arr): - """ - Creates a callable model function for a mode 1 asymmetric areal density distribution - """ - - self.P1_arr = P1_arr - - # T(n,2n) - mat_dict['T'].n2n_ddx.rgrid_IE = np.trapz(mat_dict['T'].n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0) - mat_dict['T'].n2n_dNdE_mode1 = np.trapz(mat_dict['T'].n2n_ddx.rgrid_IE[:,:,None]* - (1.0+self.P1_arr[None,None,:]*mat_dict['T'].n2n_mu[:,None,None]) - ,mat_dict['T'].n2n_mu,axis=0) - - mat_dict['T'].n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,mat_dict['T'].n2n_dNdE_mode1,bounds_error=False) - - # D(n,2n) - mat_dict['D'].n2n_ddx.rgrid_IE = np.trapz(mat_dict['D'].n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0) - mat_dict['D'].n2n_dNdE_mode1 = np.trapz(mat_dict['D'].n2n_ddx.rgrid_IE[:,:,None]* - (1.0+self.P1_arr[None,None,:]*mat_dict['D'].n2n_mu[:,None,None]) - ,mat_dict['D'].n2n_mu,axis=0) - - mat_dict['D'].n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,mat_dict['D'].n2n_dNdE_mode1,bounds_error=False) - - # nT - M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*mat_dict['T'].full_scattering_mu[:,:,:,None]) - *mat_dict['T'].full_scattering_M[:,:,:,None] - *self.dNdE_DT[None,None,:,None],self.E_DTspec,axis=2) - - mat_dict['T'].M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False) - - # nD - M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*mat_dict['D'].full_scattering_mu[:,:,:,None]) - *mat_dict['D'].full_scattering_M[:,:,:,None] - *self.dNdE_DT[None,None,:,None],self.E_DTspec,axis=2) - - mat_dict['D'].M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False) - - if(self.ion_kinematics): - def model(E,rhoL,P1,vbar,dv,fT,fD,Yn): - A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT) - mat_dict['T'].M_prim = mat_dict['T'].M_mode1_interp(P1) - mat_dict['D'].M_prim = mat_dict['D'].M_mode1_interp(P1) - dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(E,vbar,dv) - dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(E,vbar,dv) - dNdE_Tn2n = mat_dict['T'].n2n_dNdE_mode1(E,P1) - dNdE_Dn2n = mat_dict['D'].n2n_dNdE_mode1(E,P1) - dNdE_tot = A_1S*(fT*dNdE_nT+fD*dNdE_nD+fD*dNdE_Dn2n+fT*dNdE_Tn2n) - # Primary - dNdE_DD = (fD/fT)*(frac_T_default/frac_D_default)*self.I_DD(E) - dNdE_TT = (fT/fD)*(frac_D_default/frac_T_default)*self.I_TT(E) - return Yn*(dNdE_tot+dNdE_DD+dNdE_TT) - else: - """ Incomplete """ - def model(): - return None - self.model = model \ No newline at end of file + """ + A class which constructs various simple models for the full spectrum in DT for fitting data + + User must provide energy grids (and velocity grids if including ion kinematics) + + One can create a model function from the following list of approximations: + + -- Symmetric areal density + -- Asymmetric Mode 1 areal density + + The primary spectra are assumed isotropic and with moments defined by a single temperature + + This class doesn't represent the full set of spectrum models which can be produced by NeSST! Just a common few... + """ + + def __init__(self, E_DTspec, E_sspec, vion_arr=None): + self.E_DTspec = E_DTspec + self.E_sspec = E_sspec + print("### Initialising data on energy grids... ###") + init_DT_scatter(E_sspec, E_DTspec) + if vion_arr is not None: + self.ion_kinematics = True + self.vion_arr = vion_arr + print("### Initialising scattering matrices on ion velocity grid... ###") + init_DT_ionkin_scatter(vion_arr, nT=True, nD=True) + else: + self.ion_kinematics = False + print("### Init Done. ###") + + def set_primary_Tion(self, Tion): + if Tion < 0.1: + print("~~ WARNING Low Tion (< 100 eV) ~~") + self.Tion = Tion + + self.DTmean, _, self.DTvar = DTprimspecmoments(Tion) + self.DDmean, _, self.DDvar = DDprimspecmoments(Tion) + + Y_DT = 1.0 + Y_DD = yield_from_dt_yield_ratio("dd", Y_DT, Tion) + Y_TT = yield_from_dt_yield_ratio("tt", Y_DT, Tion) + + self.dNdE_DT = Y_DT * QBrysk(self.E_DTspec, self.DTmean, self.DTvar) # Brysk shape i.e. Gaussian + self.dNdE_DD = Y_DD * QBrysk(self.E_sspec, self.DDmean, self.DDvar) # Brysk shape i.e. Gaussian + self.dNdE_TT = Y_TT * dNdE_TT(self.E_sspec, Tion) + + self.I_DT = interpolate_1d(self.E_DTspec, self.dNdE_DT, fill_value=0.0, bounds_error=False) + self.I_DD = interpolate_1d(self.E_sspec, self.dNdE_DD) + self.I_TT = interpolate_1d(self.E_sspec, self.dNdE_TT) + + def init_symmetric_model(self): + """ + Creates a callable model function for a symmetric areal density distribution + """ + + rhoL_func = lambda x: np.ones_like(x) + + if self.ion_kinematics: + calc_DT_ionkin_primspec_rhoL_integral(self.dNdE_DT, nT=True, nD=True) + mat_dict["D"].calc_n2n_dNdE(self.dNdE_DT, rhoL_func) + mat_dict["T"].calc_n2n_dNdE(self.dNdE_DT, rhoL_func) + else: + mat_dict["D"].calc_station_elastic_dNdE(self.dNdE_DT, rhoL_func) + mat_dict["T"].calc_station_elastic_dNdE(self.dNdE_DT, rhoL_func) + mat_dict["D"].calc_n2n_dNdE(self.dNdE_DT, rhoL_func) + mat_dict["T"].calc_n2n_dNdE(self.dNdE_DT, rhoL_func) + + dNdE_Dn2n = interpolate_1d(self.E_sspec, mat_dict["D"].n2n_dNdE, fill_value=0.0, bounds_error=False) + dNdE_Tn2n = interpolate_1d(self.E_sspec, mat_dict["T"].n2n_dNdE, fill_value=0.0, bounds_error=False) + + if self.ion_kinematics: + + def model(E, rhoL, vbar, dv, fT, fD, Yn): + """ + Symmetric areal density model with scattering ion velocity distribution with mean and std dev, vbar and dv in m/s + """ + A_1S = rhoR_2_A1s(rhoL, frac_D=fD, frac_T=fT) + dNdE_nT = mat_dict["T"].matrix_interpolate_gaussian(E, vbar, dv) + dNdE_nD = mat_dict["D"].matrix_interpolate_gaussian(E, vbar, dv) + dNdE_tot = A_1S * (fT * dNdE_nT + fD * dNdE_nD + fD * dNdE_Dn2n(E) + fT * dNdE_Tn2n(E)) + return Yn * ( + dNdE_tot + + (fD / fT) * (frac_T_default / frac_D_default) * self.I_DD(E) + + (fT / fD) * (frac_D_default / frac_T_default) * self.I_TT(E) + ) + else: + """ Incomplete """ + + def model(E, rhoL, Ts, fT, fD, Yn): + """ + Symmetric areal density model with scattering temperature Ts, in keV + """ + A_1S = rhoR_2_A1s(rhoL, frac_D=fD, frac_T=fT) + dNdE_nT = mat_dict["T"].elastic_dNdE.copy() + dNdE_nD = mat_dict["D"].elastic_dNdE.copy() + if Ts > 0.1: + T_MeV = Ts / 1e3 + E_nT0 = ((sm.A_T - 1.0) / (sm.A_T + 1.0)) ** 2 * self.DTmean + dE_nT = np.sqrt(8.0 * sm.A_T * E_nT0 / (sm.A_T + 1.0) ** 2 * T_MeV) # noqa + E_nD0 = ((sm.A_D - 1.0) / (sm.A_D + 1.0)) ** 2 * self.DTmean + dE_nD = np.sqrt(8.0 * sm.A_D * E_nD0 / (sm.A_D + 1.0) ** 2 * T_MeV) # noqa + + dNdE_nT = interpolate_1d(self.E_sspec, dNdE_nT, fill_value=0.0, bounds_error=False) + dNdE_nD = interpolate_1d(self.E_sspec, dNdE_nD, fill_value=0.0, bounds_error=False) + + dNdE_tot = A_1S * (fT * dNdE_nT(E) + fD * dNdE_nD(E) + fD * dNdE_Dn2n(E) + fT * dNdE_Tn2n(E)) + return Yn * ( + dNdE_tot + + (fD / fT) * (frac_T_default / frac_D_default) * self.I_DD(E) + + (fT / fD) * (frac_D_default / frac_T_default) * self.I_TT(E) + ) + + self.model = model + + def init_modeone_model(self, P1_arr): + """ + Creates a callable model function for a mode 1 asymmetric areal density distribution + """ + + self.P1_arr = P1_arr + + # T(n,2n) + mat_dict["T"].n2n_ddx.rgrid_IE = np.trapezoid( + mat_dict["T"].n2n_ddx.rgrid * self.dNdE_DT[:, None, None], self.E_DTspec, axis=0 + ) + mat_dict["T"].n2n_dNdE_mode1 = np.trapezoid( + mat_dict["T"].n2n_ddx.rgrid_IE[:, :, None] + * (1.0 + self.P1_arr[None, None, :] * mat_dict["T"].n2n_mu[:, None, None]), + mat_dict["T"].n2n_mu, + axis=0, + ) + + mat_dict["T"].n2n_dNdE_mode1 = interpolate_2d( + self.E_sspec, self.P1_arr, mat_dict["T"].n2n_dNdE_mode1, bounds_error=False + ) + + # D(n,2n) + mat_dict["D"].n2n_ddx.rgrid_IE = np.trapezoid( + mat_dict["D"].n2n_ddx.rgrid * self.dNdE_DT[:, None, None], self.E_DTspec, axis=0 + ) + mat_dict["D"].n2n_dNdE_mode1 = np.trapezoid( + mat_dict["D"].n2n_ddx.rgrid_IE[:, :, None] + * (1.0 + self.P1_arr[None, None, :] * mat_dict["D"].n2n_mu[:, None, None]), + mat_dict["D"].n2n_mu, + axis=0, + ) + + mat_dict["D"].n2n_dNdE_mode1 = interpolate_2d( + self.E_sspec, self.P1_arr, mat_dict["D"].n2n_dNdE_mode1, bounds_error=False + ) + + # nT + M_mode1 = np.trapezoid( + (1.0 + self.P1_arr[None, None, None, :] * mat_dict["T"].full_scattering_mu[:, :, :, None]) + * mat_dict["T"].full_scattering_M[:, :, :, None] + * self.dNdE_DT[None, None, :, None], + self.E_DTspec, + axis=2, + ) + + mat_dict["T"].M_mode1_interp = interpolate_1d(self.P1_arr, M_mode1, axis=-1, bounds_error=False) + + # nD + M_mode1 = np.trapezoid( + (1.0 + self.P1_arr[None, None, None, :] * mat_dict["D"].full_scattering_mu[:, :, :, None]) + * mat_dict["D"].full_scattering_M[:, :, :, None] + * self.dNdE_DT[None, None, :, None], + self.E_DTspec, + axis=2, + ) + + mat_dict["D"].M_mode1_interp = interpolate_1d(self.P1_arr, M_mode1, axis=-1, bounds_error=False) + + if self.ion_kinematics: + + def model(E, rhoL, P1, vbar, dv, fT, fD, Yn): + A_1S = rhoR_2_A1s(rhoL, frac_D=fD, frac_T=fT) + mat_dict["T"].M_prim = mat_dict["T"].M_mode1_interp(P1) + mat_dict["D"].M_prim = mat_dict["D"].M_mode1_interp(P1) + dNdE_nT = mat_dict["T"].matrix_interpolate_gaussian(E, vbar, dv) + dNdE_nD = mat_dict["D"].matrix_interpolate_gaussian(E, vbar, dv) + dNdE_Tn2n = mat_dict["T"].n2n_dNdE_mode1(E, P1) + dNdE_Dn2n = mat_dict["D"].n2n_dNdE_mode1(E, P1) + dNdE_tot = A_1S * (fT * dNdE_nT + fD * dNdE_nD + fD * dNdE_Dn2n + fT * dNdE_Tn2n) + # Primary + dNdE_DD = (fD / fT) * (frac_T_default / frac_D_default) * self.I_DD(E) + dNdE_TT = (fT / fD) * (frac_D_default / frac_T_default) * self.I_TT(E) + return Yn * (dNdE_tot + dNdE_DD + dNdE_TT) + else: + """ Incomplete """ + + def model(): + return None + + self.model = model diff --git a/src/NeSST/spectral_model.py b/src/NeSST/spectral_model.py index beaff08..822c07f 100644 --- a/src/NeSST/spectral_model.py +++ b/src/NeSST/spectral_model.py @@ -2,60 +2,81 @@ import numpy as np -from NeSST.constants import * -from NeSST.utils import * -from NeSST.endf_interface import retrieve_ENDF_data import NeSST.collisions as col import NeSST.cross_sections as xs +from NeSST.constants import * +from NeSST.endf_interface import retrieve_ENDF_data +from NeSST.utils import * ################## # Material class # ################## # A values needed for scattering kinematics -A_H = Mp/Mn -A_D = Md/Mn -A_T = Mt/Mn -A_C = MC/Mn -A_Be = MBe/Mn +A_H = Mp / Mn +A_D = Md / Mn +A_T = Mt / Mn +A_C = MC / Mn +A_Be = MBe / Mn + def unity(x): return np.ones_like(x) -class material_data: - def __init__(self,label,json): +class material_data: + def __init__(self, label, json): self.label = label self.json = json ENDF_data = retrieve_ENDF_data(self.json) - self.A = ENDF_data['A'] - - if(ENDF_data['interactions'].total): - self.sigma_tot = interpolate_1d(ENDF_data['total_xsec']['E'],ENDF_data['total_xsec']['sig'],method='linear',bounds_error=False,fill_value=0.0) - - if(ENDF_data['interactions'].elastic): - self.sigma = interpolate_1d(ENDF_data['elastic_xsec']['E'],ENDF_data['elastic_xsec']['sig'],method='linear',bounds_error=False,fill_value=0.0) - - self.elastic_legendre = ENDF_data['elastic_dxsec']['legendre'] - if(self.elastic_legendre): + self.A = ENDF_data["A"] + + if ENDF_data["interactions"].total: + self.sigma_tot = interpolate_1d( + ENDF_data["total_xsec"]["E"], + ENDF_data["total_xsec"]["sig"], + method="linear", + bounds_error=False, + fill_value=0.0, + ) + + if ENDF_data["interactions"].elastic: + self.sigma = interpolate_1d( + ENDF_data["elastic_xsec"]["E"], + ENDF_data["elastic_xsec"]["sig"], + method="linear", + bounds_error=False, + fill_value=0.0, + ) + + self.elastic_legendre = ENDF_data["elastic_dxsec"]["legendre"] + if self.elastic_legendre: self.legendre_dx_spline = [unity] - for i in range(ENDF_data['elastic_dxsec']['N_l']): - self.legendre_dx_spline.append(interpolate_1d(ENDF_data['elastic_dxsec']['E'],ENDF_data['elastic_dxsec']['a_l'][:,i],method='linear',bounds_error=False,fill_value=0.0)) + for i in range(ENDF_data["elastic_dxsec"]["N_l"]): + self.legendre_dx_spline.append( + interpolate_1d( + ENDF_data["elastic_dxsec"]["E"], + ENDF_data["elastic_dxsec"]["a_l"][:, i], + method="linear", + bounds_error=False, + fill_value=0.0, + ) + ) else: - self.elastic_SDX_table = ENDF_data['elastic_dxsec']['SDX'] + self.elastic_SDX_table = ENDF_data["elastic_dxsec"]["SDX"] - self.l_n2n = ENDF_data['interactions'].n2n - if(ENDF_data['interactions'].n2n): - if(ENDF_data['n2n_dxsec']['LAW'] == 6): - self.n2n_ddx = xs.doubledifferentialcrosssection_LAW6(ENDF_data['n2n_xsec'],ENDF_data['n2n_dxsec']) - elif(ENDF_data['n2n_dxsec']['LAW'] == 7): - self.n2n_ddx = xs.doubledifferentialcrosssection_data(ENDF_data['n2n_xsec'],ENDF_data['n2n_dxsec']) + self.l_n2n = ENDF_data["interactions"].n2n + if ENDF_data["interactions"].n2n: + if ENDF_data["n2n_dxsec"]["LAW"] == 6: + self.n2n_ddx = xs.doubledifferentialcrosssection_LAW6(ENDF_data["n2n_xsec"], ENDF_data["n2n_dxsec"]) + elif ENDF_data["n2n_dxsec"]["LAW"] == 7: + self.n2n_ddx = xs.doubledifferentialcrosssection_data(ENDF_data["n2n_xsec"], ENDF_data["n2n_dxsec"]) - self.l_inelastic = ENDF_data['interactions'].inelastic - if(ENDF_data['interactions'].inelastic): - self.n_inelastic = ENDF_data['n_inelastic'] + self.l_inelastic = ENDF_data["interactions"].inelastic + if ENDF_data["interactions"].inelastic: + self.n_inelastic = ENDF_data["n_inelastic"] self.isigma = [] self.inelasticQ = [] @@ -64,24 +85,36 @@ def __init__(self,label,json): self.inelastic_SDX_table = [] for i_inelastic in range(self.n_inelastic): - xsec_table = ENDF_data[f'inelastic_xsec_n{i_inelastic+1}'] - self.isigma.append(interpolate_1d(xsec_table['E'],xsec_table['sig'],method='linear',bounds_error=False,fill_value=0.0)) - - dxsec_table = ENDF_data[f'inelastic_dxsec_n{i_inelastic+1}'] - self.inelasticQ.append(dxsec_table['Q']) - - self.inelastic_legendre.append(dxsec_table['legendre']) - if(dxsec_table['legendre']): + xsec_table = ENDF_data[f"inelastic_xsec_n{i_inelastic + 1}"] + self.isigma.append( + interpolate_1d( + xsec_table["E"], xsec_table["sig"], method="linear", bounds_error=False, fill_value=0.0 + ) + ) + + dxsec_table = ENDF_data[f"inelastic_dxsec_n{i_inelastic + 1}"] + self.inelasticQ.append(dxsec_table["Q"]) + + self.inelastic_legendre.append(dxsec_table["legendre"]) + if dxsec_table["legendre"]: idx_spline = [unity] - for i in range(dxsec_table['N_l']): - idx_spline.append(interpolate_1d(dxsec_table['E'],dxsec_table['a_l'][:,i],method='linear',bounds_error=False,fill_value=0.0)) + for i in range(dxsec_table["N_l"]): + idx_spline.append( + interpolate_1d( + dxsec_table["E"], + dxsec_table["a_l"][:, i], + method="linear", + bounds_error=False, + fill_value=0.0, + ) + ) self.legendre_idx_spline.append(idx_spline) self.inelastic_SDX_table.append(None) else: self.legendre_idx_spline.append(None) - self.inelastic_SDX_table.append(dxsec_table['SDX']) + self.inelastic_SDX_table.append(dxsec_table["SDX"]) - self.Ein = None + self.Ein = None self.Eout = None self.vvec = None @@ -89,94 +122,96 @@ def __init__(self,label,json): # Stationary ion scattered spectral shapes # ############################################ - def init_energy_grids(self,Eout,Ein): + def init_energy_grids(self, Eout, Ein): self.Eout = Eout - self.Ein = Ein + self.Ein = Ein - def init_station_scatter_matrices(self,Nm=100): + def init_station_scatter_matrices(self, Nm=100): self.init_station_elastic_scatter() - if(self.l_n2n): + if self.l_n2n: self.init_n2n_ddxs(Nm) - if(self.l_inelastic): + if self.l_inelastic: self.init_station_inelastic_scatter() # Elastic scatter matrix def init_station_elastic_scatter(self): - Ei,Eo = np.meshgrid(self.Ein,self.Eout) - muc = col.muc(self.A,Ei,Eo,1.0,-1.0,0.0) - sigma = self.sigma(self.Ein) - self.elastic_mu0 = col.mu_out(self.A,Ei,Eo,0.0) - if(self.elastic_legendre): - Tlcoeff,Nl = xs.interp_Tlcoeff(self.legendre_dx_spline,self.Ein) - Tlcoeff_interp = 0.5*(2*np.arange(0,Nl)+1)*Tlcoeff - dsdO = xs.diffxsec_legendre_eval(sigma,muc,Tlcoeff_interp) + Ei, Eo = np.meshgrid(self.Ein, self.Eout) + muc = col.muc(self.A, Ei, Eo, 1.0, -1.0, 0.0) + sigma = self.sigma(self.Ein) + self.elastic_mu0 = col.mu_out(self.A, Ei, Eo, 0.0) + if self.elastic_legendre: + Tlcoeff, Nl = xs.interp_Tlcoeff(self.legendre_dx_spline, self.Ein) + Tlcoeff_interp = 0.5 * (2 * np.arange(0, Nl) + 1) * Tlcoeff + dsdO = xs.diffxsec_legendre_eval(sigma, muc, Tlcoeff_interp) else: - dsdO = xs.diffxsec_table_eval(sigma,muc,Ei,self.elastic_SDX_table) - jacob = col.g(self.A,Ei,Eo,1.0,-1.0,0.0) - self.elastic_dNdEdmu = jacob*dsdO + dsdO = xs.diffxsec_table_eval(sigma, muc, Ei, self.elastic_SDX_table) + jacob = col.g(self.A, Ei, Eo, 1.0, -1.0, 0.0) + self.elastic_dNdEdmu = jacob * dsdO # Inelastic scatter matrix # Currently uses classical kinematics def init_station_inelastic_scatter(self): - Ei,Eo = np.meshgrid(self.Ein,self.Eout) + Ei, Eo = np.meshgrid(self.Ein, self.Eout) self.inelastic_mu0 = [] self.inelastic_dNdEdmu = [] for i_inelastic in range(self.n_inelastic): - kin_a2 = (self.A/(self.A+1))**2*(1.0+(self.A+1)/self.A*self.inelasticQ[i_inelastic]/Ei) + kin_a2 = (self.A / (self.A + 1)) ** 2 * (1.0 + (self.A + 1) / self.A * self.inelasticQ[i_inelastic] / Ei) kin_a2_safe = kin_a2.copy() kin_a2_safe[kin_a2_safe < 0.0] = 1.0 - kin_a = np.sqrt(kin_a2_safe) - kin_b = 1.0/(self.A+1) - muc = ((Eo/Ei)-kin_a**2-kin_b**2)/(2*kin_a*kin_b) - sigma = self.isigma[i_inelastic](self.Ein) - inelastic_mu0 = (np.sqrt(Eo/Ei)-(kin_a**2-kin_b**2)*np.sqrt(Ei/Eo))/(2*kin_b) + kin_a = np.sqrt(kin_a2_safe) + kin_b = 1.0 / (self.A + 1) + muc = ((Eo / Ei) - kin_a**2 - kin_b**2) / (2 * kin_a * kin_b) + sigma = self.isigma[i_inelastic](self.Ein) + inelastic_mu0 = (np.sqrt(Eo / Ei) - (kin_a**2 - kin_b**2) * np.sqrt(Ei / Eo)) / (2 * kin_b) inelastic_mu0[kin_a2 < 0.0] = 0.0 self.inelastic_mu0.append(inelastic_mu0) - if(self.inelastic_legendre[i_inelastic]): - Tlcoeff,Nl = xs.interp_Tlcoeff(self.legendre_idx_spline[i_inelastic],self.Ein) - Tlcoeff_interp = 0.5*(2*np.arange(0,Nl)+1)*Tlcoeff - dsdO = xs.diffxsec_legendre_eval(sigma,muc,Tlcoeff_interp) + if self.inelastic_legendre[i_inelastic]: + Tlcoeff, Nl = xs.interp_Tlcoeff(self.legendre_idx_spline[i_inelastic], self.Ein) + Tlcoeff_interp = 0.5 * (2 * np.arange(0, Nl) + 1) * Tlcoeff + dsdO = xs.diffxsec_legendre_eval(sigma, muc, Tlcoeff_interp) else: - dsdO = xs.diffxsec_table_eval(sigma,muc,Ei,self.inelastic_SDX_table[i_inelastic]) + dsdO = xs.diffxsec_table_eval(sigma, muc, Ei, self.inelastic_SDX_table[i_inelastic]) + + jacob = 2.0 / ((kin_a + kin_b) ** 2 - (kin_a - kin_b) ** 2) / Ei + inelastic_dNdEdmu = jacob * dsdO - jacob = 2.0/((kin_a+kin_b)**2-(kin_a-kin_b)**2)/Ei - inelastic_dNdEdmu = jacob*dsdO - inelastic_dNdEdmu[kin_a2 < 0.0] = 0.0 self.inelastic_dNdEdmu.append(inelastic_dNdEdmu) - def init_n2n_ddxs(self,Nm=100): - self.n2n_mu = np.linspace(-1.0,1.0,Nm) - self.n2n_ddx.regular_grid(self.Ein,self.n2n_mu,self.Eout) + def init_n2n_ddxs(self, Nm=100): + self.n2n_mu = np.linspace(-1.0, 1.0, Nm) + self.n2n_ddx.regular_grid(self.Ein, self.n2n_mu, self.Eout) - def calc_dNdEs(self,I_E,rhoL_func): - self.calc_station_elastic_dNdE(I_E,rhoL_func) - if(self.l_n2n): - self.calc_n2n_dNdE(I_E,rhoL_func) - if(self.l_inelastic): - self.calc_station_inelastic_dNdE(I_E,rhoL_func) + def calc_dNdEs(self, I_E, rhoL_func): + self.calc_station_elastic_dNdE(I_E, rhoL_func) + if self.l_n2n: + self.calc_n2n_dNdE(I_E, rhoL_func) + if self.l_inelastic: + self.calc_station_inelastic_dNdE(I_E, rhoL_func) # Spectrum produced by scattering of incoming isotropic neutron source I_E with normalised areal density asymmetry rhoR_asym_func - def calc_station_elastic_dNdE(self,I_E,rhoL_func): + def calc_station_elastic_dNdE(self, I_E, rhoL_func): rhoL_asym = rhoL_func(self.elastic_mu0) - self.elastic_dNdE = np.trapz(self.elastic_dNdEdmu*rhoL_asym*I_E[None,:],self.Ein,axis=1) + self.elastic_dNdE = np.trapezoid(self.elastic_dNdEdmu * rhoL_asym * I_E[None, :], self.Ein, axis=1) - def calc_station_inelastic_dNdE(self,I_E,rhoL_func): + def calc_station_inelastic_dNdE(self, I_E, rhoL_func): self.inelastic_dNdE = np.zeros(self.Eout.shape[0]) for i_inelastic in range(self.n_inelastic): rhoL_asym = rhoL_func(self.inelastic_mu0[i_inelastic]) - self.inelastic_dNdE += np.trapz(self.inelastic_dNdEdmu[i_inelastic]*rhoL_asym*I_E[None,:],self.Ein,axis=1) + self.inelastic_dNdE += np.trapezoid( + self.inelastic_dNdEdmu[i_inelastic] * rhoL_asym * I_E[None, :], self.Ein, axis=1 + ) - def calc_n2n_dNdE(self,I_E,rhoL_func): + def calc_n2n_dNdE(self, I_E, rhoL_func): rhoL_asym = rhoL_func(self.n2n_mu) - grid_dNdE = np.trapz(self.n2n_ddx.rgrid*rhoL_asym[None,:,None],self.n2n_mu,axis=1) - self.n2n_dNdE = np.trapz(I_E[:,None]*grid_dNdE,self.Ein,axis=0) + grid_dNdE = np.trapezoid(self.n2n_ddx.rgrid * rhoL_asym[None, :, None], self.n2n_mu, axis=1) + self.n2n_dNdE = np.trapezoid(I_E[:, None] * grid_dNdE, self.Ein, axis=0) - def rhoR_2_A1s(self,rhoR): - mbar = self.A*Mn_kg - A_1S = rhoR*(sigmabarn/mbar) + def rhoR_2_A1s(self, rhoR): + mbar = self.A * Mn_kg + A_1S = rhoR * (sigmabarn / mbar) return A_1S # # Spectrum produced by scattering of incoming neutron source with anisotropic birth spectrum @@ -193,61 +228,68 @@ def rhoR_2_A1s(self,rhoR): # I_E_aniso = b_spec(Ei,prim_mean,var_iso) # dsdO = diffxsec_legendre_eval(sigma,muc,Tlcoeff_interp) # jacob = col.g(self.A,Ei,Eo,1.0,-1.0,0.0) - # res = np.trapz(jacob*dsdO*I_E_aniso*rhoR_asym,Ein,axis=-1) + # res = np.trapezoid(jacob*dsdO*I_E_aniso*rhoR_asym,Ein,axis=-1) # return res ##################################################### # Inclusion of ion velocities to scattering kernels # ##################################################### - def full_scattering_matrix_create(self,vvec): + def full_scattering_matrix_create(self, vvec): self.vvec = vvec - Eo,vv,Ei = np.meshgrid(self.Eout,vvec,self.Ein,indexing='ij') + Eo, vv, Ei = np.meshgrid(self.Eout, vvec, self.Ein, indexing="ij") # Reverse velocity direction so +ve vf is implosion # Choose this way round so vf is +ve if shell coming TOWARDS detector - vf = -vv - muout = col.mu_out(self.A,Ei,Eo,vf) - jacob = col.g(self.A,Ei,Eo,1.0,muout,vf) - flux_change = col.flux_change(Ei,1.0,vf) + vf = -vv + muout = col.mu_out(self.A, Ei, Eo, vf) + jacob = col.g(self.A, Ei, Eo, 1.0, muout, vf) + flux_change = col.flux_change(Ei, 1.0, vf) # Integrand of Eq. 8 in A. J. Crilly 2019 PoP - dsigdOmega = xs.dsigdOmega(self.A,Ei,Eo,self.Ein,1.0,muout,vf,self) + dsigdOmega = xs.dsigdOmega(self.A, Ei, Eo, self.Ein, 1.0, muout, vf, self) - self.full_scattering_M = flux_change*dsigdOmega*jacob + self.full_scattering_M = flux_change * dsigdOmega * jacob self.full_scattering_mu = muout - self.rhoL_mult = np.ones_like(muout) + self.rhoL_mult = np.ones_like(muout) - def scattering_matrix_apply_rhoLfunc(self,rhoL_func): + def scattering_matrix_apply_rhoLfunc(self, rhoL_func): # Find multiplicative factor for areal density asymmetries - self.rhoL_mult = rhoL_func(self.full_scattering_mu) + self.rhoL_mult = rhoL_func(self.full_scattering_mu) # Integrate out the birth neutron spectrum - def matrix_primspec_int(self,I_E): - self.M_prim = np.trapz(self.rhoL_mult*self.full_scattering_M*I_E[None,None,:],self.Ein,axis=2) + def matrix_primspec_int(self, I_E): + self.M_prim = np.trapezoid(self.rhoL_mult * self.full_scattering_M * I_E[None, None, :], self.Ein, axis=2) # Integrate out the ion velocity distribution - def matrix_interpolate_gaussian(self,E,vbar,dv): + def matrix_interpolate_gaussian(self, E, vbar, dv): # Integrating over Gaussian - gauss = np.exp(-(self.vvec-vbar)**2/2.0/(dv**2))/np.sqrt(2*np.pi)/dv - M_v = np.trapz(self.M_prim*gauss[None,:],self.vvec,axis=1) + gauss = np.exp(-((self.vvec - vbar) ** 2) / 2.0 / (dv**2)) / np.sqrt(2 * np.pi) / dv + M_v = np.trapezoid(self.M_prim * gauss[None, :], self.vvec, axis=1) # Interpolate to energy points E - interp = interpolate_1d(self.Eout,M_v,method='linear',bounds_error=False) + interp = interpolate_1d(self.Eout, M_v, method="linear", bounds_error=False) return interp(E) + # Load in TT spectrum # Based on Appelbe, stationary emitter, temperature range between 1 and 10 keV # https://www.sciencedirect.com/science/article/pii/S1574181816300295 # N.B. requires some unit conversion to uniform eV -TT_data = np.loadtxt(data_dir + "TT_spec_temprange.txt") -TT_spec_E = TT_data[:,0]*1e6 # MeV to eV -TT_spec_T = np.linspace(1.0,20.0,40)*1e3 # keV to eV -TT_spec_dNdE = TT_data[:,1:]/1e6 # 1/MeV to 1/eV -TT_2dinterp = interpolate_2d(TT_spec_E,TT_spec_T,TT_spec_dNdE,method='linear',bounds_error=False,fill_value=0.0) +TT_data = np.loadtxt(data_dir + "TT_spec_temprange.txt") +TT_spec_E = TT_data[:, 0] * 1e6 # MeV to eV +TT_spec_T = np.linspace(1.0, 20.0, 40) * 1e3 # keV to eV +TT_spec_dNdE = TT_data[:, 1:] / 1e6 # 1/MeV to 1/eV +TT_2dinterp = interpolate_2d(TT_spec_E, TT_spec_T, TT_spec_dNdE, method="linear", bounds_error=False, fill_value=0.0) # TT reactivity -TT_reac_McNally_data = np.loadtxt(data_dir + "TT_reac_McNally.dat") # sigmav im m^3/s # From https://www.osti.gov/servlets/purl/5992170 - N.B. not in agreement with experimental measurements -TT_reac_McNally_spline = interpolate_1d(TT_reac_McNally_data[:,0],TT_reac_McNally_data[:,1],method='linear',bounds_error=False,fill_value=0.0) -TT_reac_Hale_data = np.loadtxt(data_dir + "TT_reac_Hale.dat") # T in MeV, sigmav im cm^3/s # From Hale -TT_reac_Hale_spline = interpolate_1d(TT_reac_Hale_data[:,0]*1e3,TT_reac_Hale_data[:,1]*1e-6,method='linear',bounds_error=False,fill_value=0.0) +TT_reac_McNally_data = np.loadtxt( + data_dir + "TT_reac_McNally.dat" +) # sigmav im m^3/s # From https://www.osti.gov/servlets/purl/5992170 - N.B. not in agreement with experimental measurements +TT_reac_McNally_spline = interpolate_1d( + TT_reac_McNally_data[:, 0], TT_reac_McNally_data[:, 1], method="linear", bounds_error=False, fill_value=0.0 +) +TT_reac_Hale_data = np.loadtxt(data_dir + "TT_reac_Hale.dat") # T in MeV, sigmav im cm^3/s # From Hale +TT_reac_Hale_spline = interpolate_1d( + TT_reac_Hale_data[:, 0] * 1e3, TT_reac_Hale_data[:, 1] * 1e-6, method="linear", bounds_error=False, fill_value=0.0 +) # TT_reac_data = np.loadtxt(data_dir + "TT_reac_ENDF.dat") # sigmav im m^3/s # From ENDF # TT_reac_spline = interpolate_1d(TT_reac_data[:,0],TT_reac_data[:,1],method='linear',bounds_error=False,fill_value=0.0) @@ -260,54 +302,63 @@ def matrix_interpolate_gaussian(self,E,vbar,dv): # Caughlan & Fowler: https://doi.org/10.1146/annurev.aa.13.090175.000441 # McNally: https://www.osti.gov/servlets/purl/5992170 + # Output in m3/s, Ti in eV -def reac_DT(Ti,model='BoschHale'): - Ti_kev = Ti/1e3 - if(model == 'BoschHale'): +def reac_DT(Ti, model="BoschHale"): + Ti_kev = Ti / 1e3 + if model == "BoschHale": # Bosch Hale DT and DD reactivities # Taken from Atzeni & Meyer ter Vehn page 19 C1 = 643.41e-22 - xi = 6.6610*Ti_kev**(-0.333333333) - eta = 1-np.polyval([-0.10675e-3,4.6064e-3,15.136e-3,0.0e0],Ti_kev)/np.polyval([0.01366e-3,13.5e-3,75.189e-3,1.0e0],Ti_kev) - return C1*eta**(-0.833333333)*xi**2*np.exp(-3*eta**(0.333333333)*xi) - elif(model == 'CaughlanFowler'): - T9 = (Ti_kev*sc.e*1e3/sc.k)/1e9 - T9_1third = T9**(1./3.) - poly = np.polyval([17.24,10.52,1.16,1.80,0.092,1.0],T9_1third) - return (1/sc.N_A)*(8.09e4*poly*np.exp(-4.524/T9**(1./3.)-(T9/0.120)**2) + 8.73e2*np.exp(-0.523/T9))/T9**(2./3.) + xi = 6.6610 * Ti_kev ** (-0.333333333) + eta = 1 - np.polyval([-0.10675e-3, 4.6064e-3, 15.136e-3, 0.0e0], Ti_kev) / np.polyval( + [0.01366e-3, 13.5e-3, 75.189e-3, 1.0e0], Ti_kev + ) + return C1 * eta ** (-0.833333333) * xi**2 * np.exp(-3 * eta ** (0.333333333) * xi) + elif model == "CaughlanFowler": + T9 = (Ti_kev * sc.e * 1e3 / sc.k) / 1e9 + T9_1third = T9 ** (1.0 / 3.0) + poly = np.polyval([17.24, 10.52, 1.16, 1.80, 0.092, 1.0], T9_1third) + return ( + (1 / sc.N_A) + * (8.09e4 * poly * np.exp(-4.524 / T9 ** (1.0 / 3.0) - (T9 / 0.120) ** 2) + 8.73e2 * np.exp(-0.523 / T9)) + / T9 ** (2.0 / 3.0) + ) else: - print(f'WARNING: DT model name ({model}) not recognised! Default to 0') + print(f"WARNING: DT model name ({model}) not recognised! Default to 0") return np.zeros_like(Ti) - -def reac_DD(Ti,model='BoschHale'): - Ti_kev = Ti/1e3 - if(model == 'BoschHale'): + + +def reac_DD(Ti, model="BoschHale"): + Ti_kev = Ti / 1e3 + if model == "BoschHale": # Bosch Hale DT and DD reactivities # Taken from Atzeni & Meyer ter Vehn page 19 C1 = 3.5741e-22 - xi = 6.2696*Ti_kev**(-0.333333333) - eta = 1-np.polyval([5.8577e-3,0.0e0],Ti_kev)/np.polyval([-0.002964e-3,7.6822e-3,1.0e0],Ti_kev) - return C1*eta**(-0.833333333)*xi**2*np.exp(-3*eta**(0.333333333)*xi) - elif(model == 'CaughlanFowler'): - T9 = (Ti_kev*sc.e*1e3/sc.k)/1e9 - T9_1third = T9**(1./3.) - poly = np.polyval([-0.071,-0.041,0.6,0.876,0.098,1.0],T9_1third) - return (1/sc.N_A)*3.97e2/T9**(2./3.)*np.exp(-4.258/T9**(1./3.))*poly + xi = 6.2696 * Ti_kev ** (-0.333333333) + eta = 1 - np.polyval([5.8577e-3, 0.0e0], Ti_kev) / np.polyval([-0.002964e-3, 7.6822e-3, 1.0e0], Ti_kev) + return C1 * eta ** (-0.833333333) * xi**2 * np.exp(-3 * eta ** (0.333333333) * xi) + elif model == "CaughlanFowler": + T9 = (Ti_kev * sc.e * 1e3 / sc.k) / 1e9 + T9_1third = T9 ** (1.0 / 3.0) + poly = np.polyval([-0.071, -0.041, 0.6, 0.876, 0.098, 1.0], T9_1third) + return (1 / sc.N_A) * 3.97e2 / T9 ** (2.0 / 3.0) * np.exp(-4.258 / T9 ** (1.0 / 3.0)) * poly else: - print(f'WARNING: DD model name ({model}) not recognised! Default to 0') + print(f"WARNING: DD model name ({model}) not recognised! Default to 0") return np.zeros_like(Ti) - -def reac_TT(Ti,model='Hale'): - Ti_kev = Ti/1e3 - if(model == 'Hale'): + + +def reac_TT(Ti, model="Hale"): + Ti_kev = Ti / 1e3 + if model == "Hale": return TT_reac_Hale_spline(Ti_kev) - elif(model == 'McNally'): + elif model == "McNally": return TT_reac_McNally_spline(Ti_kev) - elif(model == 'CaughlanFowler'): - T9 = (Ti_kev*sc.e*1e3/sc.k)/1e9 - T9_1third = T9**(1./3.) - poly = np.polyval([0.225,0.148,-0.272,-0.455,0.086,1.0],T9_1third) - return (1/sc.N_A)*1.67e3/T9**(2./3.)*np.exp(-4.872/T9**(1./3.))*poly + elif model == "CaughlanFowler": + T9 = (Ti_kev * sc.e * 1e3 / sc.k) / 1e9 + T9_1third = T9 ** (1.0 / 3.0) + poly = np.polyval([0.225, 0.148, -0.272, -0.455, 0.086, 1.0], T9_1third) + return (1 / sc.N_A) * 1.67e3 / T9 ** (2.0 / 3.0) * np.exp(-4.872 / T9 ** (1.0 / 3.0)) * poly else: - print(f'WARNING: TT model name ({model}) not recognised! Default to 0') - return np.zeros_like(Ti) \ No newline at end of file + print(f"WARNING: TT model name ({model}) not recognised! Default to 0") + return np.zeros_like(Ti) diff --git a/src/NeSST/time_of_flight.py b/src/NeSST/time_of_flight.py index ee38b74..ace0f93 100644 --- a/src/NeSST/time_of_flight.py +++ b/src/NeSST/time_of_flight.py @@ -1,62 +1,70 @@ -from NeSST.constants import * -from NeSST.collisions import * -from NeSST.utils import * -from NeSST.endf_interface import retrieve_total_cross_section_from_ENDF_file +from dataclasses import dataclass, field +from typing import List + import numpy as np from scipy.integrate import cumulative_trapezoid as cumtrapz from scipy.integrate import quad from scipy.special import erf -from dataclasses import dataclass,field -from typing import List +from NeSST.collisions import * +from NeSST.constants import * +from NeSST.endf_interface import retrieve_total_cross_section_from_ENDF_file +from NeSST.utils import * + @dataclass class LOS_material_component: - number_fraction : float - A : float - ENDF_file : str + number_fraction: float + A: float + ENDF_file: str + @dataclass class LOS_material: - density : float - ntot : float = field(init=False) - mavg : float = field(init=False) - length : float - components : List[LOS_material_component] + density: float + ntot: float = field(init=False) + mavg: float = field(init=False) + length: float + components: List[LOS_material_component] def __post_init__(self): self.mavg = 0.0 for comp in self.components: - self.mavg += comp.number_fraction*comp.A + self.mavg += comp.number_fraction * comp.A self.mavg *= sc.atomic_mass - self.ntot = self.density/self.mavg + self.ntot = self.density / self.mavg -def get_LOS_attenuation(LOS_materials : List[LOS_material]): + +def get_LOS_attenuation(LOS_materials: List[LOS_material]): tau_interp_list = [] for LOS_material in LOS_materials: - ntot_barn = 1e-28*LOS_material.ntot - L = LOS_material.length + ntot_barn = 1e-28 * LOS_material.ntot + L = LOS_material.length for LOS_component in LOS_material.components: - ncomp = ntot_barn*LOS_component.number_fraction - E,sigma_tot = retrieve_total_cross_section_from_ENDF_file(LOS_component.ENDF_file) - tau_interp = interpolate_1d(E,L*ncomp*sigma_tot,method='linear') + ncomp = ntot_barn * LOS_component.number_fraction + E, sigma_tot = retrieve_total_cross_section_from_ENDF_file(LOS_component.ENDF_file) + tau_interp = interpolate_1d(E, L * ncomp * sigma_tot, method="linear") tau_interp_list.append(tau_interp) def LOS_attenuation(E): total_tau = tau_interp_list[0](E) - for i in range(1,len(tau_interp_list)): + for i in range(1, len(tau_interp_list)): total_tau += tau_interp_list[i](E) transmission = np.exp(-total_tau) return transmission - + return LOS_attenuation -def get_power_law_NLO(p,Enorm=E0_DT): - A = mat_dict['H'].sigma(Enorm)*Enorm**p + +def get_power_law_NLO(p, Enorm=E0_DT): + A = mat_dict["H"].sigma(Enorm) * Enorm**p + def power_law_NLO(E): - return mat_dict['H'].sigma(E)*E**p/A + return mat_dict["H"].sigma(E) * E**p / A + return power_law_NLO + def get_Verbinski_NLO(Enorm=E0_DT): """ Using equation (3) from @@ -70,22 +78,23 @@ def get_Verbinski_NLO(Enorm=E0_DT): ISSN 0168-9002, https://doi.org/10.1016/j.nima.2024.169779. """ - V_E,V_L = np.loadtxt(data_dir+'VerbinskiLproton.csv',delimiter=',',unpack=True) - cumulative_L = cumtrapz(y=np.insert(V_L,0,0.0),x=np.insert(V_E,0,0.0)) - L_integral = interpolate_1d(np.insert(V_E,0,0.0)*1e6,np.insert(cumulative_L,0,0.0),method='cubic') + V_E, V_L = np.loadtxt(data_dir + "VerbinskiLproton.csv", delimiter=",", unpack=True) + cumulative_L = cumtrapz(y=np.insert(V_L, 0, 0.0), x=np.insert(V_E, 0, 0.0)) + L_integral = interpolate_1d(np.insert(V_E, 0, 0.0) * 1e6, np.insert(cumulative_L, 0, 0.0), method="cubic") def Verbinski_NLO(E): - return mat_dict['H'].sigma(E)*L_integral(E)/E - + return mat_dict["H"].sigma(E) * L_integral(E) / E + A = Verbinski_NLO(Enorm) def norm_Verbinski_NLO(E): - return Verbinski_NLO(E)/A - + return Verbinski_NLO(E) / A + return norm_Verbinski_NLO -def get_BirksBetheBloch_NLO(akB,Enorm=E0_DT): - """ + +def get_BirksBetheBloch_NLO(akB, Enorm=E0_DT): + r""" akB in eV Combining: @@ -110,23 +119,27 @@ def get_BirksBetheBloch_NLO(akB,Enorm=E0_DT): """ def BirksBetheBloch_NLO(E): - L_integral = 0.5*E**2-akB*E-akB*(akB+E)*np.log(1.0+E/akB) - return mat_dict['H'].sigma(E)*L_integral/E - + L_integral = 0.5 * E**2 - akB * E - akB * (akB + E) * np.log(1.0 + E / akB) + return mat_dict["H"].sigma(E) * L_integral / E + A = BirksBetheBloch_NLO(Enorm) def norm_BirksBetheBloch_NLO(E): - return BirksBetheBloch_NLO(E)/A - + return BirksBetheBloch_NLO(E) / A + return norm_BirksBetheBloch_NLO -def get_BirksBethe_NLO(akB,I,mp = sc.m_p,Enorm=E0_DT): - """ + +def get_BirksBethe_NLO(akB, excitation_energy, mp=sc.m_p, Enorm=E0_DT): + r""" akB in eV Combining: kB in m/eV a in eV^2/m + + I = excitation_energy in eV + Bethe formula for stopping power: dE/dx = a/E * ln(4 me E / mp I) @@ -144,206 +157,246 @@ def get_BirksBethe_NLO(akB,I,mp = sc.m_p,Enorm=E0_DT): ISSN 0168-9002, https://doi.org/10.1016/j.nima.2024.169779. """ - Istar = I * mp / sc.m_e / 4.0 + Istar = excitation_energy * mp / sc.m_e / 4.0 def kB_dEdx(Ep): - Ep_lim = max([Ep,np.e*Istar]) - return akB/Ep_lim*np.log(Ep_lim/Istar) + Ep_lim = max([Ep, np.e * Istar]) + return akB / Ep_lim * np.log(Ep_lim / Istar) def dLdE(Ep): - return 1.0/(1.0+kB_dEdx(Ep)) - + return 1.0 / (1.0 + kB_dEdx(Ep)) + def L(Ep): - return quad(dLdE,0.0,Ep)[0] - + return quad(dLdE, 0.0, Ep)[0] + def L_integral_scalar(En): - return quad(L,0.0,En)[0] - + return quad(L, 0.0, En)[0] + L_integral = np.vectorize(L_integral_scalar) def BirksBethe_NLO(E): - return mat_dict['H'].sigma(E)*L_integral(E)/E - + return mat_dict["H"].sigma(E) * L_integral(E) / E + A = BirksBethe_NLO(Enorm) def norm_BirksBethe_NLO(E): - return BirksBethe_NLO(E)/A - - return norm_BirksBethe_NLO,kB_dEdx,L + return BirksBethe_NLO(E) / A + + return norm_BirksBethe_NLO, kB_dEdx, L + def get_unity_sensitivity(): def unity_sensitivity(En): return np.ones_like(En) + return unity_sensitivity + def combine_detector_sensitivities(model_list): def total_sensitivity(E): sensitivity = model_list[0](E) - for i in range(1,len(model_list)): + for i in range(1, len(model_list)): sensitivity *= model_list[i](E) return sensitivity + return total_sensitivity + def get_transit_time_tophat_IRF(scintillator_thickness): - def transit_time_tophat_IRF(t_detected,En): - vn = Ekin_2_beta(En,Mn)*c - t_transit = scintillator_thickness/vn - tt_d,tt_a = np.meshgrid(t_detected,t_detected,indexing='ij') - _,tt_t = np.meshgrid(t_detected,t_transit,indexing='ij') - R = np.eye(t_detected.shape[0])+np.heaviside(tt_d-tt_a,0.0)-np.heaviside(tt_d-(tt_a+tt_t),1.0) - R_norm = np.sum(R,axis=1) + def transit_time_tophat_IRF(t_detected, En): + vn = Ekin_2_beta(En, Mn) * c + t_transit = scintillator_thickness / vn + tt_d, tt_a = np.meshgrid(t_detected, t_detected, indexing="ij") + _, tt_t = np.meshgrid(t_detected, t_transit, indexing="ij") + R = np.eye(t_detected.shape[0]) + np.heaviside(tt_d - tt_a, 0.0) - np.heaviside(tt_d - (tt_a + tt_t), 1.0) + R_norm = np.sum(R, axis=1) # Catch the zeros R_norm[R_norm == 0] = 1 - R /= R_norm[:,None] + R /= R_norm[:, None] return R + return transit_time_tophat_IRF -def get_transit_time_tophat_w_decayingGaussian_IRF(scintillator_thickness,gaussian_FWHM,decay_time,sig_shift=2.0): + +def get_transit_time_tophat_w_decayingGaussian_IRF(scintillator_thickness, gaussian_FWHM, decay_time, sig_shift=2.0): """ - + See: https://pubs.aip.org/aip/rsi/article/68/1/610/1070804/Interpretation-of-neutron-time-of-flight-signals """ - sig = gaussian_FWHM/2.355 - shift_t = sig_shift*sig + sig = gaussian_FWHM / 2.355 + shift_t = sig_shift * sig + def filter(t): - t_shift = t-shift_t - erf_arg = (t_shift-sig**2/decay_time)/np.sqrt(2*sig**2) - gauss = np.exp(-t_shift/decay_time)*np.exp(0.5*sig**2/decay_time**2)/(2*decay_time)*(1+erf(erf_arg)) + t_shift = t - shift_t + erf_arg = (t_shift - sig**2 / decay_time) / np.sqrt(2 * sig**2) + gauss = ( + np.exp(-t_shift / decay_time) * np.exp(0.5 * sig**2 / decay_time**2) / (2 * decay_time) * (1 + erf(erf_arg)) + ) gauss[t < 0] = 0.0 return gauss + _tophat_IRF = get_transit_time_tophat_IRF(scintillator_thickness) - def transit_time_tophat_w_decayingGaussian_IRF(t_detected,En): - R = _tophat_IRF(t_detected,En) - t_filter = t_detected-0.5*(t_detected[-1]+t_detected[0]) + + def transit_time_tophat_w_decayingGaussian_IRF(t_detected, En): + R = _tophat_IRF(t_detected, En) + t_filter = t_detected - 0.5 * (t_detected[-1] + t_detected[0]) filt = filter(t_filter) - R = np.apply_along_axis(lambda m : np.convolve(m,filt,mode='same'), axis=0, arr=R) - R_norm = np.sum(R,axis=1) + R = np.apply_along_axis(lambda m: np.convolve(m, filt, mode="same"), axis=0, arr=R) + R_norm = np.sum(R, axis=1) # Catch the zeros R_norm[R_norm == 0] = 1 - R /= R_norm[:,None] + R /= R_norm[:, None] return R + return transit_time_tophat_w_decayingGaussian_IRF -def get_transit_time_tophat_w_gateddecayingGaussian_IRF(scintillator_thickness,sig,decay_time,shift_t,sig_turnon): + +def get_transit_time_tophat_w_gateddecayingGaussian_IRF(scintillator_thickness, sig, decay_time, shift_t, sig_turnon): """ - + See: https://pubs.aip.org/aip/rsi/article/68/1/610/1070804/Interpretation-of-neutron-time-of-flight-signals """ def gate(x): - g = 2.0/(1.0+np.exp(-x))-1.0 + g = 2.0 / (1.0 + np.exp(-x)) - 1.0 g[x < 0.0] = 0.0 return g - + def filter(t): - t_shift = t-shift_t - erf_arg = (t_shift-sig**2/decay_time)/np.sqrt(2*sig**2) - gauss = np.exp(-t_shift/decay_time)*np.exp(0.5*sig**2/decay_time**2)/(2*decay_time)*(1+erf(erf_arg)) - gauss *= gate(t/sig_turnon) - return gauss/np.trapz(gauss,x=t) - + t_shift = t - shift_t + erf_arg = (t_shift - sig**2 / decay_time) / np.sqrt(2 * sig**2) + gauss = ( + np.exp(-t_shift / decay_time) * np.exp(0.5 * sig**2 / decay_time**2) / (2 * decay_time) * (1 + erf(erf_arg)) + ) + gauss *= gate(t / sig_turnon) + return gauss / np.trapezoid(gauss, x=t) + _tophat_IRF = get_transit_time_tophat_IRF(scintillator_thickness) - def transit_time_tophat_w_doubledecayingGaussian_IRF(t_detected,En): - R = _tophat_IRF(t_detected,En) - t_filter = t_detected-0.5*(t_detected[-1]+t_detected[0]) + + def transit_time_tophat_w_doubledecayingGaussian_IRF(t_detected, En): + R = _tophat_IRF(t_detected, En) + t_filter = t_detected - 0.5 * (t_detected[-1] + t_detected[0]) filt = filter(t_filter) - R = np.apply_along_axis(lambda m : np.convolve(m,filt,mode='same'), axis=0, arr=R) - R_norm = np.sum(R,axis=1) + R = np.apply_along_axis(lambda m: np.convolve(m, filt, mode="same"), axis=0, arr=R) + R_norm = np.sum(R, axis=1) # Catch the zeros R_norm[R_norm == 0] = 1 - R /= R_norm[:,None] + R /= R_norm[:, None] return R + return transit_time_tophat_w_doubledecayingGaussian_IRF -def get_transit_time_tophat_w_doubledecayingGaussian_IRF(scintillator_thickness,gaussian_FWHM,decay_times,comp_1_frac,sig_shift=2.0): + +def get_transit_time_tophat_w_doubledecayingGaussian_IRF( + scintillator_thickness, gaussian_FWHM, decay_times, comp_1_frac, sig_shift=2.0 +): """ - + See: https://pubs.aip.org/aip/rsi/article/68/1/610/1070804/Interpretation-of-neutron-time-of-flight-signals """ - sig = gaussian_FWHM/2.355 - shift_t = sig_shift*sig + sig = gaussian_FWHM / 2.355 + shift_t = sig_shift * sig + + def single_decay_comp(t, decay_time): + t_shift = t - shift_t + erf_arg = (t_shift - sig**2 / decay_time) / np.sqrt(2 * sig**2) + gauss = ( + np.exp(-t_shift / decay_time) * np.exp(0.5 * sig**2 / decay_time**2) / (2 * decay_time) * (1 + erf(erf_arg)) + ) + gauss *= 0.5 * (1 + erf(t / sig)) + return gauss - def single_decay_comp(t,decay_time): - t_shift = t-shift_t - erf_arg = (t_shift-sig**2/decay_time)/np.sqrt(2*sig**2) - gauss = np.exp(-t_shift/decay_time)*np.exp(0.5*sig**2/decay_time**2)/(2*decay_time)*(1+erf(erf_arg)) - gauss *= 0.5*(1+erf(t/sig)) - return gauss - def filter(t): - total = comp_1_frac*single_decay_comp(t,decay_times[0])+(1-comp_1_frac)*single_decay_comp(t,decay_times[1]) - return total/np.trapz(total,x=t) - + total = comp_1_frac * single_decay_comp(t, decay_times[0]) + (1 - comp_1_frac) * single_decay_comp( + t, decay_times[1] + ) + return total / np.trapezoid(total, x=t) + _tophat_IRF = get_transit_time_tophat_IRF(scintillator_thickness) - def transit_time_tophat_w_doubledecayingGaussian_IRF(t_detected,En): - R = _tophat_IRF(t_detected,En) - t_filter = t_detected-0.5*(t_detected[-1]+t_detected[0]) + + def transit_time_tophat_w_doubledecayingGaussian_IRF(t_detected, En): + R = _tophat_IRF(t_detected, En) + t_filter = t_detected - 0.5 * (t_detected[-1] + t_detected[0]) filt = filter(t_filter) - R = np.apply_along_axis(lambda m : np.convolve(m,filt,mode='same'), axis=0, arr=R) - R_norm = np.sum(R,axis=1) + R = np.apply_along_axis(lambda m: np.convolve(m, filt, mode="same"), axis=0, arr=R) + R_norm = np.sum(R, axis=1) # Catch the zeros R_norm[R_norm == 0] = 1 - R /= R_norm[:,None] + R /= R_norm[:, None] return R + return transit_time_tophat_w_doubledecayingGaussian_IRF -def get_transit_time_tophat_w_tGaussian_IRF(scintillator_thickness,gaussian_FWHM,tgaussian_peak_pos): - sig = gaussian_FWHM/2.355 - mu = (tgaussian_peak_pos**2-sig**2)/tgaussian_peak_pos + +def get_transit_time_tophat_w_tGaussian_IRF(scintillator_thickness, gaussian_FWHM, tgaussian_peak_pos): + sig = gaussian_FWHM / 2.355 + mu = (tgaussian_peak_pos**2 - sig**2) / tgaussian_peak_pos + def filter(t): - gauss = t*np.exp(-0.5*((t-mu)/sig)**2) + gauss = t * np.exp(-0.5 * ((t - mu) / sig) ** 2) gauss[t < 0] = 0.0 return gauss + _tophat_IRF = get_transit_time_tophat_IRF(scintillator_thickness) - def transit_time_tophat_w_tGaussian_IRF(t_detected,En): - R = _tophat_IRF(t_detected,En) - t_filter = t_detected-0.5*(t_detected[-1]+t_detected[0]) + + def transit_time_tophat_w_tGaussian_IRF(t_detected, En): + R = _tophat_IRF(t_detected, En) + t_filter = t_detected - 0.5 * (t_detected[-1] + t_detected[0]) filt = filter(t_filter) - R = np.apply_along_axis(lambda m : np.convolve(m,filt,mode='same'), axis=0, arr=R) - R_norm = np.sum(R,axis=1) + R = np.apply_along_axis(lambda m: np.convolve(m, filt, mode="same"), axis=0, arr=R) + R_norm = np.sum(R, axis=1) # Catch the zeros R_norm[R_norm == 0] = 1 - R /= R_norm[:,None] + R /= R_norm[:, None] return R + return transit_time_tophat_w_tGaussian_IRF -class nToF: - def __init__(self,distance,sensitivity,instrument_response_function - ,normtime_start=5.0,normtime_end=20.0,normtime_N=2048,detector_normtime=None): +class nToF: + def __init__( + self, + distance, + sensitivity, + instrument_response_function, + normtime_start=5.0, + normtime_end=20.0, + normtime_N=2048, + detector_normtime=None, + ): self.distance = distance self.sensitivity = sensitivity self.instrument_response_function = instrument_response_function - if(detector_normtime is None): - self.detector_normtime = np.linspace(normtime_start,normtime_end,normtime_N) + if detector_normtime is None: + self.detector_normtime = np.linspace(normtime_start, normtime_end, normtime_N) else: self.detector_normtime = detector_normtime - self.detector_time = self.detector_normtime*self.distance/c + self.detector_time = self.detector_normtime * self.distance / c # Init instrument response values self.compute_instrument_response() def compute_instrument_response(self): - self.En_det = beta_2_Ekin(1.0/self.detector_normtime,Mn) - - self.dEdt = Jacobian_dEdnorm_t(self.En_det,Mn) + self.En_det = beta_2_Ekin(1.0 / self.detector_normtime, Mn) + + self.dEdt = Jacobian_dEdnorm_t(self.En_det, Mn) self.sens = self.sensitivity(self.En_det) - self.R = self.instrument_response_function(self.detector_time,self.En_det) + self.R = self.instrument_response_function(self.detector_time, self.En_det) + + def get_dNdt(self, En, dNdE): + dNdE_interp = np.interp(self.En_det, En, dNdE, left=0.0, right=0.0) + return dNdE_interp * self.dEdt - def get_dNdt(self,En,dNdE): - dNdE_interp = np.interp(self.En_det,En,dNdE,left=0.0,right=0.0) - return dNdE_interp*self.dEdt + def get_signal(self, En, dNdE): + dNdt = self.get_dNdt(En, dNdE) - def get_signal(self,En,dNdE): - dNdt = self.get_dNdt(En,dNdE) + return self.detector_time, self.detector_normtime, np.matmul(self.R, self.sens * dNdt) - return self.detector_time,self.detector_normtime,np.matmul(self.R,self.sens*dNdt) - - def get_signal_no_IRF(self,En,dNdE): - dNdt = self.get_dNdt(En,dNdE) + def get_signal_no_IRF(self, En, dNdE): + dNdt = self.get_dNdt(En, dNdE) - return self.detector_time,self.detector_normtime,dNdt + return self.detector_time, self.detector_normtime, dNdt diff --git a/src/NeSST/utils.py b/src/NeSST/utils.py index dbcbe18..4075cf3 100644 --- a/src/NeSST/utils.py +++ b/src/NeSST/utils.py @@ -1,5 +1,5 @@ -from scipy.interpolate import RegularGridInterpolator,interp1d import numpy as np +from scipy.interpolate import RegularGridInterpolator, interp1d """ @@ -7,60 +7,68 @@ """ + def safe_array(x): """Converts floats or lists to appropriate sized arrays""" - if not hasattr(x,'ndim'): + if not hasattr(x, "ndim"): x = np.array(x) return np.atleast_1d(x) -def interpolate_1d(x_points,values,method='linear',bounds_error=True,fill_value=np.nan,axis=None): + +def interpolate_1d(x_points, values, method="linear", bounds_error=True, fill_value=np.nan, axis=None): x_points = safe_array(x_points) assert x_points.ndim == 1 - if(axis is None): + if axis is None: axis = -1 - return interp1d(x_points,values,kind=method,bounds_error=bounds_error,fill_value=fill_value,axis=axis) + return interp1d(x_points, values, kind=method, bounds_error=bounds_error, fill_value=fill_value, axis=axis) """ Below is some future proofing where we attempt to reclaim some of interp1d's behaviour with RegularGridInterpolator It is less performant so is not used, but if interp1d becomes deprecated it might be needed... """ - if(values.ndim == 1): - if(np.all(np.isclose(x_points, x_points[0]))): + if values.ndim == 1: + if np.all(np.isclose(x_points, x_points[0])): # Sometimes the array is just a series of the same points... def point_interpolant(x): - y = fill_value*np.ones_like(x) + y = fill_value * np.ones_like(x) y[np.isclose(x, x_points[0])] = values[0] return y + return point_interpolant else: # Sometimes ENDF interpreted format has repeated elements... - if(not np.all(x_points[1:] > x_points[:-1])): - x_points,unique_indices = np.unique(x_points,return_index=True) + if not np.all(x_points[1:] > x_points[:-1]): + x_points, unique_indices = np.unique(x_points, return_index=True) values = values[unique_indices] # We often want to call 1D interpolators with non 1D shaped inputs - RGI = RegularGridInterpolator((x_points,),values,method,bounds_error,fill_value) + RGI = RegularGridInterpolator((x_points,), values, method, bounds_error, fill_value) + def dimension_matching_interp(x): x = np.atleast_1d(x) xshape = x.shape y = RGI(x.flatten()) return y.reshape(xshape) + return dimension_matching_interp else: # Currently doesn't seem to be axis argument support for RegularGridInterpolator... - return interp1d(x_points,values,kind=method,bounds_error=bounds_error,fill_value=fill_value,axis=axis) + return interp1d(x_points, values, kind=method, bounds_error=bounds_error, fill_value=fill_value, axis=axis) + -def interpolate_2d(x_points,y_points,values,method='linear',bounds_error=True,fill_value=np.nan): +def interpolate_2d(x_points, y_points, values, method="linear", bounds_error=True, fill_value=np.nan): x_points = safe_array(x_points) y_points = safe_array(y_points) assert x_points.ndim == 1 assert y_points.ndim == 1 - RGI = RegularGridInterpolator((x_points,y_points),values,method,bounds_error,fill_value) - def dimension_matching_interp(x,y): + RGI = RegularGridInterpolator((x_points, y_points), values, method, bounds_error, fill_value) + + def dimension_matching_interp(x, y): x = safe_array(x) y = safe_array(y) assert x.ndim == 1 assert y.ndim == 1 - xx,yy = np.meshgrid(x,y,indexing='ij') - f = RGI((xx,yy)) + xx, yy = np.meshgrid(x, y, indexing="ij") + f = RGI((xx, yy)) return np.squeeze(f) + return dimension_matching_interp diff --git a/tests/test_core.py b/tests/test_core.py index 0a0c564..6f6a63d 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -1,6 +1,6 @@ -import pytest -import numpy as np import NeSST as nst +import numpy as np +import pytest def test_DTprimspecmoments_mean(): @@ -36,14 +36,16 @@ def test_DTprimspecmoments_mean_with_tion(): assert DTmean_cold < DTmean_hot # units eV + def test_DTprimspecmoments_variance_with_tion(): # checks the relative magnitude of the var _, _, DTvar_cold = nst.DTprimspecmoments(Tion=5.0e3) # units eV - _, _,DTvar_hot = nst.DTprimspecmoments(Tion=10.0e3) # units eV + _, _, DTvar_hot = nst.DTprimspecmoments(Tion=10.0e3) # units eV assert DTvar_cold < DTvar_hot # units eV**2 + def test_DDprimspecmoments_variance_with_tion(): # checks the relative magnitude of the var @@ -52,22 +54,24 @@ def test_DDprimspecmoments_variance_with_tion(): assert DDvar_cold < DDvar_hot # units eV**2 + def test_DDprimspecmoments_variance_relative_size(): # checks the relative magnitude of the var DDmean, DDstddev, DDvar = nst.DDprimspecmoments(Tion=5.0e3) # units eV # Check that the standard deviation is about 3% of the mean value - assert np.isclose((100/DDmean)*DDstddev, 3, atol=0.3) + assert np.isclose((100 / DDmean) * DDstddev, 3, atol=0.3) # Check variance is standard deviation squared assert np.isclose(DDvar, DDstddev**2) + def test_DTprimspecmoments_variance_relative_size(): # checks the relative magnitude of the var DTmean, DTstddev, DTvar = nst.DTprimspecmoments(Tion=5.0e3) # units eV # Check that the standard deviation is about 3% of the mean value - assert np.isclose((100/DTmean)*DTstddev, 1, atol=0.3) + assert np.isclose((100 / DTmean) * DTstddev, 1, atol=0.3) # Check variance is standard deviation squared - assert np.isclose(DTvar, DTstddev**2) \ No newline at end of file + assert np.isclose(DTvar, DTstddev**2) diff --git a/tests/test_endf.py b/tests/test_endf.py index e91c8bd..01f5d24 100644 --- a/tests/test_endf.py +++ b/tests/test_endf.py @@ -1,14 +1,14 @@ -import pytest -import numpy as np import NeSST as nst +import numpy as np + def test_endf(): - ENDF_data = nst.sm.retrieve_ENDF_data('H2.json') + ENDF_data = nst.sm.retrieve_ENDF_data("H2.json") - assert ENDF_data['interactions'].total - assert ENDF_data['interactions'].elastic - assert ENDF_data['interactions'].n2n + assert ENDF_data["interactions"].total + assert ENDF_data["interactions"].elastic + assert ENDF_data["interactions"].n2n - assert ENDF_data['elastic_dxsec']['legendre'] + assert ENDF_data["elastic_dxsec"]["legendre"] - assert np.isclose(ENDF_data['n2n_dxsec']['Q_react'],-2.224640e+6) \ No newline at end of file + assert np.isclose(ENDF_data["n2n_dxsec"]["Q_react"], -2.224640e6) diff --git a/tests/test_time_of_flight.py b/tests/test_time_of_flight.py index 13e39bc..5ed2ddd 100644 --- a/tests/test_time_of_flight.py +++ b/tests/test_time_of_flight.py @@ -1,31 +1,33 @@ -import pytest -import numpy as np import NeSST as nst +import numpy as np + def test_IRF_normalisation(): """ - + For uniform sensitivity, the transform to ntof space should be conservative """ - nToF_distance = 20.0 # m + nToF_distance = 20.0 # m flat_sens = nst.time_of_flight.get_unity_sensitivity() tophat_10cmthickscintillator = nst.time_of_flight.get_transit_time_tophat_IRF(scintillator_thickness=0.1) - decayingGaussian_10cmthickscintillator = nst.time_of_flight.get_transit_time_tophat_w_decayingGaussian_IRF(scintillator_thickness=0.1,gaussian_FWHM=2.5e-9,decay_time=2e-9) + decayingGaussian_10cmthickscintillator = nst.time_of_flight.get_transit_time_tophat_w_decayingGaussian_IRF( + scintillator_thickness=0.1, gaussian_FWHM=2.5e-9, decay_time=2e-9 + ) - E_ntof = np.linspace(1.0e6,10.0e6,500) - DDmean,_,DDvar = nst.DDprimspecmoments(Tion=5.0e3) - dNdE = nst.QBrysk(E_ntof,DDmean,DDvar) + E_ntof = np.linspace(1.0e6, 10.0e6, 500) + DDmean, _, DDvar = nst.DDprimspecmoments(Tion=5.0e3) + dNdE = nst.QBrysk(E_ntof, DDmean, DDvar) - test_20m_nToF_1 = nst.time_of_flight.nToF(nToF_distance,flat_sens,tophat_10cmthickscintillator) - t_det,normt_det,signal = test_20m_nToF_1.get_signal(E_ntof,dNdE) + test_20m_nToF_1 = nst.time_of_flight.nToF(nToF_distance, flat_sens, tophat_10cmthickscintillator) + t_det, normt_det, signal = test_20m_nToF_1.get_signal(E_ntof, dNdE) - integral_signal_1 = np.trapz(y=signal,x=normt_det) + integral_signal_1 = np.trapezoid(y=signal, x=normt_det) - test_20m_nToF_2 = nst.time_of_flight.nToF(nToF_distance,flat_sens,decayingGaussian_10cmthickscintillator) - t_det,normt_det,signal = test_20m_nToF_2.get_signal(E_ntof,dNdE) + test_20m_nToF_2 = nst.time_of_flight.nToF(nToF_distance, flat_sens, decayingGaussian_10cmthickscintillator) + t_det, normt_det, signal = test_20m_nToF_2.get_signal(E_ntof, dNdE) - integral_signal_2 = np.trapz(y=signal,x=normt_det) + integral_signal_2 = np.trapezoid(y=signal, x=normt_det) - assert np.isclose(integral_signal_1,1.0,rtol=1e-3) - assert np.isclose(integral_signal_2,1.0,rtol=1e-3) + assert np.isclose(integral_signal_1, 1.0, rtol=1e-3) + assert np.isclose(integral_signal_2, 1.0, rtol=1e-3) diff --git a/tests/test_utils.py b/tests/test_utils.py index 361130f..5339c7c 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,53 +1,57 @@ -import pytest -import numpy as np import NeSST as nst +import numpy as np + def test_linear_1D_interpolate(): - x = np.array([0.0,1.0]) - y = np.array([0.0,1.0]) - f = nst.utils.interpolate_1d(x,y) + x = np.array([0.0, 1.0]) + y = np.array([0.0, 1.0]) + f = nst.utils.interpolate_1d(x, y) x_test = np.random.rand(100) y_test = f(x_test) - assert np.all(np.isclose(x_test,y_test)) + assert np.all(np.isclose(x_test, y_test)) + def test_linear_1D_interpolate_nonunique(): - x = np.array([0.0,0.0,1.0,1.0]) - y = np.array([0.0,0.0,1.0,1.0]) - f = nst.utils.interpolate_1d(x,y) + x = np.array([0.0, 0.0, 1.0, 1.0]) + y = np.array([0.0, 0.0, 1.0, 1.0]) + f = nst.utils.interpolate_1d(x, y) x_test = np.random.rand(100) y_test = f(x_test) - assert np.all(np.isclose(x_test,y_test)) - + assert np.all(np.isclose(x_test, y_test)) + + def test_linear_1D_interpolate_dimension_matching(): - x = np.array([0.0,1.0]) - y = np.array([0.0,1.0]) - f = nst.utils.interpolate_1d(x,y) - x_test = np.random.rand(100).reshape(50,2) + x = np.array([0.0, 1.0]) + y = np.array([0.0, 1.0]) + f = nst.utils.interpolate_1d(x, y) + x_test = np.random.rand(100).reshape(50, 2) y_test = f(x_test) - assert np.all(np.isclose(x_test,y_test)) + assert np.all(np.isclose(x_test, y_test)) assert x_test.shape == y_test.shape + def test_linear_2D_interpolate(): - def z_func(x,y): - return x-0.5*y - x = np.array([0.0,1.0]) - y = np.array([0.0,1.0]) - xx,yy = np.meshgrid(x,y, indexing='ij') - z = z_func(xx,yy) - f = nst.utils.interpolate_2d(x,y,z) + def z_func(x, y): + return x - 0.5 * y + + x = np.array([0.0, 1.0]) + y = np.array([0.0, 1.0]) + xx, yy = np.meshgrid(x, y, indexing="ij") + z = z_func(xx, yy) + f = nst.utils.interpolate_2d(x, y, z) # Equal sizes x_test = np.random.rand(50) y_test = np.random.rand(50) - z_test = f(x_test,y_test) - xx,yy = np.meshgrid(x_test,y_test, indexing='ij') - zz = z_func(xx,yy) + z_test = f(x_test, y_test) + xx, yy = np.meshgrid(x_test, y_test, indexing="ij") + zz = z_func(xx, yy) assert zz.shape == z_test.shape - assert np.all(np.isclose(zz,z_test)) + assert np.all(np.isclose(zz, z_test)) # Unequal sizes x_test = np.random.rand(50) y_test = np.random.rand(3) - z_test = f(x_test,y_test) - xx,yy = np.meshgrid(x_test,y_test, indexing='ij') - zz = z_func(xx,yy) - assert np.all(np.isclose(zz,z_test)) - assert zz.shape == z_test.shape \ No newline at end of file + z_test = f(x_test, y_test) + xx, yy = np.meshgrid(x_test, y_test, indexing="ij") + zz = z_func(xx, yy) + assert np.all(np.isclose(zz, z_test)) + assert zz.shape == z_test.shape diff --git a/tools/table_build.py b/tools/table_build.py index 3a97775..baa1c4e 100644 --- a/tools/table_build.py +++ b/tools/table_build.py @@ -1,10 +1,10 @@ -import numpy as np import NeSST as nst +import numpy as np -Ein = np.linspace(13.25,14.75,100) -Eout = np.linspace(1.0,15.0,400) -vvec = np.linspace(-1200.0e3,1200.0e3,80) -P1 = np.linspace(-1.0,1.0,20) +Ein = np.linspace(13.25, 14.75, 100) +Eout = np.linspace(1.0, 15.0, 400) +vvec = np.linspace(-1200.0e3, 1200.0e3, 80) +P1 = np.linspace(-1.0, 1.0, 20) -nT_table = nst.dsigdE_table("../dsigdE_tables/initial_test/","nT") -nT_table.matrix_create_and_save(Ein,Eout,vvec,P1) \ No newline at end of file +nT_table = nst.dsigdE_table("../dsigdE_tables/initial_test/", "nT") +nT_table.matrix_create_and_save(Ein, Eout, vvec, P1) From 89d475b05382baff4e37bd9043d5b237927ee391 Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Thu, 24 Jul 2025 12:25:15 +0100 Subject: [PATCH 2/5] Adding linting workflow --- .github/workflows/pre-commit.yaml | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 .github/workflows/pre-commit.yaml diff --git a/.github/workflows/pre-commit.yaml b/.github/workflows/pre-commit.yaml new file mode 100644 index 0000000..869e5f4 --- /dev/null +++ b/.github/workflows/pre-commit.yaml @@ -0,0 +1,19 @@ +name: Code linting + +on: + pull_request: + + push: + branches: + - main + +jobs: + pre-commit: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version-file: "pyproject.toml" + - uses: pre-commit/action@v3.0.1 From 86c7e82fa85f53bbefdbe46f6727f2bb45ed5918 Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Thu, 24 Jul 2025 12:27:53 +0100 Subject: [PATCH 3/5] Missing numpy typing --- src/NeSST/cross_sections.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/NeSST/cross_sections.py b/src/NeSST/cross_sections.py index d0938b9..746e438 100644 --- a/src/NeSST/cross_sections.py +++ b/src/NeSST/cross_sections.py @@ -2,6 +2,7 @@ from dataclasses import dataclass import numpy as np +import numpy.typing as npt from numpy.polynomial.legendre import legval from scipy.interpolate import griddata From 59fa744720f973dd820a0de2c280a4f249954a11 Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Thu, 24 Jul 2025 12:35:22 +0100 Subject: [PATCH 4/5] Fixing imports --- src/NeSST/__init__.py | 2 ++ src/NeSST/fitting.py | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index f1002f4..8b92b01 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,3 +1,5 @@ ################## from .constants import * from .core import * +from .fitting import * +from .time_of_flight import * \ No newline at end of file diff --git a/src/NeSST/fitting.py b/src/NeSST/fitting.py index 1f08967..645400c 100644 --- a/src/NeSST/fitting.py +++ b/src/NeSST/fitting.py @@ -1,4 +1,3 @@ -from NeSST.constants import * from NeSST.core import * from NeSST.utils import * From cf77dab7eab6cca508e4ad4abc7359c0ca550bfd Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Thu, 24 Jul 2025 12:36:17 +0100 Subject: [PATCH 5/5] linting --- src/NeSST/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index 8b92b01..69621e8 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -2,4 +2,4 @@ from .constants import * from .core import * from .fitting import * -from .time_of_flight import * \ No newline at end of file +from .time_of_flight import *