{ "cells": [ { "cell_type": "markdown", "id": "6c0dbe0a", "metadata": {}, "source": [ "# Energy Packet Initialization\n", "\n", "## Theory\n", "\n", "While it is instructive to think about tracking the propagation history of\n", "individual photons when illustrating the basic idea behind Monte Carlo radiative transfer\n", "techniques, there are important numerical reasons for using a different\n", "discretization scheme. Instead of thinking in the photon picture, it brings\n", "significant advantages to follow the idea of [] and\n", "[] and consider parcels of radiant energy as the fundamental\n", "building blocks of the Monte Carlo calculation. These basic Monte Carlo quanta\n", "are commonly referred to as \"energy packets\" or simply \"packets\", and are composed of many photons with the same frequency.\n", "\n", "During a Monte Carlo calculation, $N$ (a large number) packets, all with a certain\n", "energy $\\varepsilon$, are created at the inner boundary of the computational domain (which is discussed in the [model section](../setup/model.rst)) known as the photosphere. Currently, the photosphere is modeled as a spherical [blackbody](https://en.wikipedia.org/wiki/Black-body_radiation) with a radius $r_\\mathrm{boundary\\_inner}$ and temperature $T_\\mathrm{inner}$. Both of these quantities are calculated as a part of the [model](../setup/model.ipynb), and $T_\\mathrm{inner}$ is additionally updated throughout the simulation as a part of the [convergence](../est_and_conv/convergence.rst) process.\n", "\n", "In TARDIS, all packets are assigned identical energies **in the lab frame** (see [Reference Frames](propagation.rst#reference-frames)), and the total (lab-frame) energy of the packets is 1 erg (and thus each packet has an energy of $\\frac{1}{N}$ ergs).\n", "\n", "
\n", " \n", "Note\n", "\n", "The indivisible energy packet scheme does not require that all packets have the same energy. This is just a convenient and simple choice adopted in TARDIS.\n", "\n", "
\n", "\n", "
\n", " \n", "Note\n", "\n", "Since all packets have the same total, and photon energy is proportional to frequency, higher-frequency packets will represent less real photons than lower-frequency packets.\n", "\n", "
\n", "\n", "Since the photosphere is modeled as a blackbody, its total luminosity $L_\\mathrm{inner}$ (recall that luminosity is energy emitted divided by the time in which it is emitted) is\n", "$$L_\\mathrm{inner}=\\frac{N\\varepsilon}{\\Delta t}=4 \\pi r_\\mathrm{boundary\\_inner}^2 \\sigma_{\\mathrm{R}} T_{\\mathrm{inner}}^4$$\n", "where $\\sigma_\\mathrm{R}$ is the Stefan-Boltzmann constant and $\\Delta t$ is the physical duration of the simulation. In order to make this relationship hold (remembering that $N\\varepsilon = 1$ erg), we use\n", "$$\\Delta t = \\frac{1}{L_\\mathrm{inner}}=\\frac{1}{4 \\pi r_\\mathrm{boundary\\_inner}^2 \\sigma_{\\mathrm{R}} T_{\\mathrm{inner}}^4}.$$\n", "\n", "
\n", " \n", "Note\n", "\n", "As will be shown in the code example, this will lead to unphysically small values for $\\Delta t$. It may be easier to think of the Monte Carlo packets not as packets of energy $\\epsilon$ going through a simulation of duration $\\Delta t$, but as packets of luminosity that carry an energy $\\epsilon$ over a time $\\Delta t$ (and thus truly being luminosity packets of luminosity $\\frac{\\epsilon}{\\Delta t}$). Indeed, this view of the packets will be useful when deriving the [Monte Carlo Estimators](../est_and_conv/estimators.ipynb).\n", "\n", "
\n", "\n", "During packet initialization, each packet is assigned an initial propagation direction $\\mu$ which is the cosine of the angle $\\theta$ which the packet's path makes with the radial direction (see the image below). Using a pseudo-random number generator which generates numbers $z$ uniformly distributed on the interval $[0,1]$, the propagation direction is determined (due to physical considerations beyond the scope of this documentation) according to\n", "$$\\mu = \\sqrt{z}.$$\n", "This sampling is shown in the \"Code Example\" section.\n", "\n", "\n", "\n", "Finally, each packet is assigned an initial frequency (or more precisely, the initial frequency of its constituent photons). Note that since each packet has the same energy, each packet will represent a different number of real photons. The sampling on packet frequencies is more involved than that of the propagation direction, as it involves sampling the Planck distribution (see below). TARDIS uses the technique described in [] and summarized in [] for this purpose.\n", "\n", "During the simulation, the energy of the packet remains constant in the local\n", "co-moving frame (see [Reference Frames](propagation.rst#reference-frames)). This naturally ensures energy\n", "conservation and constitutes the main advantage of this discretization scheme. **However, while the energy of the packets is conserved in the co-moving frame, the co-moving frequency of the packet (and thus the photons represented by the packet) may vary over the course of the simulation. Thus, a packet may represent several different numbers of real photons throughout their lifetimes.**\n", "\n", "## Code Example\n", "\n", "We now demonstrate the TARDIS packet initialization framework:" ] }, { "cell_type": "code", "execution_count": null, "id": "426325e5", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from tardis.montecarlo.packet_source import BlackBodySimpleSource\n", "from astropy import units as u\n", "from tardis import constants as const\n", "import matplotlib.pyplot as plt\n", "\n", "#The random number generator that will be used\n", "rng = np.random.default_rng()" ] }, { "cell_type": "markdown", "id": "4ae02998", "metadata": {}, "source": [ "The following cell contains values that you can change to see how it affects the spectrum (in an actual simulation, the seed and number of packets are set in the [Monte Carlo configuration](../../io/configuration/components/montecarlo.rst), and the photospheric radius is calculated as a part of the [model](../setup/model.ipynb)):" ] }, { "cell_type": "code", "execution_count": null, "id": "bc34bf33", "metadata": {}, "outputs": [], "source": [ "# Seed for the pseudo-random number generator\n", "seed = 1\n", "\n", "# Number of packets generated\n", "n_packets = 40000\n", "\n", "# Radius of the supernova's photosphere in cm\n", "r_boundary_inner = 1e15 * u.cm" ] }, { "cell_type": "markdown", "id": "450faf76", "metadata": {}, "source": [ "We set the temperature of the photosphere $T_\\mathrm{inner}$, which will determine the photospheric luminosity (in an actual simulation, $T_\\mathrm{inner}$ is initially calculated as a part of the [model](../setup/model.ipynb) and updated as a part of the [convergence](../est_and_conv/convergence.rst) process):" ] }, { "cell_type": "code", "execution_count": null, "id": "3fb3ca8c", "metadata": {}, "outputs": [], "source": [ "# Temperature in K\n", "temperature_inner = 10000 * u.K\n", "\n", "luminosity_inner = 4 * np.pi * (r_boundary_inner**2) * const.sigma_sb * (temperature_inner**4)\n", "\n", "# Makes sure the luminosity is given in erg/s\n", "luminosity_inner = luminosity_inner.to('erg/s')\n", "\n", "print('Luminosity of photosphere:', luminosity_inner)" ] }, { "cell_type": "markdown", "id": "516633c5", "metadata": {}, "source": [ "We now generate the ensemble of packets. The array of packet energies and radii are also shown." ] }, { "cell_type": "code", "execution_count": null, "id": "925e9e1b", "metadata": {}, "outputs": [], "source": [ "# We define our packet source\n", "packet_source = BlackBodySimpleSource(seed)\n", "\n", "radii, nus, mus, energies = packet_source.create_packets(\n", " temperature_inner.value, \n", " n_packets, \n", " rng, \n", " r_boundary_inner)\n", "\n", "# Sets the energies in units of ergs\n", "energies *= u.erg\n", "\n", "# Sets the frequencies in units of Hz\n", "nus *= u.Hz\n", "\n", "print('Energies:', energies)\n", "print('Radii:', radii)" ] }, { "cell_type": "markdown", "id": "16936bce", "metadata": {}, "source": [ "We set the timespan of the simulation so that each packet contributes the appropriate luminosity to the spectrum." ] }, { "cell_type": "code", "execution_count": null, "id": "fed35f47", "metadata": {}, "outputs": [], "source": [ "# Time of simulation\n", "t_simulation = 1 * u.erg / luminosity_inner\n", "print('Time of simulation:', t_simulation)\n", "\n", "# Array of luminosity contribution by each packet\n", "lumin_per_packet = energies / t_simulation\n", "print('Luminosity per packet:', lumin_per_packet)" ] }, { "cell_type": "markdown", "id": "0d839222", "metadata": {}, "source": [ "We define important constants, and for comparison's sake, we code the Planck distribution function\n", "$$L_\\nu (\\nu)=\\frac{8\\pi r_\\mathrm{boundary\\_inner}^2 h\\nu^3}{c^2}\\frac{1}{\\exp\\left(\\frac{h\\nu}{k_BT_\\mathrm{inner}}\\right)-1}$$\n", "where $L_\\nu$ is the luminosity density (see [Basic Spectrum Generation](../spectrum/basic.ipynb)) with respect to frequency, $\\nu$ is frequency, $h$ is Planck's constant, $c$ is the speed of light, and $k_B$ is Boltzmann's constant:" ] }, { "cell_type": "code", "execution_count": null, "id": "916a5e22", "metadata": { "scrolled": false }, "outputs": [], "source": [ "h = const.h.cgs\n", "c2 = const.c.cgs**2\n", "kB = const.k_B.cgs\n", "\n", "def planck_function(nu):\n", " return 8 * np.pi**2 * r_boundary_inner**2 * h * nu**3 / (c2 * (np.exp(h * nu / (kB * temperature_inner)) - 1))" ] }, { "cell_type": "markdown", "id": "78230177", "metadata": {}, "source": [ "We plot the Planck distribution and a histogram of the generated packet distribution:" ] }, { "cell_type": "code", "execution_count": null, "id": "913fcdbb", "metadata": {}, "outputs": [], "source": [ "# We set important quantites for making our histogram\n", "bins = 200\n", "nus_planck = np.linspace(min(nus), max(nus), bins)\n", "bin_width = nus_planck[1] - nus_planck[0]\n", "\n", "# In the histogram plot below, the weights argument is used \n", "# to make sure our plotted spectrum has the correct y-axis scale\n", "plt.hist(nus.value,\n", " bins=bins,\n", " weights=lumin_per_packet/bin_width)\n", "\n", "# We plot the planck function for comparison\n", "plt.plot(nus_planck, planck_function(nus_planck))\n", "\n", "plt.xlabel('Frequency (Hz)')\n", "plt.ylabel('Luminosity density w.r.t. frequency (erg/s/Hz)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ad4f0f0e", "metadata": {}, "source": [ "We finally plot the generated $\\mu$ density distribution, followed by the generated $\\theta=\\arccos (\\mu)$ density distribution, compared with the respective curves $\\rho = 2\\mu$ and $\\rho = \\sin(2\\theta)$:" ] }, { "cell_type": "code", "execution_count": null, "id": "ae4c97de", "metadata": {}, "outputs": [], "source": [ "x = np.linspace(0, 1, 1000)\n", "\n", "plt.hist(mus, bins=bins, density=True)\n", "plt.plot(x, 2*x)\n", "plt.xlabel('Propagation direction')\n", "plt.ylabel('Probability density')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "41daa433", "metadata": {}, "outputs": [], "source": [ "thetas = np.linspace(0, np.pi/2, 1000)\n", "\n", "plt.hist(np.arccos(mus), bins=bins, density=True)\n", "plt.plot(thetas, np.sin(2*thetas))\n", "plt.xlabel('Angle with normal (rad)')\n", "plt.ylabel('Probability density')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2f661592", "metadata": {}, "source": [ "## Custom Packet Source\n", "\n", "TARDIS allows for the user to input a custom function that generates energy packets instead of the basic blackbody source described here. See [Running TARDIS with a Custom Packet Source](../../io/optional/custom_source.ipynb) for more information." ] } ], "metadata": { "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.7.12" } }, "nbformat": 4, "nbformat_minor": 5 }