Gamma ray propagation

docs/index.rst

@@ -78,6 +78,7 @@ Mission Statement
     physics/montecarlo/index
     physics/est_and_conv/index
     physics/spectrum/index
+    physics/energy_input/index
 
 
 .. toctree::

docs/physics/energy_input/gammaray_deposition.ipynb

@@ -0,0 +1,306 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Gamma-ray energy deposition\n",
+    "\n",
+    "This notebook provides the initial implementation of Gamma-ray energy deposition into an arbitrary ejecta.\n",
+    "It is a WORK IN PROGRESS and should NOT be used for any science work until further notice."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Main loop\n",
+    "\n",
+    "Generates a simple 1D ejecta and a list of gamma-ray objects.\n",
+    "\n",
+    "Runs packets of gamma-rays through the ejecta. Handles interactions by calling the appropriate function. \n",
+    "\n",
+    "Adds deposited energy and output energy to 2 different dataframes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false,
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import astropy as apy\n",
+    "\n",
+    "from tardis.energy_input.indivisible_packets import main_gamma_ray_loop\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "from tardis.model import Radial1DModel\n",
+    "from tardis.io.config_reader import Configuration\n",
+    "\n",
+    "# Set up packet count\n",
+    "num_packets = 50000\n",
+    "\n",
+    "# Lock seed\n",
+    "np.random.seed(1)\n",
+    "\n",
+    "# Adjust model\n",
+    "config = Configuration.from_yaml(\"../../tardis/io/tests/data/tardis_configv1_density_exponential_nebular.yml\")\n",
+    "config.model.structure.velocity.start = 1 * apy.units.km / apy.units.s\n",
+    "config.model.structure.density.rho_0 = 5e2 * apy.units.g / (apy.units.cm ** 3)\n",
+    "config.supernova.time_explosion = 2.0 * apy.units.d\n",
+    "\n",
+    "config.atom_data = \"kurucz_cd23_chianti_H_He.h5\"\n",
+    "\n",
+    "# Create model\n",
+    "model = Radial1DModel.from_config(config)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate plasma"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "from tardis.plasma.properties import Density, Abundance, IsotopeAbundance, IsotopeNumberDensity, AtomicData, AtomicMass, IsotopeMass, NumberDensity, SelectedAtoms\n",
+    "from tardis.plasma.base import BasePlasma\n",
+    "from tardis.io.atom_data import AtomData\n",
+    "\n",
+    "input = [Density, Abundance, IsotopeAbundance, AtomicData, AtomicMass, IsotopeNumberDensity, NumberDensity, SelectedAtoms, IsotopeMass]\n",
+    "\n",
+    "plasma = BasePlasma(\n",
+    "        plasma_properties=input,\n",
+    "        density = model.density,\n",
+    "        abundance = model.abundance,\n",
+    "        isotope_abundance = model.raw_isotope_abundance,\n",
+    "        atomic_data = AtomData.from_hdf(config.atom_data)\n",
+    "    )\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "# Compute energy deposition rate\n",
+    "# ejecta_energy_df is the deposited energy\n",
+    "# ejecta_plot_energy_df is information for plotting\n",
+    "# escape_energy is the escaping energy\n",
+    "# decayed_packet_count is the number of packets created per shell\n",
+    "# energy_plot_positrons is the deposited energy from positrons\n",
+    "# estimated_deposition is the deposited energy from the Kasen (2006) estimator (currently not functional)\n",
+    "(\n",
+    "    energy_df,\n",
+    "    energy_plot_df,\n",
+    "    escape_energy,\n",
+    "    decayed_packet_count,\n",
+    "    energy_plot_positrons,\n",
+    "    estimated_deposition\n",
+    ") = main_gamma_ray_loop(\n",
+    "    num_packets,\n",
+    "    model,\n",
+    "    plasma,\n",
+    "    time_steps=50,\n",
+    "    time_end=50.0,\n",
+    ")\n",
+    "\n",
+    "ejecta_energy = energy_plot_df[\"energy_input\"]\n",
+    "ejecta_energy_r = energy_plot_df[\"energy_input_r\"]\n",
+    "energy_input_time = energy_plot_df[\"energy_input_time\"]\n",
+    "energy_input_type = energy_plot_df[\"energy_input_type\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Plotting results\n",
+    "\n",
+    "Energy deposited at a given radius"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(dpi=150, facecolor='w')\n",
+    "ax = fig.add_subplot(111)\n",
+    "\n",
+    "scatter = ax.scatter(np.array(ejecta_energy_r)/np.max(model.v_outer.value), np.array(ejecta_energy), s=1, alpha=0.1)\n",
+    "ax.set_xlabel(\"r/R\")\n",
+    "ax.set_ylabel(\"E (keV)\")\n",
+    "ax.semilogy();"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Interactions binned by radius"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(dpi=150, facecolor='w')\n",
+    "ax = fig.add_subplot(111)\n",
+    "ax.hist(np.array(ejecta_energy_r), bins=200)\n",
+    "#ax.set_xlim(0, 1)\n",
+    "ax.set_xlabel(\"r (cm)\")\n",
+    "ax.set_ylabel(\"Interaction count\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Density Profile"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(dpi=150, facecolor='w')\n",
+    "plt.semilogy(model.r_middle/np.max(model.r_outer), model.density)\n",
+    "plt.plot(0,0)\n",
+    "plt.ylabel(\"Density g cm$^{-3}$\")\n",
+    "plt.xlabel(\"r/R\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fraction of energy escaping from the ejecta"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.sum(np.sum(escape_energy) / (np.sum(escape_energy) + np.sum(energy_df)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Spectrum of escape energy at the final time step"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from tardis.energy_input.energy_source import read_artis_lines\n",
+    "\n",
+    "ni56_lines = read_artis_lines(\"ni56\")\n",
+    "co56_lines = read_artis_lines(\"co56\")\n",
+    "\n",
+    "plt.figure(figsize=(12, 6), dpi=150)\n",
+    "plt.step(escape_energy.index, escape_energy.iloc[:,49], label=\"$\\gamma$-spectrum\", where=\"post\")\n",
+    "plt.xlabel(\"Energy (keV)\")\n",
+    "plt.ylabel(\"keV/keV/s\");\n",
+    "plt.loglog()\n",
+    "plt.xlim(100, 4000)\n",
+    "#plt.ylim(0.01, 100)\n",
+    "plt.vlines(ni56_lines.energy*1000, 0.01, 0.02, color=\"k\", label=\"Ni56\")\n",
+    "plt.vlines(co56_lines.energy*1000, 0.01, 0.02, color=\"r\", label=\"Co56\")\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Energy deposition rate\n",
+    "\n",
+    "Dataframe index is the radial grid location. Columns are time steps in seconds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "energy_df"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Energy deposition rate versus radius at time_explosion"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(dpi=150, facecolor='w')\n",
+    "plt.semilogy(model.v_middle/np.max(model.v_outer), energy_df.iloc[:, 0])\n",
+    "plt.xlabel(\"r/R\")\n",
+    "plt.ylabel(\"Energy deposition rate [eV/s$]\");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "1d9f0af959d2c458edd906b0042305fe807dd53effdb40eb4325690d1852d0c6"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.13 ('tardis_new')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

docs/physics/energy_input/index.rst

@@ -0,0 +1,6 @@
+***************************
+Gamma Ray Energy Deposition
+***************************
+
+.. toctree::
+    gammaray_deposition

tardis/energy_input/GXPacket.py

@@ -0,0 +1,70 @@
+from enum import IntEnum
+import numpy as np
+from numba import int64, float64
+from numba.experimental import jitclass
+
+
+class GXPacketStatus(IntEnum):
+    BETA_DECAY = -1
+    COMPTON_SCATTER = 0
+    PHOTOABSORPTION = 1
+    PAIR_CREATION = 2
+    IN_PROCESS = 3
+    END = 4
+
+
+gxpacket_spec = [
+    ("location", float64[:]),
+    ("direction", float64[:]),
+    ("energy_rf", float64),
+    ("energy_cmf", float64),
+    ("nu_rf", float64),
+    ("nu_cmf", float64),
+    ("status", int64),
+    ("shell", int64),
+    ("time_current", float64),
+    ("tau", float64),
+]
+
+
+@jitclass(gxpacket_spec)
+class GXPacket(object):
+    """
+    Indivisible gamma-ray packet
+    """
+
+    def __init__(
+        self,
+        location,
+        direction,
+        energy_rf,
+        energy_cmf,
+        nu_rf,
+        nu_cmf,
+        status,
+        shell,
+        time_current,
+    ):
+        self.location = location
+        self.direction = direction
+        self.energy_rf = energy_rf
+        self.energy_cmf = energy_cmf
+        self.nu_rf = nu_rf
+        self.nu_cmf = nu_cmf
+        self.status = status
+        self.shell = shell
+        self.time_current = time_current
+        # TODO: rename to tau_event
+        self.tau = -np.log(np.random.random())
+
+    def get_location_r(self):
+        """Calculate radius of the packet
+
+        Returns:
+            float: packet radius
+        """
+        return np.sqrt(
+            self.location[0] ** 2.0
+            + self.location[1] ** 2.0
+            + self.location[2] ** 2.0
+        )