{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating galaxy density maps using the method described in Duivenvoorden et al. 2016 (this method is based on Darvish et al. 2015)\n",
    "\n",
    "\n",
    "![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=100&v=4>)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import os\n",
    "from astropy.io import fits\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "import pprint\n",
    "from scipy.special import erf\n",
    "from scipy.integrate import simps, trapz\n",
    "from scipy.stats.mstats import gmean\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline  \n",
    "\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3)\n",
    "\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "cdict = {'red': ((0., 0.8, 0.8),\n",
    "                 (0.08, 0.8, 0.8),\n",
    "                 (0.3, 0.0, 0.0),\n",
    "                 (1, 1, 1)),\n",
    "         'green': ((0., 1.0, 1.0),\n",
    "                 (0.08, 0.8, 0.8),\n",
    "                   (0.3, 0.0, 0.0),\n",
    "                   (0.73, 0.0, 0.0),\n",
    "                   (1, 0, 0)),\n",
    "         'blue': ((0., 0.8, 0.8),\n",
    "                 (0.08, 1, 1),\n",
    "                  (0.45, 0.5, 0.5),\n",
    "                  (0.65, 0.0, 0.0),\n",
    "                  (1, 0, 0))}\n",
    "\n",
    "my_cmap = LinearSegmentedColormap('my_colormap',cdict,120)\n",
    "name_field = 'COSMOS'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Codes uses the HELP masterlist to select objects with photometric redshifts. \n",
    "### The data needs to be homogeneous, to avoid peaks in the density maps caused by deeper data. \n",
    "### A mass cut (or IRAC/K-band) can be used to select galaxies with accurate photometric redshifts. \n",
    "### In this example we assume a Gaussian shape of the redshift PDF, the full redshift PDFs can be found\n",
    "### soon at HeDaM to get a more accurate (but slower) result. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/ipykernel_launcher.py:16: RuntimeWarning: invalid value encountered in less\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "310367\n",
      "0.0632861860603\n"
     ]
    }
   ],
   "source": [
    "loc = ''\n",
    "hdulist = fits.open(loc+'cosmos_allz.fits') # test file whihc contains the phot-z in COSMOS\n",
    "img = hdulist[1].data\n",
    "cols = hdulist[1].columns\n",
    "hdulist.close()\n",
    "z_med = np.array(img['zm_eazy'])\n",
    "z_low = np.array(img['l68_eazy']) # 1s lower limit\n",
    "z_high = np.array(img['u68_eazy']) # 1s uper limit\n",
    "ra = np.array(img['RA_2'])\n",
    "dec = np.array(img['DEC_2'])\n",
    "mag1 = np.array(img['m_vista_ks']) # K band data\n",
    "mag2 = np.array(img['m_wircam_ks'])\n",
    "mag3 = np.array(img['m_ukidss_k'])\n",
    "\n",
    "# make the selection, in this user case a simple K band < 24 selection.\n",
    "use = (z_med >= 0.01) & (z_high <= 3.5) & ( (mag1 < 24.0) | (mag2 < 24) | (mag3 < 24)) \n",
    "\n",
    "err = (z_high[use] - z_low[use])/2. # Gaussian assumption for PDF\n",
    "z_med = z_med[use]\n",
    "ra = ra[use]\n",
    "dec = dec[use]\n",
    "print(np.size(ra))\n",
    "print(np.mean(err/(1+z_med)))\n",
    "\n",
    "\n",
    "length = 3900 # create a grid between z = 0.1 and z = 4\n",
    "rd = 3\n",
    "z_grid = np.linspace(0.1,4.0, length+1) \n",
    "idx = np.arange(length+1)\n",
    "\n",
    "zmin=0\n",
    "zmax=1\n",
    "z_slice_low = []\n",
    "z_slice_high = []\n",
    "z_slice_mean= []\n",
    "# create redshift slics with a width of 2 times the median error\n",
    "while(z_grid[zmax] < 4): \n",
    "    use = (z_med > z_grid[zmin]) & (z_med <= z_grid[zmax])\n",
    "    if np.size(err[use]) < 200:\n",
    "        zmax = zmax + 1\n",
    "    else:            \n",
    "        med = np.median(err[use])\n",
    "        if 2*med < (z_grid[zmax]-z_grid[zmin]):\n",
    "            z_slice_low = np.append(z_slice_low,z_grid[zmin])\n",
    "            z_slice_high = np.append(z_slice_high,z_grid[zmax])\n",
    "            mn = z_grid[zmin]+(z_grid[zmax]-z_grid[zmin])/2\n",
    "            delt =  np.abs(z_grid - round(mn,rd)) == np.min(np.abs(z_grid - round(mn,rd)))\n",
    "            z_slice_mean = np.append(z_slice_mean,z_grid[delt])\n",
    "            zmin = idx[delt][0]\n",
    "        zmax = zmax + 1\n",
    "        \n",
    "# saves redshift slices        \n",
    "DAT = np.array([z_slice_low,z_slice_high]).T\n",
    "header = 'z_low z_high'\n",
    "np.savetxt('z_slice_'+name_field+'.cat', DAT, delimiter=\" \", fmt='%s',  header= header, newline='\\n')        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of mock maps\n",
    "# 40 were used in the original Duivenvoorden et al. 2016 paper\n",
    "# this verion uses 5\n",
    "mock_maps = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-317a108b9783>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     37\u001b[0m         \u001b[0mweights_use\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     38\u001b[0m         \u001b[0mc2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSkyCoord\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mra_use_fs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdeg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdec_use_fs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdeg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 39\u001b[0;31m         \u001b[0msep\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mc1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mseparation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     40\u001b[0m         \u001b[0mdist_i\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msep\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marcmin\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdist\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     41\u001b[0m         \u001b[0mK_weigt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mweights_use\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdist_i\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mweights_use\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/coordinates/sky_coordinate.py\u001b[0m in \u001b[0;36mseparation\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    939\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    940\u001b[0m         \u001b[0;31m# Get the separation as a Quantity, convert to Angle in degrees\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 941\u001b[0;31m         \u001b[0msep\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mangular_separation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlon2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlat2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    942\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mAngle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msep\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdegree\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    943\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/coordinates/angle_utilities.py\u001b[0m in \u001b[0;36mangular_separation\u001b[0;34m(lon1, lat1, lon2, lat2)\u001b[0m\n\u001b[1;32m    699\u001b[0m     \u001b[0mslat1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlat1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    700\u001b[0m     \u001b[0mslat2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlat2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 701\u001b[0;31m     \u001b[0mclat1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlat1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    702\u001b[0m     \u001b[0mclat2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlat2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    703\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/astropy/coordinates/angles.py\u001b[0m in \u001b[0;36m__array_wrap__\u001b[0;34m(self, obj, context)\u001b[0m\n\u001b[1;32m    540\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    541\u001b[0m     \u001b[0;31m# Any calculation should drop to Angle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 542\u001b[0;31m     \u001b[0;32mdef\u001b[0m \u001b[0m__array_wrap__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcontext\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    543\u001b[0m         \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__array_wrap__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcontext\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcontext\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    544\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_no_angle_subclass\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "lengt = 50 #proper kpc, grid for calculating density\n",
    "for g in range(1,np.size(z_slice_low)):  \n",
    "  z_l = z_slice_low[g]\n",
    "  z_h = z_slice_high[g]\n",
    "\n",
    "  weights = (0.5*(1+erf((z_h-z_med)/(err*np.sqrt(2))))-0.5*(1+erf((z_l-z_med)/(err*np.sqrt(2)))))\n",
    "  dist = np.array(cosmo.kpc_proper_per_arcmin(z_l/2+z_h/2))\n",
    "  dist_kpc_dec = (np.max(dec)-np.min(dec))*60*dist # dont use if dec crosses -90,90 \n",
    "  dist_kpc_ra = (np.max(ra)-np.min(ra))*60*dist # dont use if ra crosses 0, 360\n",
    "  x_edg = np.linspace(np.min(ra),np.max(ra), ((dist_kpc_ra)/lengt).astype(int)+1) \n",
    "  y_edg = np.linspace(np.min(dec),np.max(dec), ((dist_kpc_dec)/lengt).astype(int)+1)\n",
    "\n",
    "  use_for_speed = (weights > 0.001)\n",
    "\n",
    "  ra_use_fs = ra[use_for_speed]\n",
    "  dec_use_fs = dec[use_for_speed]\n",
    "\n",
    "  x_p = np.zeros(np.size(x_edg) - 1)\n",
    "  y_p = np.zeros(np.size(y_edg) - 1)\n",
    "  for i in range(np.size(x_edg) - 1):\n",
    "    x_p[i] = x_edg[i]/2 + x_edg[i+1]/2\n",
    "  for i in range(np.size(y_edg) - 1):\n",
    "    y_p[i] = y_edg[i]/2 + y_edg[i+1]/2\n",
    "\n",
    "  for n in range(mock_maps+1):#41\n",
    "    if n > 0:\n",
    "      random_list = np.int_(np.rint((np.size(ra)-1)*np.random.random([np.size(ra_use_fs)])))\n",
    "      ra_use_fs = ra[random_list]\n",
    "      dec_use_fs = dec[random_list]\n",
    "\n",
    "      \n",
    "    dens_p_g = np.zeros(np.size(weights[use_for_speed]))\n",
    "    c1 = SkyCoord(ra_use_fs*u.deg, dec_use_fs*u.deg)\n",
    "\n",
    "    for i in range(np.size(dens_p_g)):\n",
    "        weights_use = weights[use_for_speed]\n",
    "        weights_use[i] = 0\n",
    "        c2 = SkyCoord(ra_use_fs[i]*u.deg, dec_use_fs[i]*u.deg)\n",
    "        sep = c1.separation(c2)\n",
    "        dist_i = sep.arcmin*dist\n",
    "        K_weigt = np.sum(weights_use/(2*np.pi*500**2)*np.exp(-(dist_i**2)/(2*500**2)))/(np.sum(weights_use))\n",
    "        dens_p_g[i] = K_weigt\n",
    "    \n",
    "    h_use = 500 * np.sqrt(gmean(dens_p_g)/dens_p_g)\n",
    "\n",
    "    size = ((np.size(x_edg)-1) * (np.size(y_edg)-1))\n",
    "    val = np.zeros(size)\n",
    "    dec_use = np.zeros(size)\n",
    "    ra_use = np.zeros(size)\n",
    "    count = 0\n",
    "\n",
    "    sum_w = np.sum(weights[use_for_speed])\n",
    "    for_f = weights[use_for_speed]/(2*np.pi*h_use**2)\n",
    "    for i in range(0,np.size(x_p)):\n",
    "      for j in range(0,np.size(y_p)):\n",
    "        ra_use[count] = x_p[i]\n",
    "        dec_use[count] = y_p[j]\n",
    "        count = count + 1\n",
    "\n",
    "    h_val = 2*h_use**2\n",
    "\n",
    "    for i in range(0,size):\n",
    "        c2 = SkyCoord(ra_use[i]*u.deg, dec_use[i]*u.deg)\n",
    "        sep = c1.separation(c2)\n",
    "        dist_i = sep.arcmin*dist\n",
    "\n",
    "        val[i] = np.sum(for_f*np.exp(-(dist_i**2)/(h_val)))/sum_w\n",
    "        #print(100*(i+0.0)/size)\n",
    "\n",
    "      \n",
    "\n",
    "\n",
    "\n",
    "    S_0, xedges, yedges = np.histogram2d(dec_use,ra_use , bins=(y_edg, x_edg), weights = val)\n",
    "    DAT = S_0\n",
    "    header = 'Matix'\n",
    "    np.savetxt('map%03d'%g+'numb%03d'%n+'_'+name_field+'.cat', DAT, delimiter=\" \", fmt='%s',  header= header, newline='\\n')\n",
    "\n",
    "    print(g, 'mean of output  = ', np.mean(S_0), np.sum(S_0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The density map is used in plt_dens.ipynb to assign the desnsity to the galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=75&v=4)\n",
    "\n",
    "**Authors**: Steven Duivenvoorden \n",
    "\n",
    "For a full description of the database and how it is organised in to `dmu_products` please the top level [readme](../readme.md).\n",
    "\n",
    "The Herschel Extragalactic Legacy Project, ([HELP](http://herschel.sussex.ac.uk/)), is a [European Commission Research Executive Agency](https://ec.europa.eu/info/departments/research-executive-agency_en)\n",
    "funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme, Grant Agreement Number 607254.\n",
    "\n",
    "[Acknowledgements](http://herschel.sussex.ac.uk/acknowledgements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}