{ "cells": [ { "cell_type": "markdown", "id": "1a955596", "metadata": {}, "source": [ "# Basic Spectrum Generation\n", "\n", "The first and most simple way in which TARDIS calculates spectra is by calculating it directly from the Monte Carlo packets after the final [Monte Carlo Iteration](../montecarlo/index.rst). This simply requires knowledge of each packet's energy and frequency in the lab frame (see [Reference Frames](../montecarlo/propagation.rst#reference-frames)) at the end of the iteration. The only other quantity needed is the time duration of the simulation $\\Delta t$, which is calculated based off of the luminosity of the supernova's photosphere (see [Energy Packet Initialization](../montecarlo/initialization.ipynb)).\n", "\n", "
\n", "\n", "Note\n", "\n", "The only packets which are used for this calculation are the packets which escape the outer boundary of the\n", "computational domain -- those reabsorbed into the photosphere are not included (see [Packet Propagation](../montecarlo/propagation.rst)).\n", "\n", "
\n", "\n", "The spectrum calculation is very straightforward. A packet of energy $E_\\mathrm{packet}$ contributes a\n", "luminosity\n", "\n", "$$L_\\mathrm{packet} = \\frac{E_\\mathrm{packet}}{\\Delta t}$$\n", "\n", "to the spectrum at its frequency.\n", "\n", "In the code below, we will see an issue with merely relying on luminosity to give us a spectrum, which will allow us to develop the concept of luminosity *density*, and then correctly plot the TARDIS spectrum." ] }, { "cell_type": "markdown", "id": "f5196349", "metadata": {}, "source": [ "We start by importing the necessary packages, loading a configuration, and setting up a simulation object (see [Setting up the Simulation](../setup/index.rst)):" ] }, { "cell_type": "code", "execution_count": null, "id": "b4e90662", "metadata": {}, "outputs": [], "source": [ "from tardis.io.config_reader import Configuration\n", "from tardis.simulation import Simulation\n", "from tardis.montecarlo import TARDISSpectrum\n", "from tardis.io.atom_data.util import download_atom_data\n", "from astropy import units as u\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "id": "5a3adb2e", "metadata": {}, "outputs": [], "source": [ "# We download the atomic data needed to run the simulation\n", "download_atom_data('kurucz_cd23_chianti_H_He')" ] }, { "cell_type": "code", "execution_count": null, "id": "c96cfbd2", "metadata": {}, "outputs": [], "source": [ "tardis_config = Configuration.from_yaml('tardis_example.yml')\n", "\n", "sim = Simulation.from_config(tardis_config)" ] }, { "cell_type": "markdown", "id": "f169c41b", "metadata": {}, "source": [ "We now select a number of packets to run through the Monte Carlo simulation, and then run one Monte Carlo iteration (see [Monte Carlo Iteration](../montecarlo/index.rst)):" ] }, { "cell_type": "code", "execution_count": null, "id": "82279c43", "metadata": {}, "outputs": [], "source": [ "N_packets = 5000\n", "\n", "# Using the commented out code below, we can also get the number of packets\n", "# from the configuration -- try it out:\n", "#N_packets = tardis_config.no_of_packets\n", "\n", "sim.iterate(N_packets)" ] }, { "cell_type": "markdown", "id": "de4cd415", "metadata": {}, "source": [ "We call the arrays of each packet's frequency and how much energy the packet has:" ] }, { "cell_type": "code", "execution_count": null, "id": "c60863af", "metadata": {}, "outputs": [], "source": [ "nus = sim.runner.output_nu\n", "nus" ] }, { "cell_type": "code", "execution_count": null, "id": "ac231055", "metadata": {}, "outputs": [], "source": [ "energies = sim.runner.output_energy\n", "energies" ] }, { "cell_type": "markdown", "id": "6fc8a089", "metadata": {}, "source": [ "Notice that some energies are negative. This means that the packet ended up being reabsorbed into the photosphere (we will separate out these packets later). Also note that the number of elements of our arrays of frequencies and energies is the same as the number of packets that we ran (as it should be):" ] }, { "cell_type": "code", "execution_count": null, "id": "1b628b04", "metadata": {}, "outputs": [], "source": [ "len(nus), len(energies)" ] }, { "cell_type": "markdown", "id": "70c979d2", "metadata": {}, "source": [ "
\n", "\n", "Note\n", "\n", "The energies used for the spectrum are in the lab frame (see [Reference Frames](../montecarlo/propagation.rst#reference-frames)). Recall that while packets are all initialized with the same energy in the lab frame, throughout the simulation energy is conserved in the co-moving frame, allowing the lab frame energies of the packets to change though interactions (see [Performing an Interaction](../montecarlo/propagation.rst#performing-an-interaction)), creating the varied lab-frame energies that we see at the end of the simulation.\n", " \n", "
" ] }, { "cell_type": "markdown", "id": "54466d70", "metadata": {}, "source": [ "TARDIS will then calculate a list of the packet luminosities by dividing each element in the array of energies by the time of the simulation $\\Delta t$:" ] }, { "cell_type": "code", "execution_count": null, "id": "ba7a7d93", "metadata": {}, "outputs": [], "source": [ "luminosities = energies / sim.runner.time_of_simulation\n", "luminosities" ] }, { "cell_type": "markdown", "id": "e10d9567", "metadata": {}, "source": [ "Now, as mentioned before, we only want to include the packets that make it through to the outer boundary. To do this, TARDIS creates an array of booleans (either True or False) called a mask that tells us if the packet should be counted in the spectrum. We then can use that mask to get an array of the frequencies and energies of only the packets which we are interested in:" ] }, { "cell_type": "code", "execution_count": null, "id": "b69e8c92", "metadata": {}, "outputs": [], "source": [ "emitted_mask = sim.runner.emitted_packet_mask\n", "emitted_mask" ] }, { "cell_type": "code", "execution_count": null, "id": "c856a85e", "metadata": {}, "outputs": [], "source": [ "emitted_nus = nus[emitted_mask]\n", "emitted_nus" ] }, { "cell_type": "code", "execution_count": null, "id": "38141ad8", "metadata": {}, "outputs": [], "source": [ "emitted_luminosities = luminosities[emitted_mask]\n", "emitted_luminosities" ] }, { "cell_type": "markdown", "id": "f7cc4807", "metadata": {}, "source": [ "The length of these lists is the number of packets that made it out of the supernova:" ] }, { "cell_type": "code", "execution_count": null, "id": "b4c6f180", "metadata": {}, "outputs": [], "source": [ "len(emitted_nus), len(emitted_luminosities)" ] }, { "cell_type": "markdown", "id": "2a29c4c9", "metadata": {}, "source": [ "Now, let's plot frequency versus luminosity. We will see a very strange graph, which will lead us into developing a new strategy for plotting the spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "782099c5", "metadata": {}, "outputs": [], "source": [ "plt.scatter(emitted_nus, emitted_luminosities)\n", "plt.xlabel('Frequency (Hz)')\n", "plt.ylabel('Luminosity (erg/s)');" ] }, { "cell_type": "markdown", "id": "ee595852", "metadata": {}, "source": [ "This is not the type of plot that we are looking for. We also cannot solve this problem by adding up the luminosities of all packets with the same frequency -- in fact, since frequency is a continuum, it would be extremely unlikely that any two packets have the same *exact* frequency. To solve this problem, we will need the concept of luminosity density.\n", "\n", "## Luminosity Density\n", "\n", "What we will have to do instead is plot a histogram where we bin up different frequencies that are close to each other into a certain number of bins. We then add up the luminosities of each packet in the bin and divide by the \"width\" of the bins. For example, if we are plotting between 0 Hz and 50 Hz in 5 bins, we would add up the luminosities of the packets between 0 Hz and 10 Hz, between 10 Hz and 20 Hz, between 20 Hz and 30 Hz, etc., and then divide each value by 10 Hz, which is the width of each bin (note that the bin widths need not be uniform). This will give us the luminosity density with respect to frequency, denoted by $L_\\nu$, and measured in ergs per second per Hertz. This can be interpreted as the luminosity per unit Hertz, i.e. how much luminosity will be in an interval with a width of 1 Hz.\n", "\n", "The division step here is crucial, as otherwise the values on the y-axis will, for example, approximately double if we double the widths of our bins (as we would be adding about double the luminosity contributions into that bin). We clearly want the values on the y-axis to be independent of the number of bins we break the luminosity into, thus making luminosity density the best way to plot a spectrum.\n", "\n", "Note that we can also have luminosity density with respect to wavelength, $L_\\lambda$, which we get by binning the packets by wavelength instead of frequency. Since the width of the bins would now have the dimensions of length, $L_\\lambda$ will have units of ergs per second per Angstrom." ] }, { "cell_type": "markdown", "id": "dbc69e96", "metadata": {}, "source": [ "Now, to generate our spectrum, we select the bounds for our spectrum and the number of bins that we group the packets into. Feel free to change the number of bins to see how it affects the spectrum!" ] }, { "cell_type": "code", "execution_count": null, "id": "b43dcb8d", "metadata": {}, "outputs": [], "source": [ "# The lowest frequency we plot\n", "freq_start = 1.5e14 * u.Hz\n", "\n", "# The highest frequency we plot\n", "freq_stop = 3e15 * u.Hz\n", "\n", "# The number of bins\n", "N = 500" ] }, { "cell_type": "markdown", "id": "2e375d04", "metadata": {}, "source": [ "The above information can also be retrieved from the configuration using the commented-out code below. Note that the configuration has the bounds specified in terms of wavelengths. Therefore we must convert these to frequencies, noting that the frequency is the speed of light divided by the wavelength (Astropy has a built-in way to do this, which we shall use). Additionally, since wavelength and frequency are inversely related, the lower bound for the wavelengths is the upper bound for the frequencies and vice versa." ] }, { "cell_type": "code", "execution_count": null, "id": "e763bbf1", "metadata": {}, "outputs": [], "source": [ "#freq_start = tardis_config.spectrum.stop.to('Hz', u.spectral())\n", "\n", "#freq_stop = tardis_config.spectrum.start.to('Hz', u.spectral())\n", "\n", "#N = tardis_config.spectrum.num" ] }, { "cell_type": "markdown", "id": "30e8090a", "metadata": {}, "source": [ "Next, TARDIS generates the list of frequency bins. The array shown contain the boundaries between successive bins as well as the lower and upper bounds of the spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "226c63c8", "metadata": {}, "outputs": [], "source": [ "spectrum_frequency = np.linspace(freq_start, freq_stop, N+1)\n", "spectrum_frequency" ] }, { "cell_type": "markdown", "id": "12663ddd", "metadata": {}, "source": [ "Then, TARDIS creates a histogram where we add up the luminosity in each bin:" ] }, { "cell_type": "code", "execution_count": null, "id": "2a86d791", "metadata": {}, "outputs": [], "source": [ "emitted_luminosity_hist = u.Quantity(np.histogram(emitted_nus,\n", " weights=emitted_luminosities,\n", " bins=spectrum_frequency,\n", " )[0], \"erg / s\",)\n", "emitted_luminosity_hist" ] }, { "cell_type": "markdown", "id": "59c0fbfd", "metadata": {}, "source": [ "Finally, we input this information into the TARDISSpectrum class, which will generate the luminosity density with respect to both wavelength and frequency and allow us to plot both of these." ] }, { "cell_type": "code", "execution_count": null, "id": "392be1e4", "metadata": {}, "outputs": [], "source": [ "spectrum = TARDISSpectrum(spectrum_frequency, emitted_luminosity_hist)" ] }, { "cell_type": "code", "execution_count": null, "id": "10e3fe34", "metadata": {}, "outputs": [], "source": [ "spectrum.plot(mode='frequency')" ] }, { "cell_type": "code", "execution_count": null, "id": "599ff77a", "metadata": {}, "outputs": [], "source": [ "spectrum.plot(mode='wavelength')" ] }, { "cell_type": "markdown", "id": "a53e2b79", "metadata": {}, "source": [ "
\n", "\n", "Note\n", "\n", "Most of this process is done internally by TARDIS. Given a simulation object `sim` that has been run, calling `sim.runner.spectrum` will give you a `TARDISSpectrum` object that can then be plotted using the `.plot()` method. See, for example, our [Quickstart Guide](../../quickstart.ipynb). This notebook just demonstrates how TARDIS generates this spectrum when `sim.runner.spectrum` is called.\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "49b60fed", "metadata": {}, "source": [ "You may notice that the bins are not uniformly spaced with respect to wavelength. This is because we use the same bins for both wavelength and frequency, and thus the bins which are uniformly spaced with respect to frequency are not uniformly spaced with respect to wavelength. This is okay though, since luminosity density allows us to alter the bins without significantly altering the graph!\n", "\n", "Another thing you may notice in these graphs is the lack of a smooth curve. This is due to **noise**: the effects of the random nature of Monte Carlo simulations. This makes it very difficult to get a precise spectrum without drastically increasing the number of Monte Carlo packets. To solve this problem, TARDIS uses [virtual packets](virtualpackets.rst) and [the formal integral method](formal_integral.rst) to generate a spectrum with less noise." ] } ], "metadata": { "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 }