add neutrino cooling routines

docs/source/index.rst

@@ -42,16 +42,14 @@ codes.
 
 .. toctree::
    :maxdepth: 1
-   :caption: Advanced Usage
+   :caption: Advanced Rate Operations
 
    screening-examples
    modify-example.ipynb
    approx-rates-examples.ipynb
    custom-rates.ipynb
    unimportant-rates.ipynb
    partition-function
-   NSE-example.ipynb
-   julia
 
 .. toctree::
    :maxdepth: 1
@@ -60,6 +58,24 @@ codes.
    networks
    cxx-networks
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Nuclear Statistical Equilibrium
+
+   NSE-example.ipynb
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Plasma Neutrinos
+
+   neutrino-cooling.ipynb
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Interfacing With Other Languages
+
+   julia
+
 .. toctree::
    :maxdepth: 1
    :caption: Reference

pynucastro/neutrino_cooling/__init__.py

@@ -0,0 +1,6 @@
+"""Neutrino cooling routines"""
+
+__all__ = ["neutrino_cooling", "sneut5"]
+
+from .neutrino_cooling import NeutrinoCooling
+from .sneut5 import sneut5

pynucastro/neutrino_cooling/neutrino_cooling.py

@@ -0,0 +1,46 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from .sneut5 import sneut5
+
+
+class NeutrinoCooling:
+    """A class to provide an interface to explore neutrino cooling.
+    This calls a specific implementation of the cooling that includes
+    contributions from pairs, plasma, recombination, bremsstrahlung,
+    and photoneutrinos."""
+
+    def __init__(self, neutrino_function=sneut5):
+
+        self.func = neutrino_function
+
+    def plot(self, *, Tmin=1.e7, Tmax=1.e10, rhomin=1.e3, rhomax=1.e10,
+             abar=20, zbar=10, npts_temp=50, npts_rho=50):
+        """given a fixed abar/zbar, make a plot of neutrino cooling over
+        a range of density and temperature"""
+
+        Ts = np.logspace(np.log10(Tmin), np.log10(Tmax), npts_temp)
+        rhos = np.logspace(np.log10(rhomin), np.log10(rhomax), npts_rho)
+
+        data = np.zeros((npts_temp, npts_rho))
+
+        for i, T in enumerate(Ts):
+            for j, rho in enumerate(rhos):
+                data[i, j] = np.log10(self.func(rho, T, abar=abar, zbar=zbar))
+
+        fig, ax = plt.subplots()
+
+        # by default, imshow puts the first index on the y axis
+        # and the second on the x axis, so label accordingly
+
+        extent = (np.log10(rhomin), np.log10(rhomax),
+                  np.log10(Tmin), np.log10(Tmax))
+
+        im = ax.imshow(data, cmap="viridis",
+                       extent=extent, aspect="auto", origin="lower")
+        fig.colorbar(im, ax=ax)
+
+        ax.set_ylabel(r"$\log(T)$ [K]")
+        ax.set_xlabel(r"$\log(\rho)$ [g/cm$^3$]")
+
+        return fig