{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7f4a6a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b02f8348",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units\n",
    "import numpy as np"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "37fbb2ca",
   "metadata": {},
   "source": [
    "# Flux Definition"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "fad6f642",
   "metadata": {},
   "source": [
    "SkyLLH provides a sophisticated class collection to describe a differential particle flux function."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4ee1511d",
   "metadata": {},
   "source": [
    "## General Flux Function"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e4765264",
   "metadata": {},
   "source": [
    "The most general differential flux function is:\n",
    "\n",
    "$\\frac{\\mathrm{d}^4\\Phi(\\alpha,\\delta,E,t | \\vec{p}_\\mathrm{s})}{\\mathrm{d}A \\mathrm{d}\\Omega \\mathrm{d}E \\mathrm{d}t}$,\n",
    "\n",
    "which is a function of celestial coordinates right-ascention, $\\alpha$, and declination, $\\delta$, energy $E$, and time $t$, given source parameters $\\vec{p}_\\mathrm{s}$, e.g. source location and spectral index $\\gamma$ for a power-law energy profile."
   ]
  },
  {
   "cell_type": "raw",
   "id": "1e06fcf2",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The abstract base class for all flux function models is \n",
    ":py:class:`~skyllh.core.flux_model.FluxModel`, which is derived from \n",
    ":py:class:`~skyllh.core.math.MathFunction` and \n",
    ":py:class:`~skyllh.core.model.Model`."
   ]
  },
  {
   "cell_type": "raw",
   "id": "007a68ec",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The :py:class:`~skyllh.core.flux_model.FluxModel` has the abstract \n",
    ":py:meth:`~skyllh.core.flux_model.FluxModel.__call__` method defined, which will\n",
    "evaluate the flux function for values of :math:`\\alpha`, :math:`\\delta`, \n",
    ":math:`E`, and :math:`t` given in specific units."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e94984fc",
   "metadata": {},
   "source": [
    "## Units for flux models"
   ]
  },
  {
   "cell_type": "raw",
   "id": "f5a5c595",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The :py:class:`~skyllh.core.flux_model.FluxModel` class defines the units used \n",
    "for length, angle, energy, and time. The default units are configured through \n",
    "a :py:class:`skyllh.core.config.Config` instance, in particular through the\n",
    "``cfg['units']['defaults']['fluxes']`` dictionary of the \n",
    ":py:class:`~skyllh.core.config.Config` instance. Units must be derived classes \n",
    "from the :py:class:`astropy.units.UnitBase` class."
   ]
  },
  {
   "cell_type": "raw",
   "id": "aee5873d",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The :py:class:`~skyllh.core.flux_model.FluxModel` class has the properties \n",
    ":py:attr:`~skyllh.core.flux_model.FluxModel.unit_str` and \n",
    ":py:attr:`~skyllh.core.flux_model.FluxModel.unit_latex_str` for a representation\n",
    "of the units as a ``str`` object in plain text and latex code."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8f4e7358",
   "metadata": {},
   "source": [
    "## Factorized Flux Function"
   ]
  },
  {
   "cell_type": "raw",
   "id": "98ebd803",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Usually the flux function can be split into a spatial, energy, and time profile \n",
    "with an overall normalization constant in differential flux unit, i.e. \n",
    ":math:`\\mathrm{area}^{-1} \\mathrm{solid-angle}^{-1} \\mathrm{energy}^{-1} \\mathrm{time}^{-1}`. \n",
    "Hence, SkyLLH provides the class \n",
    ":py:class:`~skyllh.core.flux_model.FactorizedFluxModel`, which describes a \n",
    "differential flux function as the product of individual flux profiles:\n",
    "\n",
    ".. math::\n",
    "\n",
    "    \\Phi(\\alpha,\\delta,E,t | \\vec{p}_\\mathrm{s}) =\n",
    "            \\Phi_0\n",
    "            \\Psi(\\alpha,\\delta|\\vec{p}_\\mathrm{s})\n",
    "            \\epsilon(E|\\vec{p}_\\mathrm{s})\n",
    "            T(t|\\vec{p}_\\mathrm{s})\n",
    "\n",
    "The abstract base class for any flux profile is \n",
    ":py:class:`~skyllh.core.flux_model.FluxProfile`, which is derived from \n",
    ":py:class:`~skyllh.core.math.MathFunction`.\n",
    "\n",
    "The abstract base class for a spatial, energy, and time flux profile is \n",
    ":py:class:`~skyllh.core.flux_model.SpatialFluxProfile`, \n",
    ":py:class:`~skyllh.core.flux_model.EnergyFluxProfile`, and \n",
    ":py:class:`~skyllh.core.flux_model.TimeFluxProfile`, respectively, and are \n",
    "derived from :py:class:`~skyllh.core.flux_model.FluxProfile`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3242d1d9",
   "metadata": {},
   "source": [
    "### Steady Point-Like Flux"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9f489769",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "A very common source hypothesis is a steadily emitting point-like source. Hence, \n",
    "SkyLLH provides the class \n",
    ":py:class:`~skyllh.core.flux_model.SteadyPointlikeFFM`. It takes a flux \n",
    "normalization :math:`\\Phi_0`, and an energy profile as constructor arguments. \n",
    "As spatial profile it uses the \n",
    ":py:class:`~skyllh.core.flux_model.PointSpatialFluxProfile` class.\n",
    "\n",
    "As an example we create a steady point-like factorized flux model with a \n",
    "power-law energy flux profile.\n",
    "\n",
    "But first we need to create a configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3ffd47bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.config import Config\n",
    "cfg=Config()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "raw",
   "id": "1c80f05e",
   "metadata": {},
   "source": [
    "Now we see what the default units are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3a336a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.flux_model import FluxModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9b799b48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'angle': Unit(\"rad\"), 'energy': Unit(\"GeV\"), 'length': Unit(\"cm\"), 'time': Unit(\"s\")}\n"
     ]
    }
   ],
   "source": [
    "print(FluxModel.get_default_units(cfg=cfg))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "raw",
   "id": "4cc57768",
   "metadata": {},
   "source": [
    "Now we need to create the energy flux profile. As reference energy and spectral\n",
    "index we choose $E_0=10^3$~GeV and $\\gamma=2$, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "74d36971",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.flux_model import PowerLawEnergyFluxProfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2c104c38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(E / (1000 GeV))^-2\n"
     ]
    }
   ],
   "source": [
    "energy_profile = PowerLawEnergyFluxProfile(\n",
    "    E0=1e3, \n",
    "    gamma=2, \n",
    "    energy_unit=units.GeV,\n",
    "    cfg=cfg)\n",
    "print(energy_profile)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "ac4b6bb9",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The next step is to create the \n",
    ":py:class:`~skyllh.core.flux_model.SteadyPointlikeFFM` class instance. As \n",
    "normalization constant we choose \n",
    ":math:`\\Phi_0 = 10^{-8} \\text{GeV}^{-1}\\text{cm}^{-2}\\text{sr}^{-1}\\text{s}^{-1}`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "13335852",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skyllh.core.flux_model import SteadyPointlikeFFM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e06b2222",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000e-08 * (E / (1000 GeV))^-2 * 1 (GeV cm^2 s)^-1\n"
     ]
    }
   ],
   "source": [
    "fluxmodel = SteadyPointlikeFFM(\n",
    "    Phi0=1e-8,\n",
    "    energy_profile=energy_profile,\n",
    "    angle_unit=units.radian,\n",
    "    time_unit=units.s,\n",
    "    length_unit=units.cm,\n",
    "    cfg=cfg\n",
    ")\n",
    "print(fluxmodel)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9a282441",
   "metadata": {},
   "source": [
    "## Evaluating the flux model function\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bdbb055a",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The flux model function can be evaluated by calling its \n",
    ":py:meth:`~skyllh.core.flux_model.FluxModel.__call__` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a22a30bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flux.shape = (1, 1, 1)\n",
      "flux = [[[1.11111111e-09]]]\n"
     ]
    }
   ],
   "source": [
    "flux = fluxmodel(E=3e3)\n",
    "print(f'flux.shape = {flux.shape}')\n",
    "print(f'flux = {flux}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "raw",
   "id": "9371a4cf",
   "metadata": {},
   "source": [
    "It returns a 3-dimensional numpy array, where the first, second, and third \n",
    "dimension represents the spatial, energy, and time axes, respectively. Hence, \n",
    "we can evaluate the flux model for different spatial coordinates, energies, and \n",
    "times by a single call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6b6596b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flux.shape = (1, 2, 3)\n",
      "flux = [[[1.11111111e-09 1.11111111e-09 1.11111111e-09]\n",
      "  [2.50000000e-07 2.50000000e-07 2.50000000e-07]]]\n"
     ]
    }
   ],
   "source": [
    "flux = fluxmodel(E=[3e3, 2e2], t=[1, 2, 3])\n",
    "print(f'flux.shape = {flux.shape}')\n",
    "print(f'flux = {flux}')"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}