add a constants module

pynucastro/constants/__init__.py

@@ -0,0 +1,3 @@
+"""Fundamental constants in CGS (unless otherwise noted)"""
+
+__all__ = ["constants"]

pynucastro/constants/constants.py

@@ -0,0 +1,17 @@
+import scipy.constants as physical_constants
+
+# atomic unit mass in g
+m_u = physical_constants.value("unified atomic mass unit") * physical_constants.kilo
+m_u_MeV = physical_constants.value('atomic mass constant energy equivalent in MeV')
+k = physical_constants.value("Boltzmann constant") / physical_constants.erg
+k_MeV = physical_constants.value("Boltzmann constant in eV/K") / physical_constants.mega
+h = physical_constants.value("Planck constant") / physical_constants.erg
+hbar = physical_constants.value("reduced Planck constant") / physical_constants.erg
+N_A = physical_constants.value("Avogadro constant")
+q_e = physical_constants.value("elementary charge") * (physical_constants.c * 100) / 10  # C to statC (esu)
+c_light = physical_constants.value("speed of light in vacuum") / physical_constants.centi
+m_p = physical_constants.value("proton mass") * physical_constants.kilo
+m_n = physical_constants.value("neutron mass") * physical_constants.kilo
+
+erg2MeV = physical_constants.erg / (physical_constants.eV * physical_constants.mega)
+MeV2erg = (physical_constants.eV * physical_constants.mega) / physical_constants.erg

pynucastro/networks/nse_network.py

@@ -1,10 +1,10 @@
 import warnings
 
 import numpy as np
-from scipy import constants
 from scipy.optimize import fsolve
 
 from pynucastro._version import version
+from pynucastro.constants import constants
 from pynucastro.networks.rate_collection import Composition, RateCollection
 from pynucastro.nucdata import Nucleus
 from pynucastro.rates import TabularRate
@@ -75,11 +75,9 @@ class NSENetwork(RateCollection):
 
     def _evaluate_n_e(self, state, Xs):
 
-        m_u = constants.value("unified atomic mass unit") * constants.kilo  # atomic unit mass in g
-
         n_e = 0.0
         for nuc in self.unique_nuclei:
-            n_e += nuc.Z * state.dens * Xs[nuc] / (nuc.A * m_u)
+            n_e += nuc.Z * state.dens * Xs[nuc] / (nuc.A * constants.m_u)
 
         return n_e
 
@@ -96,10 +94,6 @@ def _evaluate_mu_c(self, n_e, state, use_coulomb_corr=True):
         ye_max = max(nuc.Z/nuc.A for nuc in self.unique_nuclei)
         assert state.ye >= ye_low and state.ye <= ye_max, "input electron fraction goes outside of scope for current network"
 
-        # Setting up the constants needed to compute mu_c
-        k = constants.value("Boltzmann constant") / constants.erg           # boltzmann in erg/K
-        Erg2MeV = constants.erg / (constants.eV * constants.mega)
-
         # u_c is the coulomb correction term for NSE
         # Calculate the composition at NSE, equations found in appendix of Calder paper
         u_c = {}
@@ -112,7 +106,7 @@ def _evaluate_mu_c(self, n_e, state, use_coulomb_corr=True):
                 A_3 = -0.5 * np.sqrt(3.0) - A_1 / np.sqrt(A_2)
 
                 Gamma = state.gamma_e_fac * n_e ** (1.0/3.0) * nuc.Z ** (5.0 / 3.0) / state.temp
-                u_c[nuc] = Erg2MeV * k * state.temp * (A_1 * (np.sqrt(Gamma * (A_2 + Gamma)) - A_2 * np.log(np.sqrt(Gamma / A_2) +
+                u_c[nuc] = constants.erg2MeV * constants.k * state.temp * (A_1 * (np.sqrt(Gamma * (A_2 + Gamma)) - A_2 * np.log(np.sqrt(Gamma / A_2) +
                                       np.sqrt(1.0 + Gamma / A_2))) + 2.0 * A_3 * (np.sqrt(Gamma) - np.arctan(np.sqrt(Gamma))))
             else:
                 u_c[nuc] = 0.0
@@ -121,12 +115,6 @@ def _evaluate_mu_c(self, n_e, state, use_coulomb_corr=True):
 
     def _nucleon_fraction_nse(self, u, u_c, state):
 
-        # Define constants: amu, boltzmann, planck, and electron charge
-        m_u = constants.value("unified atomic mass unit") * constants.kilo  # atomic unit mass in g
-        k = constants.value("Boltzmann constant") / constants.erg           # boltzmann in erg/K
-        h = constants.value("Planck constant") / constants.erg              # in cgs
-        Erg2MeV = constants.erg / (constants.eV * constants.mega)
-
         Xs = {}
         up_c = u_c[Nucleus("p")]
 
@@ -139,10 +127,10 @@ def _nucleon_fraction_nse(self, u, u_c, state):
             if not nuc.spin_states:
                 raise ValueError(f"The spin of {nuc} is not implemented for now.")
 
-            nse_exponent = (nuc.Z * u[0] + nuc.N * u[1] - u_c[nuc] + nuc.Z * up_c + nuc.nucbind * nuc.A) / (k * state.temp * Erg2MeV)
+            nse_exponent = (nuc.Z * u[0] + nuc.N * u[1] - u_c[nuc] + nuc.Z * up_c + nuc.nucbind * nuc.A) / (constants.k * state.temp * constants.erg2MeV)
             nse_exponent = min(500.0, nse_exponent)
 
-            Xs[nuc] = m_u * nuc.A_nuc * pf * nuc.spin_states / state.dens * (2.0 * np.pi * m_u * nuc.A_nuc * k * state.temp / h**2) ** 1.5 \
+            Xs[nuc] = constants.m_u * nuc.A_nuc * pf * nuc.spin_states / state.dens * (2.0 * np.pi * constants.m_u * nuc.A_nuc * constants.k * state.temp / constants.h**2) ** 1.5 \
                     * np.exp(nse_exponent)
 
         return Xs

pynucastro/networks/python_network.py

@@ -5,8 +5,7 @@
 import shutil
 import sys
 
-from scipy import constants
-
+from pynucastro.constants import constants
 from pynucastro.networks.rate_collection import RateCollection
 from pynucastro.rates.rate import ApproximateRate
 from pynucastro.screening import get_screening_map
@@ -204,15 +203,12 @@ def _write_network(self, outfile=None):
 
         # we'll compute the masses here in erg
 
-        m_u = constants.value('atomic mass constant energy equivalent in MeV')
-        MeV2erg = (constants.eV * constants.mega) / constants.erg
-
         of.write("\n")
 
         of.write("# masses in ergs\n")
         of.write("mass = np.zeros((nnuc), dtype=np.float64)\n\n")
         for n in self.unique_nuclei:
-            mass = n.A_nuc * m_u * MeV2erg
+            mass = n.A_nuc * constants.m_u_MeV * constants.MeV2erg
             of.write(f"mass[j{n.raw}] = {mass}\n")
 
         of.write("\n")

pynucastro/nucdata/nucleus.py

@@ -5,8 +5,7 @@
 import os
 import re
 
-from scipy.constants import physical_constants
-
+from pynucastro.constants import constants
 from pynucastro.nucdata.binding_table import BindingTable
 from pynucastro.nucdata.elements import PeriodicTable
 from pynucastro.nucdata.mass_table import MassTable
@@ -17,9 +16,6 @@
 _pynucastro_rates_dir = os.path.join(_pynucastro_dir, 'library')
 _pynucastro_tabular_dir = os.path.join(_pynucastro_rates_dir, 'tabular')
 
-#set the atomic mass unit constant in MeV
-m_u, _, _ = physical_constants['atomic mass constant energy equivalent in MeV']
-
 #read the mass excess table once and store it at the module-level
 _mass_table = MassTable()
 
@@ -147,7 +143,7 @@ def __init__(self, name, dummy=False):
 
         # Now we will define the Nuclear Mass,
         try:
-            self.A_nuc = float(self.A) + _mass_table.get_mass_diff(a=self.A, z=self.Z) / m_u
+            self.A_nuc = float(self.A) + _mass_table.get_mass_diff(a=self.A, z=self.Z) / constants.m_u_MeV
         except NotImplementedError:
             self.A_nuc = None
 

pynucastro/rates/rate.py

@@ -10,7 +10,6 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy.constants import physical_constants
 
 try:
     import numba
@@ -43,16 +42,9 @@ def wrapper(*args, **kwargs):
         return wrap(cls_or_spec)
 
 
+from pynucastro.constants import constants
 from pynucastro.nucdata import Nucleus, UnsupportedNucleus
 
-amu_mev, _, _ = physical_constants['atomic mass constant energy equivalent in MeV']
-hbar, _, _ = physical_constants['reduced Planck constant']
-amu, _, _ = physical_constants['atomic mass constant']
-k_B_mev_k, _, _ = physical_constants['Boltzmann constant in eV/K']
-k_B_mev_k /= 1.0e6
-k_B, _, _ = physical_constants['Boltzmann constant']
-N_a, _, _ = physical_constants['Avogadro constant']
-
 _pynucastro_dir = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
 _pynucastro_rates_dir = os.path.join(_pynucastro_dir, 'library')
 _pynucastro_tabular_dir = os.path.join(_pynucastro_rates_dir, 'tabular')
@@ -1751,7 +1743,7 @@ def __init__(self, rate, compute_Q=False, use_pf=False):
             a = ssets.a
             prefactor = 0.0
             Q = 0.0
-            prefactor += -np.log(N_a) * (len(self.rate.reactants) - len(self.rate.products))
+            prefactor += -np.log(constants.N_A) * (len(self.rate.reactants) - len(self.rate.products))
 
             for nucr in self.rate.reactants:
                 prefactor += 1.5*np.log(nucr.A) + np.log(nucr.spin_states)
@@ -1761,7 +1753,7 @@ def __init__(self, rate, compute_Q=False, use_pf=False):
                 Q -= nucp.A_nuc
 
             if self.compute_Q:
-                Q = Q * amu_mev
+                Q = Q * constants.m_u_MeV
             else:
                 Q = self.rate.Q
 
@@ -1770,12 +1762,12 @@ def __init__(self, rate, compute_Q=False, use_pf=False):
             if len(self.rate.reactants) == len(self.rate.products):
                 prefactor += 0.0
             else:
-                F = (amu * k_B * 1.0e5 / (2.0*np.pi*hbar**2))**(1.5*(len(self.rate.reactants) - len(self.rate.products)))
+                F = (constants.m_u * constants.k * 1.0e9 / (2.0*np.pi*constants.hbar**2))**(1.5*(len(self.rate.reactants) - len(self.rate.products)))
                 prefactor += np.log(F)
 
             a_rev = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
             a_rev[0] = prefactor + a[0]
-            a_rev[1] = a[1] - Q / (1.0e9 * k_B_mev_k)
+            a_rev[1] = a[1] - Q / (1.0e9 * constants.k_MeV)
             a_rev[2] = a[2]
             a_rev[3] = a[3]
             a_rev[4] = a[4]

pynucastro/reduction/reduction.py

@@ -8,6 +8,7 @@
 import numpy as np
 
 from pynucastro import Composition, Nucleus
+from pynucastro.constants import constants
 from pynucastro.reduction import drgep, mpi_importer, sens_analysis
 from pynucastro.reduction.generate_data import dataset
 from pynucastro.reduction.load_network import load_network
@@ -45,29 +46,22 @@ def get_net_info(net, comp, rho, T):
     ebind = np.zeros(len(net.unique_nuclei), dtype=np.float64)
     m = np.zeros(len(net.unique_nuclei), dtype=np.float64)
 
-    mass_neutron = 1.67492721184e-24
-    mass_proton = 1.67262163783e-24
-    c_light = 2.99792458e10
-
     for i, n in enumerate(net.unique_nuclei):
 
         y[i] = y_dict[n]
         ydot[i] = ydots_dict[n]
         z[i] = n.Z
         a[i] = n.A
         ebind[i] = n.nucbind or 0.0
-        m[i] = mass_proton * n.Z + mass_neutron * n.N - ebind[i] / c_light**2
+        m[i] = constants.m_p * n.Z + constants.m_n * n.N - ebind[i] / constants.c_light**2
 
     return NetInfo(y, ydot, z, a, ebind, m)
 
 
 def enuc_dot(net_info):
     """Calculate the nuclear energy generation rate."""
 
-    avo = 6.0221417930e23
-    c_light = 2.99792458e10
-
-    return -np.sum(net_info.ydot * net_info.m) * avo * c_light**2
+    return -np.sum(net_info.ydot * net_info.m) * constants.N_A * constants.c_light**2
 
 
 def ye_dot(net_info):

pynucastro/screening/screen.py

@@ -2,8 +2,8 @@
 Python implementations of screening routines.
 """
 import numpy as np
-from scipy import constants
 
+from pynucastro.constants import constants
 # use the jitclass placeholder from rate.py
 from pynucastro.rates.rate import jitclass, numba
 
@@ -17,12 +17,6 @@ def njit(func):
            "make_plasma_state", "make_screen_factors", "potekhin_1998"]
 
 
-amu = constants.value("atomic mass constant") / constants.gram  # kg to g
-q_e = constants.value("elementary charge") * (constants.c * 100) / 10  # C to statC (esu)
-hbar = constants.value("reduced Planck constant") / constants.erg  # J*s to erg*s
-k_B = constants.value("Boltzmann constant") / constants.erg  # J/K to erg/K
-
-
 @jitclass()
 class PlasmaState:
     """
@@ -62,14 +56,14 @@ def __init__(self, temp, dens, Ys, Zs):
         self.z2bar = np.sum(Zs ** 2 * Ys) / ytot
 
         # Average mass and total number density
-        mbar = self.abar * amu
+        mbar = self.abar * constants.m_u
         ntot = self.dens / mbar
         # Electron number density
         # zbar * ntot works out to sum(z[i] * n[i]), after cancelling terms
         self.n_e = self.zbar * ntot
 
         # temperature-independent part of Gamma_e, from Chugunov 2009 eq. 6
-        self.gamma_e_fac = q_e ** 2 / k_B * np.cbrt(4 * np.pi / 3) * np.cbrt(self.n_e)
+        self.gamma_e_fac = constants.q_e ** 2 / constants.k * np.cbrt(4 * np.pi / 3) * np.cbrt(self.n_e)
 
 
 @jitclass()
@@ -107,7 +101,7 @@ def __init__(self, temp, dens, ye):
         self.temp = temp
         self.dens = dens
         self.ye = ye
-        self.gamma_e_fac = q_e ** 2 / k_B * np.cbrt(4.0 * np.pi / 3.0)
+        self.gamma_e_fac = constants.q_e ** 2 / constants.k * np.cbrt(4.0 * np.pi / 3.0)
 
 
 def make_plasma_state(temp, dens, molar_fractions):
@@ -231,9 +225,9 @@ def chugunov_2007(state, scn_fac):
     mu12 = scn_fac.a1 * scn_fac.a2 / (scn_fac.a1 + scn_fac.a2)
     z_factor = scn_fac.z1 * scn_fac.z2
     n_i = state.n_e / scn_fac.ztilde ** 3
-    m_i = 2 * mu12 * amu
+    m_i = 2 * mu12 * constants.m_u
 
-    T_p = hbar / k_B * q_e * np.sqrt(4 * np.pi * z_factor * n_i / m_i)
+    T_p = constants.hbar / constants.k * constants.q_e * np.sqrt(4 * np.pi * z_factor * n_i / m_i)
 
     # Normalized temperature
     T_norm = state.temp / T_p
@@ -350,7 +344,7 @@ def chugunov_2009(state, scn_fac):
     Gamma_12 = Gamma_e * z1z2 / scn_fac.ztilde
 
     # Coulomb barrier penetrability, eq. 10
-    tau_factor = np.cbrt(27 / 2 * (np.pi * q_e ** 2 / hbar) ** 2 * amu / k_B)
+    tau_factor = np.cbrt(27 / 2 * (np.pi * constants.q_e ** 2 / constants.hbar) ** 2 * constants.m_u / constants.k)
     tau_12 = tau_factor * scn_fac.aznut / np.cbrt(state.temp)
 
     # eq. 12