:orphan:

Creating New Source Populations
===============================
If you're interested in using *ARES* to evolve radiation backgrounds at high-:math:`z` but it doesn't have the source populations you're interested in, it shouldn't be too hard to add your own. This section provides a brief tutorial on how to do that.

Before diving in, however, you may want to consider whether the current machinery has all the functionality you need. For example, though the main set of *ARES* source populations were designed to model galaxies, and so their parameters refer to galaxy properties, fundamentally these populations are all connected to dark matter halo populations. So, if that's true of your source-population-of-interest as well, you may be able to trick *ARES* into doing what you want.

The main source populations that you might manipulate are:

* ``GalaxyAggregate``
    Assumes star formation rate density (SFRD) proportional to rate at which matter collapses into halos. Radiation backgrounds proportional to SFRD. Key parameters include ``pop_fstar``, ``pop_rad_yield``, ``pop_fesc``, which are all effectively population-averaged quantities, hence the "aggregate" designation.

* ``GalaxyCohort``
    Assumes 1:1 galaxy:halo mapping, and assumes star formation rate (SFR) in galaxies is proportional to the mass accretion rate onto their parent halo. Luminosity of galaxies is proportional to their SFR, so radiation backgrounds are generated by performing an integral over the luminosity function of galaxies. Key parameters include ``pop_fstar`` (now a parameterizable function; see :doc:`uth_pq`), ``pop_rad_yield``, etc. Unlike ``GalaxyAggregate``, galaxy properties now vary as a function of mass and redshift. But, because all galaxies formed at the same time with the same mass are identical, they still evolve as cohorts, i.e., galaxies don't break mass-rank ordering.

If neither of these will accommodate your goals, continue onward! Perhaps someday your new source population will be listed in the list above.

Step One: A Simple Source Population
------------------------------------
Let's create a new source population with the following properties:

- Objects form at a constant rate (per co-moving volume element).
- Objects have a flat spectrum.

This isn't a terribly interesting source population, but it's one whose behavior we can easily test.

To get started, make a new file, let's call it ``DummyPopulation.py``, which we'll eventually add to the aptly-named ``populations`` sub-module in *ARES*, i.e., in ``$ARES/ares/populations``. For now, you can put this file anywhere.


::

    # A few important imports
    import numpy as np
    from ares.util import ParameterFile
    from ares.populations.Halo import HaloPopulation

    class DummyPopulation(HaloPopulation):
        def __init__(self, **kwargs):
            HaloPopulation.__init__(self, **kwargs)

A few things to notice off the bat. First, we inherited ``HaloPopulation`` (which itself inherits a generic ``Population`` object that holds all sorts of meta-data), and called its constructor. This creates a ``ParameterFile`` instance and stores it in the attribute ``pf``, which is a dictionary that contains all the parameters relevant to this source population. If you haven't supplied any parameters, this dictionary will be filled with a bunch of defaults.

.. note :: Once you move ``DummyPopulation.py`` into the *ARES* source code directory (``$ARES/ares/populations``), you'll need to convert to *relative* imports for *ARES*. For example, ``from ares.util import ParameterFile`` should instead be ``from ..util import ParameterFile``.

Inheritance of the ``HaloPopulation`` gives us immediate access to the halo mass function through the ``halos`` attribute.

To do some testing, copy-paste the above snippet into a Python terminal and then initialize an instance of your class:

::

    pop = DummyPopulation()

    print(pop.halos.tab_z.shape, pop.halos.tab_z.min(), pop.halos.tab_z.max())
    print(pop.halos.tab_M.shape, pop.halos.tab_M.min(), pop.halos.tab_M.max())
    print(pop.halos.tab_dndm.shape)

So, the halo mass function is stored as a lookup table of shape (1145, 1200), with dimensions corresponding to (redshift, halo mass).

At this stage, your source population isn't terribly useful -- you could have just initialized a ``HaloPopulation`` and been done! However, the thing that sets apart source populations from halo populations is the radiation they emit. In order to evolve the cosmic radiation backgrounds generated by sources, the solvers in *ARES* simply loop over all population objects and access their ``Emissivity`` method, which means our ``DummyPopulation`` object needs such a method if it's to be any use to us.

.. note:: Generally, this emissivity is tied to the SFRD, which (by convention) we simply call ``SFRD``, though *ARES* will only go looking for the SFRD if requested by the user.

There are a few rules for the Emissivity function.

- It must take redshift, ``z``, as its first positional argument.
- It must have three additional keyword arguments, ``E``, ``Emin``, and ``Emax``, all set to ``None`` by default.
- It should return the volume-averaged emissivity of the source population in units of :math:`\mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-3}` if ``E==None``. If ``E!=None``, then the user has requested a *specific* emissivity, and so the units should be :math:`\mathrm{erg} \ \mathrm{s}^{-1} \ \mathrm{cm}^{-3} \ \mathrm{eV}^{-1}`. In other words, in the latter case if you integrate the resulting emissivity over an array of photon energies (in eV) from ``pop_Emin`` to ``pop_Emax``, you should recover the luminosity density of the population. In both cases, the :math:`\mathrm{cm}^3` volume should be a *co-moving* volume!

So, here's a start:

::

    def Emissivity(self, z, E=None, Emin=None, Emax=None):
        """
        An informative docstring.
        """

        # Check to see if these sources are on, i.e., is the supplied redshift
        # in the range (pop_zdead, pop_zform).
        on = self.on(z)
        if not np.any(on):
            return z * on

        # Compute the emissivity
        if E is not None:
            # Compute a specific emissivity and return
            raise NotImplemented('do something!')
        else:
            # Compute a band-integrated emissivity and return
            raise NotImplemented('do something else!')


No matter what, the ``Emissivity`` method needs to know about the spectral energy distribution (SED) of the source population. The SED is defined by the user via the ``pop_sed`` parameter, which can be a variety of things (see :doc:`uth_pop_radiation`). This information is parsed by the ``Population`` class and stored in an attribute ``src``, which we have access to already because we inherited ``Population``!

The attribute ``src`` is an instance of yet another class that represents the SED of an object (these are defined in the ``sources`` sub-module of *ARES*). It's most important method is ``Spectrum``, which is a function of photon energy (in eV), and returns a normalized version of the SED. See :doc:`params_sources`` for a more detailed description of how this is done.

So, in completing your ``Emissivity`` function, you'll need to use the ``src.Spectrum`` function to properly account for the SED of your sources.


Step Two: Integrating with *ARES*
---------------------------------
In the previous section we created a stand-alone (but useless) source population. In order for it to work within *ARES*, we need to make a few changes in other areas of the code.

First, we must decide how the user is to indicate that this population is of interest, which means setting the ``pop_sfr_model`` parameter. Let's set things up so that passing in ``pop_sfr_model='dummy'`` will trigger the creation of our ``DummyPopulation`` source. To make sure *ARES* knows about this, we need to navigate to ``ares.populations.GalaxyPopulation.py``.

.. note :: Again, since the sources are generally assumed to be galaxies, whose luminosity is mostly from star formation, the main parameter is ``pop_sfr_model`` and the class that initializes a "generic" source population is ``GalaxyPopulation``. Apologies if your model does not invoke star formation or galaxies!

Near the bottom of ``ares.populations.GalaxyPopulation.py``, there is a series of if/else statements that are checking the value of ``pop_sfr_model``, and initializing the appropriate class depending on its value. For example, if ``pop_sfr_model=='sfe-func'`` we initialize a ``GalaxyCohort, if ``pop_sfr_model=='fcoll'`` we initialize a ``GalaxyAggregate``, and so on. Within this if/else block, you need only add

::

    elif model in ['dummy']:
        return DummyPopulation(**kwargs)

You'll of course need to be sure to add an ``import`` statement for ``DummyPopulation`` at the top of the file.

Finally, if you haven't already, move ``DummyPopulation.py`` into the ``ares.populations`` module, and make sure to convert to relative imports (as noted in previous section).


Step Three: Testing the new population
--------------------------------------
First, let's make sure we can initialize an instance of the new source population through *ARES*:

::

    import ares

    pop = ares.populations.GalaxyPopulation(pop_sfr_model='dummy')

and verify that its routines behave as expected.

Now, to verify that the population works within an *ARES* simulation, let's compare the results of two calculations: one standard calculation, and the same calculation with this new source population.

Under construction!