add a constants module (#688)

all modules now get physical constants from this instead of scipy directly.

Co-authored-by: Eric T. Johnson <yut23@users.noreply.github.com>

pynucastro/constants/__init__.py

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

pynucastro/constants/constants.py

@@ -0,0 +1,25 @@
+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")
+m_p = physical_constants.value("proton mass") * physical_constants.kilo
+m_p_MeV = physical_constants.value("proton mass energy equivalent in MeV")
+m_n = physical_constants.value("neutron mass") * physical_constants.kilo
+m_n_MeV = physical_constants.value("neutron mass energy equivalent in MeV")
+m_e = physical_constants.value("electron mass") * physical_constants.kilo
+m_e_MeV = physical_constants.value("electron mass energy equivalent in MeV")
+
+N_A = physical_constants.value("Avogadro constant")
+
+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
+
+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
+
+erg2MeV = physical_constants.erg / (physical_constants.eV * physical_constants.mega)
+MeV2erg = (physical_constants.eV * physical_constants.mega) / physical_constants.erg

pynucastro/constants/test/test_constants.py

@@ -0,0 +1,32 @@
+import math
+
+from pytest import approx
+
+from pynucastro.constants import constants
+
+
+class TestConstants:
+    def test_constants(self):
+
+        assert constants.m_u == approx(1.66054e-24, abs=0, rel=1.e-6)
+        assert constants.m_u_MeV == approx(constants.m_u * constants.c_light**2 * constants.erg2MeV, abs=0, rel=1.e-6)
+        assert constants.m_p == approx(1.67262192e-24, abs=0, rel=1.e-6)
+        assert constants.m_p_MeV == approx(constants.m_p * constants.c_light**2 * constants.erg2MeV, abs=0, rel=1.e-6)
+        assert constants.m_n == approx(1.674927498e-24, abs=0, rel=1.e-6)
+        assert constants.m_n_MeV == approx(constants.m_n * constants.c_light**2 * constants.erg2MeV, abs=0, rel=1.e-6)
+        assert constants.m_e == approx(9.1093837e-28, abs=0, rel=1.e-6)
+        assert constants.m_e_MeV == approx(constants.m_e * constants.c_light**2 * constants.erg2MeV)
+
+        assert constants.N_A == approx(6.02214e23)
+
+        assert constants.k == approx(1.38064e-16, abs=0, rel=1.e-5)
+        assert constants.k_MeV == approx(constants.k * constants.erg2MeV, abs=0, rel=1.e-5)
+
+        assert constants.h == approx(6.62607e-27, abs=0, rel=1.e-6)
+        assert constants.hbar == approx(6.62607e-27 / 2 / math.pi, abs=0, rel=1.e-6)
+
+        assert constants.q_e == approx(4.803e-10, abs=0, rel=1.e-4)
+        assert constants.c_light == approx(2.997924e10)
+
+        assert constants.erg2MeV == approx(1.0 / 1.60218e-6, abs=0, rel=1.e-4)
+        assert constants.MeV2erg == approx(1.60218e-6, abs=0, rel=1.e-4)

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/networks/rate_collection.py

@@ -18,9 +18,9 @@
 from matplotlib.scale import SymmetricalLogTransform
 from matplotlib.ticker import MaxNLocator
 from mpl_toolkits.axes_grid1 import make_axes_locatable
-from scipy import constants
 
 # Import Rate
+from pynucastro.constants import constants
 from pynucastro.nucdata import Nucleus, PeriodicTable
 from pynucastro.rates import (ApproximateRate, DerivedRate, Library, Rate,
                               RateFileError, RatePair, TabularRate,
@@ -1119,19 +1119,17 @@ def evaluate_energy_generation(self, rho, T, composition,
         enuc = 0.
 
         # compute constants and units
-        m_n_MeV = constants.value('neutron mass energy equivalent in MeV')
-        m_p_MeV = constants.value('proton mass energy equivalent in MeV')
-        m_e_MeV = constants.value('electron mass energy equivalent in MeV')
-        MeV2erg = (constants.eV * constants.mega) / constants.erg
 
         # ion binding energy contributions. basically e=mc^2
         for nuc in self.unique_nuclei:
             # add up mass in MeV then convert to erg
-            mass = ((nuc.A - nuc.Z) * m_n_MeV + nuc.Z * (m_p_MeV + m_e_MeV) - nuc.A * nuc.nucbind) * MeV2erg
+            mass = ((nuc.A - nuc.Z) * constants.m_n_MeV +
+                    nuc.Z * (constants.m_p_MeV + constants.m_e_MeV) -
+                    nuc.A * nuc.nucbind) * constants.MeV2erg
             enuc += ydots[nuc] * mass
 
         # convert from molar value to erg/g/s
-        enuc *= -1*constants.Avogadro
+        enuc *= -1*constants.N_A
 
         # subtract neutrino losses for tabular weak reactions
         enu = 0.0
@@ -1143,7 +1141,7 @@ def evaluate_energy_generation(self, rho, T, composition,
 
                 # need to get reactant nucleus
                 nuc = r.reactants[0]
-                enu += constants.Avogadro * ys[nuc] * r.get_nu_loss(T, rho * y_e)
+                enu += constants.N_A * ys[nuc] * r.get_nu_loss(T, rho * y_e)
 
         enuc -= enu
         if return_enu:

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] = n.A_nuc * constants.m_u
 
     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):