{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ELAIS-N1 Luminosity Function\n",
    "\n",
    "Use the depth maps to get a histogram of areas with a given depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n",
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "import glob\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table, join\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "from pymoc import MOC\n",
    "import healpy as hp\n",
    "#import pandas as pd #Astropy has group_by function so apandas isn't required.\n",
    "import seaborn as sns\n",
    "\n",
    "import warnings\n",
    "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from herschelhelp_internal.utils import inMoc, coords_to_hpidx, flux_to_mag, mag_to_flux\n",
    "from herschelhelp_internal.masterlist import find_last_ml_suffix, nb_ccplots\n",
    "\n",
    "from astropy.io.votable import parse_single_table\n",
    "\n",
    "from pcigale.sed import SED\n",
    "from pcigale.sed_modules import get_module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "os.environ['GAMA_DATA'] = 'We are not using GAMA data'\n",
    "#from luminosity_function.gal_sample import CosmoLookup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "FIELD = 'ELAIS-N1'\n",
    "FILTERS_DIR = \"/opt/herschelhelp_python/database_builder/filters/\"\n",
    "DMU_DIR = '/mnt/hedam/dmu_products/'\n",
    "#DMU_DIR = '/Users/rs548/GitHub/dmu_products/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths = Table.read(\"{}dmu1/dmu1_ml_ELAIS-N1/data/depths_elais-n1_20180216.fits\".format(DMU_DIR))\n",
    "final_cat = Table.read(\"{}dmu32/dmu32_ELAIS-N1/data/ELAIS-N1_20171016.fits\".format(DMU_DIR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Histogram of areas\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "depths = depths[\"hp_idx_O_13\", \n",
    "                \"hp_idx_O_10\", \n",
    "                \"ferr_ap_irac_i1_mean\", \n",
    "                \"f_ap_irac_i1_p90\", \n",
    "                \"ferr_irac_i1_mean\", \n",
    "                \"f_irac_i1_p90\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_hist_plot = sns.distplot(depths[\"ferr_ap_irac_i1_mean\"][~np.isnan(depths[\"ferr_ap_irac_i1_mean\"])])\n",
    "depth_hist_plot.set_xlim(0,5.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0.,2.,1000)\n",
    "depth_histogram = np.histogram(depths[\"ferr_ap_irac_i1_mean\"][~np.isnan(depths[\"ferr_ap_irac_i1_mean\"])], bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(depths[\"ferr_ap_irac_i1_mean\"][~np.isnan(depths[\"ferr_ap_irac_i1_mean\"])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_histogram[1][:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ax = sns.regplot(depth_histogram[1][:-1], depth_histogram[0])\n",
    "ax.set(xlabel='ferr_ap_irac_i1_mean (uJ)', ylabel='N')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Vmax = Table.read(DMU_DIR + 'dmu28/dmu28_ELAIS-N1/data/zphot/HELP_final_results.fits')#['id']\n",
    "Vmax['id'].name = 'help_id'\n",
    "Vmax = Vmax['help_id','UVoptIR_bayes.dust.luminosity', 'UVoptIR_bayes.dust.luminosity_err']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vmax[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. First calaculate the irac_i1 flux as a function of redshift for each object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#linearly spaced z - should this be logspace?\n",
    "redshifts = np.linspace(0, 4, 100)\n",
    "Vmax.add_column(Column(data=np.full((len(Vmax), len(redshifts)), \n",
    "                                    np.full(len(redshifts), np.nan)\n",
    "                                   ) , \n",
    "                       name='f_z_relation'\n",
    "                      )\n",
    "               )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vmax[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redshifts = np.linspace(0, 4, 100)\n",
    "n_absent = 0\n",
    "n_processed = 0\n",
    "\n",
    "for gal in Vmax['help_id']:\n",
    "\n",
    "    try:\n",
    "        orig_spec = Table.read(\"{}{}{}_best_model.fits\".format(DMU_DIR, \n",
    "                                              'dmu28/dmu28_ELAIS-N1/data/zphot/best_model_fits/',\n",
    "                                              gal\n",
    "                                             ))\n",
    "        \n",
    "    except FileNotFoundError:\n",
    "        n_absent += 1\n",
    "        # print('fail')\n",
    "        continue\n",
    "        \n",
    "    #print('{} no fail'.format(gal))\n",
    "    s = SED()\n",
    "    # This is wrong because the best SED we get from CIGALE is redshifted (written by Yannick)\n",
    "    s.add_contribution(\"HELP_SED\", orig_spec['wavelength'], orig_spec['L_lambda_total'])\n",
    "    \n",
    "    fluxes = []\n",
    "    for r in redshifts:\n",
    "        sed = s.copy()\n",
    "        mod = get_module(\"redshifting\", redshift=r)\n",
    "        mod.process(sed)\n",
    "        fluxes.append(sed.compute_fnu('IRAC1'))\n",
    "        \n",
    "    Vmax['f_z_relation'][Vmax['help_id'] == gal] = fluxes\n",
    "        \n",
    "    #print(\"{}:{}\".format(gal,fluxes[0]))\n",
    "    n_processed +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('{} processed and {} missing'.format(n_processed, n_absent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Then calculate the zmax for each depth bin for each object\n",
    "\n",
    "Given the array of fluxes as a function of redhsift we interpolate at the depths based on taking the redshift that the object would have a flux equal to 5$\\sigma$ for that mean error bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example object left over as the last object from the previous loop\n",
    "np.interp(5 * bins *1.e-3, np.flip(fluxes,0), np.flip(redshifts,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Vmax.add_column(Column(data=np.full((len(Vmax), len(bins)), \n",
    "                                    np.full(len(bins), np.nan)\n",
    "                                   ) , \n",
    "                       name='zmax_histogram'\n",
    "                      )\n",
    "               )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for gal in Vmax['help_id']:\n",
    "    fluxes = np.array(Vmax[Vmax['help_id']==gal]['f_z_relation'])[0]\n",
    "    z_max_f_relation = np.interp(5 * bins *1.e-3, np.flip(fluxes,0), np.flip(redshifts,0))\n",
    "    Vmax['zmax_histogram'][Vmax['help_id'] == gal] = z_max_f_relation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check the test object above is correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array(Vmax[Vmax['help_id']=='HELP_J155700.909+550039.328']['zmax_histogram'])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV. Finally calculate the Vmax on EN1 for each object\n",
    "For every object calculate the Vmax for each depth cell and multiply by the number of cells at that depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Vmax.add_column(Column(data=np.full((len(Vmax), len(bins)), \n",
    "                                    np.full(len(bins), np.nan)\n",
    "                                   ) , \n",
    "                       name='vmax_histogram'\n",
    "                      )\n",
    "               )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cosmo = FlatLambdaCDM(H0=100. , Om0 = (1-0.7))\n",
    "\n",
    "#TODO: Make a lookup table for speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for gal in Vmax['help_id']:\n",
    "    \n",
    "    z_max_f_relation = np.array(Vmax[Vmax['help_id']==gal]['zmax_histogram'])[0]\n",
    "    v_max_f_relation = cosmo.comoving_volume(z_max_f_relation)\n",
    "    Vmax['vmax_histogram'][Vmax['help_id'] == gal] = v_max_f_relation.value\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vmax.write(\"data/vmax_ELAIS-N1.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#def gal_sample(z0,z1,Ldust,Vmax):\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}