HalphaSB_vs_density.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy import constants as const\n",
    "from astropy import units as u \n",
    "\n",
    "# load constants\n",
    "c = 3e8            # speed of light m/s\n",
    "hp = 6.626e-27     # planck's constant erg s\n",
    "H_0 = 100.         # hubble's constant /h km/s/Mpc\n",
    "h = 0.73"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#  how does the emissivity vary with temperature?  Can we trust that the variation of emissivity for H-alpha is similar \n",
    "#  to that of Ly-alpha?\n",
    "\n",
    "def Halpha():\n",
    "    ## Number density /cm^3.  Load values that agree with H-alpha emission at redshift 0\n",
    "    ## Emissivity erg/s/cm^3. Took ly-alpha emis (since didn't calculate H-alpha) \n",
    "    ## and used the energy fraction from Table B1 to approximate the emis of H-alpha\n",
    "\n",
    "  #  n_H = 4*10**-5.7         # taken as a popular density from the density map of Fig 4 B&S 2010   log(rho/rhom)~1.3\n",
    "\n",
    "    n_Harray = 4*10**-7 * 10**np.arange(0,4.1,0.1)   # number density ranging from 4*10^-7 to 4*10^-4\n",
    "\n",
    "    # Box:      dot-dash,  dash  ,  solid  , white/dotted\n",
    "    # Size:     1 Mpc   , 0.5Mpc , 0.15Mpc , 0.5Mpc\n",
    "    # log Temp: 5 - 6 K , 4 - 5 K, 3 - 4 K ,  5 K   \n",
    "    # emis:     -24.3   , -23.0  ,-23 - -23.5, -23.5 , \n",
    "    # density:  2.0     , 1.5    , 3.5     , 2.0-2.5\n",
    "    densityfactor = np.array([2.0,1.5,3.5,2.5])\n",
    "    n_Hplot = 4*10**-7 * 10**densityfactor  # dot-dash, dash, solid, white/dotted ; AVERAGE within the boxes\n",
    "    emiscoeff = np.array([-24.3,-23.0,-23.,-23.3])  \n",
    "    emis_plot = 0.66/11.26 * np.array([10**emiscoeff[0]*n_Hplot[0]**2,\n",
    "                                       10**emiscoeff[1]*n_Hplot[1]**2,\n",
    "                                       10**emiscoeff[2]*n_Hplot[2]**2,\n",
    "                                       10**emiscoeff[3]*n_Hplot[3]**2])\n",
    "    \n",
    "    z = 0                    # redshift\n",
    "    wavlen = 656.3e-9        # wavlength of line m  Which line should I use?  What fraction is \"energy\"? --> ergs, so use Halpha \n",
    "    \n",
    "    L_plot = np.array([1.,0.5,0.15,0.5])   # dot-dash, dash, solid\n",
    "    vol_plot = L_plot**3\n",
    "    \n",
    "    print('Doing H-alpha...')\n",
    "    # return:  emissivity, redshift, wavelength, length scale, volume\n",
    "    return [z,wavlen,emiscoeff,n_Harray,emis_plot,L_plot,vol_plot, densityfactor]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doing H-alpha...\n"
     ]
    }
   ],
   "source": [
    "# Call the function with the constants we want, now we want to load the H-alpha values.\n",
    "z,wavlen,emiscoeff,n_Harray,emis_plot,L_plot,vol_plot,densityfactor = Halpha()\n",
    "\n",
    "# cosmic mean density\n",
    "rho_mean = 5.1*10**-6*((1+z)/3)**3\n",
    "\n",
    "# calculate the luminosity and angular diameter distances (equation from cosmology notes, for EdS universe)\n",
    "if z == 0: \n",
    "    DL = 15. #Mpc\n",
    "    DA = 15.\n",
    "else:\n",
    "    DL = 2*c*h*10**-3/H_0 * (1.+z) * (1. - 1/np.sqrt(1.+z))  # Mpc\n",
    "    DA = 2*c*h*10**-3/H_0 / (1.+z) * (1. - 1/np.sqrt(1.+z))  # Mpc\n",
    "\n",
    "# approx angular size of the emitter (use to calculate the surface brightness)\n",
    "solidangle_plot = L_plot/DA**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  4.70029946e-35   9.37833037e-35   9.37833037e-31   4.70029946e-33]\n",
      "[ 1.    0.5   0.15  0.5 ]\n",
      "[ 1.        0.125     0.003375  0.125   ]\n"
     ]
    }
   ],
   "source": [
    "print emis_plot\n",
    "print L_plot\n",
    "print vol_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[    3.8283098      1.90962057  1718.65851305    95.70774508]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATUAAAEUCAYAAABK76GXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYVNW1t99FQyOCggREQUQRBwhxQDBgMLRBGQyKUTPo\nVQQeJOhnriNqQoAWFCSJUxIUQSHOiqJJcIAGtBkuiJ82QzNcZ5ShkSnQjHZDr/tHDVZ1V3XXcKrq\n1Kn1Ps95qHPO7rP27urzYw9r7SWqimEYhleol+kKGIZhOImJmmEYnsJEzTAMT2GiZhiGpzBRMwzD\nU5ioGYbhKVwpaiLSTEReF5H1IrJORLpnuk6GYUTGbe9r/Uwar4XHgXdU9RoRqQ80znSFDMOIiqve\nV3Gb862INAVWqGr7TNfFMIzaceP76sbh56nAdhGZISIlIjJNRI7OdKUMw4iI695XN4pafaAL8ISq\ndgH2A/dltkqGYUTBde+rG+fUNgGbVPX/+89fp9ovSUTcNWY2jBxCVSXktM73Nd24rqemqluBjSJy\nhv/SJcDaCOXSdowdO9bT9m688UZPty+d9tLdtl69eqXVXqLvazpxY08N4HfAiyKSD3wBDMlwfTzN\nhg0bMl2FlOL19rkAV72vrhQ1VV0FdMt0PQzDqBu3va+uG366kYKCAk/bO+WUU9JqL93tGzx4cNps\nef27ywZM1GLA66J26qmnptWel3+fXv/usgETNcPzFBcXZ7oKRhoxUTPo1atXpqtgGI5homakfciU\nbrzePiMcEzXDMDyFiZrheWxOLbcwUTMMw1OYqBmex+bUnGHVqlWZrkJMmKgZNjwzamXPnj3cdNNN\nnHvuuZmuSkyYqBksXLgw01VIKSbaiVNUVETnzp15+umnyc/Pz3R1YsJEzTCMqKxYsYJNmzbRrVs3\nSkpKMl2dmHBlQLthOInNqSXOXXfdRcuWLRk0aBD162eHXGRHLQ3DyAj169dn6NChma5GXNjw0/A8\nNqdWN0VFRbz99tu1llm9enWaapMcJmqGxX7mMOXl5QwfPpy+ffsyZMgQduzYEXY/dMfbH/3oR5mo\nYtyYqBmen3PyevsSZe7cuXTu3Jlp06aRn5/PnXfeSbNmzcLKjBw5kpdffhkAEYn0GNdhomYYOci4\ncePo168fGzduDK5s3nfffTUWA37zm99QWFjI4cOHM1TT+DFRMzyPzanVpHfv3jRq1IiHHnqIpUuX\n8sMf/hCA7777jltvvZWDBw8C0LVrV1asWJE1K59gomYYOclPfvITvv76a+69994wwWrYsCHbt29n\n3LhxwWuNGzfORBUTJnvk1zASJNfn1KqqqqhXr2b/pWXLlgBs3bqVlStX0q9fPwAmT54c7KllI9ZT\nM2x45lHKy8u56aabuP3222stt2vXLm644YZgKsEWLVrQtm3bNNQwNZioGRb76UFCYzanTZvGli1b\nwu6XlZVRXl4OQKdOnZgyZQoNGjTIRFUdx0TNMDxEYEeNvn37Blc2P/74Y1q3bh1WbuLEiYwcOTJ4\nfvXVV9OmTZt0VzclmKgZnieX5tTGjBkT3FFj4sSJLF26lE6dOgG+lc0A48ePp7KyksrKykxVNWWY\nqBmGhxg7diwDBgyo4XdWUVHB2Wefzfr16wFo2rQp06dP98yQMxQTNcPz5NKcWvPmzZk9e3bQ7yxA\nfn4+t912GzNmzMhQzdKHuXQYFvuZhezZs4ddu3bVWuajjz5ixowZTJ48GYCbb745KZvZMlR1bU9N\nRPJEZIWIzM50XbyO1+ecvNa+QMzmNddcU2u5Tp06MW/ePBYtWgT4YjcTjd8sLS2le/fuUe+76X11\nragBtwHrAM10RQzDDQRWNvv168emTZvIy8urUea1114L7lB79NFHs3jxYi666KKEbVZWVvLAAw9w\n/vnn17XzrWveV1eKmoicBFwGPA1kx9YAhmvxwpzavHnzwnIFBGI2q3Po0CEGDx5MRUUFAK1atUq6\ndzZ69GgqKyujDl/d9r66UtSAR4GRQFWmK2IYbmDDhg3BXAErVqwIi9mcP39+cM+z66+/nj//+c9J\nrWpW7521a9eOBQsW8MQTT0T7EXe9r4FN4NxyAAOAyf7PBcDsCGXUMHKJqqoqffHFF7WysjLsemFh\noZ5zzjn67LPPOmJn9erV2qVLF8U3jNSbb75Zy8vLg/f9715c72u6Dzeufl4IXCEilwFHAceKyHOq\nOii0UGFhYfBzQUGB5yaD00lxcbH9/lyOiHDdddcBvo7Ijh07ggHp//jHP1i7dm1Sz6+srGTSpEmM\nGzeOyspK2rVrxzPPPENeXh4PP/xwbT8a0/uaVjKtqrUdQC+sp5ZyCgsLM12FlPL+++9nugoxM3fu\nXH355ZdrLbN06VLt0KGD7t+/35HvrnrvbMSIEWG9s1Co1lPTGN7XdB9unVMLJeOrKYaRagI7avTt\n25fhw4ezadOmsPshwkGPHj249NJLWbVqVVI2I82dzZ8/nyeffJJjjjkm0cdm/n3NtKomcmA9NUfx\nek/N7cyZM0fbtm2rgObn5+vEiRNrzJ2NGjVK//a3v9X42US/u3h6Z6FQS0/NLUc29NQMw7NMmjQp\nplwBN9xwA4899hiHDh1Kyl6KemeuwkTN8Dxu9lPr378/TZo0Ce6oEYjZPHz4MMOHD2fPnj0AnHnm\nmaxZs4ajjjoqYVuR/M5KS0vp3bu3I21xCyZqhsV+ZpCzzz6bb775pkbvrH79+uTl5fGHP/wheC1R\nQavN78wrvbMwMj3+TeTA5tSMLOTw4cO13i8rK9NZs2YFz8vLy/Wbb76p9WfqmlOry+8sXrA5NcMw\nAiubgwcPrrXcd999x29/+1vWrFkDwDHHHJNwroCc652FYKJmeJ5MzqmF5gqYOXMmX375Zdj9srIy\nduzYAUC7du147rnnamRJj5doc2c/+9nPknputmCiZhgpINTvLHRls3379mHlnnjiCW655Zbgef/+\n/TnppJMSspnLvbMwMj3+TeTA5tQMlzNq1KiofmcHDhwIfj548KCOGDEi7Fo8BObUnJ47iwZZMKeW\n8QokVGkTNUfJpjCibGHv3r36y1/+UtesWRN2/fDhw9qxY0f98MMPHbFTWFio48eP1wYNGiig7dq1\n0wULFjjy7Ehkg6jFNPwUkeYxHMlNBBgZw/J+Ok+TJk2YOXNmjVwBeXl5jB49mueffz5pG6WlpQBZ\nPXeWCm0Rn/jWafg7YEsdxeqralrSOouIhu7SYRhGeigsLERVHdsIMiXaEkt3DljpRBmnDmz46SgW\n+5kYc+fO1bZt2+oZZ5xR65xYRUWFdunSRd9+++2kbVafOyssLEzJ3Fk0cHj4mQptiXX1M3rGhfjK\nGEbWU31ls2nTpkG3jABvvvlmMOFJgwYNmDdvHv3790/YZrSVTSDbVzYd15aYRE1V64yijaWMYWQC\nJ+fU3nvvvbBcAYGYzepOsvn5+QwePJgDBw4AvnycTuYKyLa5s2ikQlvcuPOtkWYs9jN2du7cGfQ7\nmzFjRnAhQFV555136NevH3l5efz85z+nSZMmNGrUKGFb0Xaj9VoAuuPEMJ5tDbQJOa53ckyd4Dg8\nvokAw3CQWbNm1djv7MiRI9qzZ0999NFHHbERq99ZuudDcXBOLVXaEktPrRtwIxDYZvNM4AWnRNUw\nso2rrroK8HUIysrKaN26NfXq1WPGjBksW7YsqWdH6p1Nnz7dE0PNCKRGW2JU1BNCPh/vlFInofDx\n/xdj5CyJOBfPnTtXp06dWmuZVatWaZs2bXTXrl0J1iycRKICsrmnpinSllgXCraGfN6WtJIahksp\nLy9n+PDh9O3bl9/97nd89tlnYfdVlaoqX3rLs88+m0GDBrF69eqkbOZyzGZKtCVOVe3hpEonoe6x\n/9diGDES8DujllwBDz74oI4bN84xm6tWrUoqZjPbe2r6/TvtmLbEa7hvKhqUwC8gvm/CqBWL/VR9\n7LHHgsLSrVu3GjGbATZu3KgdOnTQvXv3JmWvoqJCx40bl3TMpodEzTFtsa2HDIv9BC6//HKaN29e\nI1dAVVUVQ4cOZetW3yjppJNOYt26dTRp0iTh+qxevZru3bszZswYz/mduQHzUzMMoH379mzYsKHG\nHFa9evU44YQTGDlyZDAIvUGDBgnZqKys5KGHHmL8+PG5sLKZMWIWNfG5Q69PYV0MIyUUFBSEnVdW\nVkYUpoCgbdu2jblz53LDDTcAMHbsWLZv355UHVavXs2QIUMoKSkB4Oabb2bSpEmeXwiIBae1Jd7h\n5ztOGTaMdBNY2fzFL34RmMeJiIgwcuRIPvjgAwAaNmyY1G6048ePp2vXrmF5NnNhZTNOHNOWmHtq\nqqoi8rGIXKCqHzpVAcNINcXFxVRUVDBs2DA2btxIfn4+a9eupXPnzsEyZWVlHDlyhJNOOomWLVvy\n6quv0qZNm6TslpaWMnjw4GDvbMSIEfzpT38yMauG09oSb0+tO7BMRL4UkVL/kZyTjpFxvBz7WV5e\nzuzZs4M7agR6TKGCBvDiiy8yZMiQYA+uV69ejmZy8loW9BTgmLbEu1DQN+RzoP/u2IZxRmaoPufk\nJaZMmcIjjzxCfn4+hYWFjBw5Mpg0eO/evUGRuf3229myZQv79u1LSnhs7ixhnNOWOH1Jfgkc6/88\nGngD6JIKv5U66hG7Y42R03z33Xc6aNCgGn5nVVVVeu655+q8efMcsRPJ72z+/PmOPLs2POSn5pi2\nxDv8HKOq5SLSE+gNTAeeTEhNoyAibUXkfRFZKyJrROS/nXy+kVvk5+czZMiQGrkCRISJEycyc+bM\npG1E8zvLhS2CHHxfHdOWeEXtiP/fAcA0VX0LyE/EcC1UAneo6g/xjbP/n4h0dNiG4THKy8tZuXJl\nrWU+/fRTBgwYwJEjvj/jfv368dRTTyVsM9LKZq7EbIbg1PvqmLbEK2qbRWQq8GvgbRE5KoFn1Iqq\nblXVlf7P+/D5r7R20obhLQJZ0AcMGMDu3btr3A/MGXbo0IH9+/cza9as4L1Ed6O1qAAfDr6vzmlL\nnOPexsDVwOn+8xOBPqkYY/uffwrwNdCk2vW45gGM2snW2M89e/boTTfdFIzZ7Nq1q3711VdhZf71\nr3/p7Nmzw36mqqoqYZtOxWw6hZvm1KK9r7EcTmpLXKufqrofmBVyXgaUxS+ldSMiTYDXgdvU9z9A\nGKEp8goKCjy9gpdqFi5cmHW/v+LiYgYNGhT0O6u+shmgefPmTJo0iZ49e9KsWTOOPfbYhG3m4spm\ncXFxTLGzdb2vdeGktrgy9lNEGuBr4Auq+s9IZUKFLPCLDz0PlLHzus8D19xSn1jOv/zyy6Df2ciR\nIzn++OOpX78+qspTTz3F6aefTu/evenZsyclJSWsWLGCiy++OCF7CxYsYMmSJTz44INUVlYycOBA\nBg4cyJAhQ1zz+0jVeejfSEFBAffffz/VieV9TSuJdO9SeeDzTXkOeLSWMvH1mY1ayda8n++8806N\n/c6qqqp0wIABOmbMGEdsJLIbbTrJ9PAzlvc13Ycbtx76CXA9cLGIrPAf/TJdKcN99O/fP9g727Bh\nA+Cb+H/qqac47bTTknp2Lu9GGyeue19dN/xU1SU4vKJqZC9FRUWUlpZy1113RS3zxRdfcOGFF7Jy\n5Upat25N69atGTRoUPB+6NA6FixmM3bc+L7WWRkRaS0ibUKO69NRMSN9uDH2MzQL+j333FPDB01V\nOXz4MOBz1bj77rtZu3ZtUjYtZjO9pExb6hqfAgPxhSyM9R8vZXrMjM2peZpYcgU88sgjescddzhm\ns/rc2YgRI1w1dxaNTM+pJXOkSltiNW4p8oy0MGXKlJhyBezYsUPPOuss3blzZ1L2KioqdPz48a7x\nO4uXbBY1TZG2WIo8w1UMHDiQE088sUauAFVlyJAhfPnllwD84Ac/YM2aNTRv3rzOZ0bzswpEBYwe\nPTqnowIySSq0Ja4JPhHp4YRRw4jGCSecwBdffMF9990X5kgrInTq1Ik77rgjeC0vLy8hGxaz6T6c\n1JZ4Vy0Sd8c2jGp89913Ea83atQIgB07dvDkk99v1HDnnXeGncdK6MqnxWy6Fse0xVVLsUZmiCUM\nxkkCuQIuueSS4I4ZkcjPz2fSpEnMmTMH8PXMWrdObG8D653lDiZqRlrzfgZ21Jg2bRoffvhh0Bcs\nQFlZGV988QUAxx57LK+99lqNvdDi5dVXX7XeWQ4Rr6iVpqQWhucJ9M4CuQK6detGSUkJ3bp1Cyv3\n1ltvce211wZ90Lp165Z0roCpU6da78z9OKYt8YpaTxE5FkBERovImyLSxanKGN7lpZdeYtq0aeTn\n59dY2fzPf/4TLDds2DD69u3L3r17k7JXWloaXNl87733rHfmfpzTljh9Skr9//YEivHtUrncSb+V\nGOsRn3ONUSvp8HU6fPiw3nLLLRH9znr06KGvvfaaI3Yi+Z2lI1dApsh2PzX9/p12TFvcuJ234UHy\n8vKYPHlymN9ZgIcffpi33347aRuhvbPKykpGjBhBaWlpwq4fRlrx7nbeRvpxMvazvLw8mNk8Gt98\n8w0FBQVBl44ePXowY8aMhG1azKYnyI3tvGupRwIdZyPVBGI2mzdvrmVlZVHLVVVV6ZVXXqnTp09P\n2uaqVauyMmbTKTw0/HRMW9IqRg7+AuL+MozUsWfPHh02bFhYzOZnn30WVubtt9/Wl156KXi+b98+\nx3MFeHnuLBpeETUnj4xXIKFKm6i5hkWLFtW5o4aq6sqVK7VFixa6devWpG1W753VtRtttiaWiQUT\ntZqHzYcZSdGwYUM2b94c9DsLjdl87bXX2L9/PwDnnHMO8+fP5/jjj0/YlkUFGLHgup1vjeziggsu\nYMGCBfTs2bNGJqfZs2ezePFi/vrXvwI+YUuUZDI5ZVumLCM5Yuqp+dNfISINRMTWxz1GsrGfBQUF\nwVwBn376afD6448/XiNiIF6sd+ZtUqEtsWznfQ8wRkQeAZoCU5wwbLiHWGI/i4qKGDt2bK1ltmzZ\nwoUXXsjnn38OwHHHHccNN9yQcL0CfmfJxmymO2DfiI1UaUssw8/l/qMSnw+JzcPlEOXl5dx9991M\nmzYNgEsvvZSePXsG76sqFRUVNGzYkDZt2vDggw/y6aef0qFDh4RtVlZWMmnSJMaNG0dlZSXt2rVj\n+vTpFuLkPVKjLXWtJABdgREh5zdmenUDW/10lGgraNVzBUyYMKHGyuaUKVP0xhtvdKwuue53Fi/Z\nvPqZKm3J+PJrgr+M+L4Jo1YivRjPPvtsUFi6du0aNVfA3r179bzzztPNmzcnVQfzO0uMbBa1VB22\nnbcRkSuuuIL27dszYcIEli1bFran2bBhwygt9e0U06RJEz766KOEN2+E6HNnvXv3TrodYHNq2YCT\n2hKvS4dt5+1BIsV+NmvWjHXr1tGwYcMa93r06MFtt93Ge++9B0C9eolNhVRWVvLQQw8xfvx4mzsz\nbDtvwzkuuOCCiNcDgrZr1y7+/Oc/B4b+DB06lJdffjkpm+nMFWB+armFiVoOE8iC3qNHDyoqKqKW\nO/roo/nHP/7BzJkzAV9mp1atWiVk0/zOjFQjgf99Yyos0lpVt6SwPrHWQwsLC2tc79WrV8T/lYuL\niyP6YuVy+aKiIoYNG8bGjRvp3bs3F110UY3yHTt25Fe/+hUAa9asoWnTprRt29ax+uzbt48xY8bU\nEDOn21tcXBz2c274/Wdr+YsvvhhVlRo3k8RRbcn0SkWkA+gH/C/wGXBvhPtxrNcYoUTaUWPt2rU1\nys2cOVPPOussPXjwYNI2I61spjMLugW0OwcRVj/rel/TfcQrNt2AN4EV+BIllAKrHa0Q5AGfA6cA\nDYCVQMdqZZL5XnKaF198MeqOGtu3bw/bDuihhx7S7du3J2Uv3h01jPjItKjF8r7GcjipLfEa/hS4\nAmjvb8QpwCmJGK7FRg9gTsj5fcB91cok/+3kKFVVVXrPPfeE+Z0FejJ9+vTRqVOnOmIn072zXMEF\nolbn+xrL4aS2xLtQsF1V/62qX6rqhsAR5zPqog2wMeR8k/+a4QAiwqRJk8JyBQTmTx5++GEWLVqU\ntA2nYjadwvzUUopT76tj2hKvqN0vIs+IyLUicrX/uCoRw7UQ+8qFEZXy8nLef//9Wsts27aNrl27\nBs87d+7M888/n7DNSLkCbGXT8zj1vjqmLfE6394InOn/uaqQ628kYjwKm4HQ7LVt8al/GKGrnwUF\nBeaLFEJgZXPnzp2UlpbSvn37iOWOP/74pPY4C6W0tJTBgwcH9zsbMWIEf/rTn1whZva3kTjFxcV1\n9XRjel9jwDltiXPc+wl+N5BUHf5GfYFvTJ2PLRTETKSVzU8++SSsTFFRkU6ZMiV4fujQoaTmZXIt\nz6bbcMGcWp3vayyHk9oS7/BzKdApbuWMA1U9DNwKzAXWAa+q6vpU2vQCS5cupXPnzjz99NNhWdDP\nOOOMsHLt27dn1KhRfP311wARw6BiJRAVUD3PplMxm05hc2qpw8H31TFtiXf42QNYKSJfAd/5r6mq\nnu1EZQKo6rvAu04+0+scd9xxbNu2jW7dujFjxoywAPRXX32VSy+9lObNm3PaaaexePFiTj755OD9\nePN+RorZfOaZZ1wnZkZ6cOh9dUxb4hW1vkB1b2Kb2HcBHTt2pLi4mK5du9bIFbB06VJmz57NCy+8\nECwbSjxzTsnkCsgUNqeWFTimLfEOP8cBu/X75dY9QO17PBtpo3v37sFcAYGtgQAmTJjAZZddFpi7\nSAiL2TRSjGPaEq+onaOquwMnqvofoEsiho3EKCoq4vbbb6+1zM6dO+nduzcrV64EoHHjxlx33XWI\nJBayl84dNVKBzallBY5pS7yiJiLSPOSkOb4wCSPFBHbU6Nu3L48//jjvvhs+haGqHDx4EIAWLVrw\n+OOPBxcDEsV6Z0YacU5b4lx2HYRv6XU88ID/8yAnlmHjrEdMy89eoXqugEhZ0J9//nm94oorwmI3\nk2H16tUWs5kFZNqlw6nDSW1JxPgPgd/hW8btlIoGxlCHxL6RLGTmzJlhfmfRcgUcOnRIe/bsqV98\n8UXcNkJ3sYjkd2Yxm+7FK6KmDmpL2gXJocbH+11kLQcOHNAf/ehHEXtnw4cP1yVLlgTPE+2lBV6M\n6r0zr2Rysq2HnCOVoubUEa9Lh5FmGjVqRElJSQ03DYC+fftyxx13sHz5ckQk4YUAgAceeCAsz6b5\nnRnZim3n7SLKy8sjXg8I2p49eygsLKSqyhcad9VVV/Huu+8mJWarV68GCEYFOJ3JyQ2Yn1puYaLm\nAgIrm+eddx779u2LWq5JkybMmTMnmC0d4Ac/+EFCNkNXNgFb2TQ8Q1zDTxE5CrgaX/Bq4GdVVcc5\nXK+cITRXQH5+PkuXLqVPnz7B+2VlZWzevJmuXbuSl5fHCy+8kFS8JtSMCgDfLhteFbPqOQoM9+Gk\ntsTbU/sXvt0pK4F9/mN/vEaNcL+zjRs3Bn3BQgUNfAlPrr76avbs2QNAhw4daNu2baRH1kk0v7Ne\nvXp5VtCMrME5bYlnVQFYk+mVDfXI6ufs2bODfmcTJkwIW9ncunVr2Pnf//533bp1a1L2zO/Mm3hl\n9dNJbYnX8FTg7FQ0Ks56xPlVuJPx48dH9Du7+uqrddKkSY7YsFwB3sZDouaYtsQ0pyYigejoPGBI\nqrceyhX++Mc/Bj8fOXKEvDxfVMhf/vIX/vjHP6KqSa9sZtuOGqnA5tTcSyq0JdaFgsurnYdu9+B4\nYlMvUV5ezuLFi/n5z38etczu3bvp2bMnxcXFtGjRglNOOSW4TVAiRNrvbPr06VkTgG7kFI5rS0wL\nBfr9diDf4luheAx4FLgK2JqI4VygqKiIzp0784tf/CJsK6DqNGvWjL59+/LMM88kbdNtmZzcgPXS\n3EtKtCXOce9rwDPAxcDPgKeB11Ixxq6jHnHMAqSfSLkC1q9fH1bmvffeC5s3q6ioSCoYPZmYTS+H\nEXkdD82pOaYt8RpeF8u1VB9uFrUPPvigzh01VFW3bNmiLVu2rCF2iZBszGa6X4x042XR9pCoOaYt\n8fqplYhIj8CJiHQHPo7zGZ7mxBNPZPfu3XTr1o2SkhLuu+++YJjTK6+8wubNm4Plli1bxplnnpmw\nrUh5NufPn8+TTz6Zc4sBRtbjmLbEG9DeFfgfEdmIr1dwMvCJfwVD1VZBOfnkk1m0aBGdO3euEYS+\nfv16XnjhBWbPno2IcNpppyVsp/rKppvybLoNm1PLChzTlnhFrV+c5XOSc889F/AN7UtKSjj//PMB\nGDVqFP/85z+TerZlcjI8imPaEtfwU/0rFdEOpyqVDRQVFTF06NDA2D8ie/fu5fLLL2fx4sUA5Ofn\n86tf/crxXAEmaLVjOQrcj5PaEpOoiUiJE2W8QGjM5owZM3j99dfD7qtqcKeNY489lqlTp7Jt27ak\nbKY6V0C8eT8NwylSoi0xrkwcBErrOL5JxapIlPrUuUqTCmLJFTBr1izt1auXHjlyxBGbq1atsphN\nIyrZvvqZCm2JdU6tY91FOBzjs7KSt956i8sv9zk/R8qCHmDgwIFMmzaNTz75pEbS4HiwqAAjR3Be\nW5xU3XQdZKCnVlFRoT169IjYO7v11lv13XffDZ4nm9HJemfOYn5qzoHlKPAODRo0YPHixcGg81Cu\nvPJKbr/9dvr06UO9evUSXgiw3plhJI+rtvMWkT+LyHoRWSUib4hI00zUY9euXRGvBwRt37593Hvv\nvRw+7OsV9+7dmyVLllCvXuK/zmzPgu5mzE8tPYjIf/nf3dUi8j8ikhG/1YTfQhFpIcnsixOZIuCH\nqnoO8Cnwe4efXyuBlc3OnTtHFTaAxo0bU1JSwiOPPBK81rRpYvrrhizo5vJgOMSXwE/V5yg7Ht8e\naXGTrLbE6tLRQ0SK/b2nLiKyBlgDbBOR/okar46qzlPVKv/pcuAkp55dF4EdNZ5++ml27tzJ0qVL\nw+6XlZUF/c1EhBkzZnDttdcmZdMtO2osXLgwrfbSjYl2elDVZaq6x38a0/ubEm2JZeINXwxWH+CX\nwG6gu//AXStGAAAPZElEQVT6WcDKVEz2AbOB66Lci3N6MzqRdtSItBvtsmXLtFWrVrpt27akbbot\nC7oFtGcvbl0oAO4GpsZQznFtiVVgVoZ8Xl/t3oq4DMI8IvuiXB5SZhQwq5ZnJP3lBFi0aFFUv7Oy\nsjI9dOhQ8HzGjBm6ZcuWpOy5MQu610XNy7hR1PBtH7QOOC6Gso5pS+CIdfUzNBboUIw/E/lBqpfW\ndl9EBgOXAbXG/hQWFgY/FxQUJDwZfNFFF/HYY49xySWX1PA7u/fee2nTpg0TJkwAYPDgwQnZAN/c\n2aRJkywLupFVFBcX1zp8F5FbgJvwacRlwPHANKCfqv4nBhOOacv3T4ytd3UE2Os/Dod83gscTkRN\no9jpB6wFWtRRLon/a2qnelan66+/PunoALf7nXm9p2bDT+eglp4avp01Psc/hIzlSIW2xLqdd56q\nHuM/6od8PkZVnfR1+xvQBJgnIitE5AmnHlxeXl4jTrM6+/fvp3PnzmzcuBGAVq1a8fzzzyfsquGG\nlc1YsNhPwyHGAMcBT/rf3w/r+oGUaEuMatoNODHk/Ebg38BfgeaJqGkyB3H21AIxmyKiS5YsqbXs\nuHHjdOzYsXE9PxJu750Z3sBNPbVEjlRoS6yGVwQMAD8FyvAlSXgAeN3JRsZYn5i+gEgrm+vWrQsr\ns2TJEh01alTw/PDhw0nnCrA8m0a68ICoOa4tsRpeFfJ5MlAY6V66jlhE7eOPP44pV8CuXbu0TZs2\nWlJSUucz68KyoLsTm1NzjhSImuPaEutkUZ6INPB/vgR4P+SeK+NH27Zty6FDh4LzWaG5Al599VU+\n++wzAI477jiWLVsW3K02EQJzZ6G5Atw4d2YYLsRxbYn1h14GForIDuAAsBhARE7H5zDnOlq2bMnC\nhQs5/fTTa+QKKCsrY+jQoSxcuJB69erRtm3bhO1YFnT3Y7GfrsZ5bYmjm9gD+AXQOOTaGUAXJ7uj\nMdYlri5zVVVV2ALBkSNH9M0337S5Mz9eHp55nWwffmoKtCVmXwX1xXW9qar7Q659qqoZ3ca7qKiI\na665JrhjRiQOHTrEoEGDeOuttwCoV68eV155ZcJbBLklZtMpLPbTyCROa4urth6Kh9BcAbNmzeK5\n554Lu6+q7N7t6702atSI6dOns3///kiPiplIeTZt7sww3IUrJ/ljIeAkm5+fz/3338+gQYPC7hcV\nFfH73/+e5cuX06BBg6QdTEtLSxk8eLDl2cxCbE4tt8jantrGjRtrZEEPGY/Tp08fzjjjDNasWZOU\nHcuCbhjZRVw9NRG5C58PVmAySoE9wMequtLhutXKxIkTufvuu8NWNkeOHMn555/Ptddei4jwyiuv\nJGXDsqB7g+LiYuutuRwntSXentr5wAigNdAG+C3QH5gmIvfG+aykCPU7C/DrX/+awsLCWhcNYiFS\nzKaXe2cW+2m4AOe0Jc6l18VAk5DzJsAi4Giq7YWUygO/S8fBgwf11ltv1QMHDgSXnPft21f7mnQd\nWMymkU14waVDHdaWeHtqLYGKkPNKoJWqHsCpvZDi4KijjuLbb7/lgQceCF5r3LhxQs/Klh01DMOj\nOKYt8a5+vggsF5F/4hv7Xg68JCKN8e10mXYmT57MoUPJ6alFBXgbm1PLChzTlrhETVXHi8gc4EL/\npd+q6kf+z/8Vz7OS5euvv6Zdu3a0bNky4WdYnk3DcAdOaksifmqVQFXI54xQfZEgXqr7nVnvzLtY\nLy1rcERb4ppTE5HbgBfwjX+PB14Qkf9O1HgytGnTJqGfs6iAmlgYkZFpnNSWeBcKhgE/VtUxqjoa\n6I4v6UJWEIjZHD16tCdiNp3CYj8NF+CYtiQSUVAV5bNriRYVkMu9M8NwIY5oS7wTUzPwrVC8gW+F\n4kpgeqLG04HNnRk2p5YVOKYt8a5+PiIixUBP/6XBqroiEcOpJlKeTVvZNAx34qS2xCRqIrKP8KSj\nofdUVY9NxHiqML8zIxTzU3MvqdCWmERNVZvE++BMYH5niWGxn0amSIW2ZO1+atWx3lnieL0X4/X2\nGeFk7X5qASxm0zCMULJa1FavXu2pXAFGajA/tdwia0XNemeGYUQia0XNemdGrNicWm7hSlETkbtE\npEpEmkcrY70z57DhmeEkItJNRA6LyFWZsO86URORtsClwNe1lbPemXNY7KfhFCKSB0wC5vB9voG0\n4jpRAx4B7qmrkPXODMOV/A54HdieqQq4StREZCCwSVVXZ7ouhnewObX0ICJtgIHAk/5LESMFUk3a\nnW9FZB5wQoRbo4DfA31Ci6elUoZhOMFjwH2qqiIiZOj9Tbuoqeqlka6LSGfgVGCV7/fBScDHInKB\nqm6rXr6wsDD4uaCgwP43NqJisZ+JU1xcXOucpIjcwvf7njUFXvG/vy2A/iJSqar/TnU9w+rkT0fl\nOkTkK+B8Vd0V4Z66td7ZiNdfei+37/7772fs2LFpsyciqGqdPTARmQHMVtU30lCtMFw1p1YNU600\n4dUXPoDX22eE49qAdlVtn+k6GIaRGKo6JFO23dxTMwxHMD+13MJEzTAMT2GiZngem1PLLUzUDBue\nGZ7CRM2w2E/DU5ioGYbhKUzUDM9jc2q5hYmaYRiewkTN8Dw2p5ZbmKgZlvfT8BQmaobn55y83j4j\nHBM1wzA8hYma4XlsTi23MFEzDMNTmKgZnsfm1HILEzXDhmeGpzBRMyz20/AUJmqGYXgKEzXD89ic\nWm5homYYhqcwUTM8j82p5RYmaobFfhqewkTN8Pyck9fbZ4RjomYYhqcwUTM8j82p5RYmaoZheAoT\nNcPz2JxabmGiZtjwzPAUrhM1EfmdiKwXkTUiMinT9ckFLPbTcAoRKRCRFf73tzgTdXCVqInIxcAV\nwNmq2hn4S4arBKT/pUi3va+++iqt9tLdvpUrV6bNlte/u9oQkWbAZOBy//t7TSbq4SpRA24GJqpq\nJYCqbs9wfQDvi9qGDRvSai/d7du9e3fabHn9u6uD64BZqroJQFV3ZKISbhO104GfisgHIlIsIl0z\nXSHDMGLmdKC5iLwvIh+JyA2ZqET9dBsUkXnACRFujcJXn+NUtbuIdANmAu3TWT/De7isN+NlGgBd\ngN7A0cAyEflAVT9LZyVEVdNpr1ZE5F3gIVVd6D//HPixqu6sVs49lTaMHENVJfBZRG4BbvKfzgQa\nqmqh/97TwBxVfT2d9XObqP0WaK2qY0XkDGC+qp6c6XoZhlE3InIW8HegL9AQWA78WlXXpbMeaR9+\n1sF0YLqIlAIVwKAM18cwjBhR1f8VkTnAaqAKmJZuQQOX9dQMwzCSxW2rn2GISD8R+V8R+UxE7o1S\n5q/++6tE5LxU2vM7Fu7xOxeuEJE/JmFruoh86++VRivjZNtqtedw29r6V8DW+p0w/ztKOUfaF4s9\nh9t3lIgsF5GVIrJORCZGKedU++q052T7Qp6Z53/W7Cj3Hfv7dBRVdeUB5AGfA6fgW1VZCXSsVuYy\n4B3/5x8DH6TYXgHwb4fadxFwHlAa5b5jbYvRnpNtOwE41/+5CfBJir+7WOw51j7/8472/1sf+ADo\nmeLvry57jrbP/8w7gRcjPdfp9jl5uLmndgHwuapuUJ8z7ivAwGplrgCeBVDV5UAzEWmVQnsAEuFa\n3KjqYuA/tRRxsm2x2APn2rZVVVf6P+8D1gOtqxVzrH0x2gOH2ue3c8D/MR/ff4i7qhVx+vuryx44\n2D4ROQmfcD0d5bmOts9J3CxqbYCNIeeb/NfqKnNSCu0pcKG/u/2OiHRK0Fai9Um0bbGQkraJyCn4\neojLq91KSftqsedo+0SknoisBL4F3teaE+KOti8Ge05/f48CI/FN+Eci3X+fMeNmUYt1BaP6/yKJ\nrnzE8nMlQFtVPQf4G/DPBG3FilNtiwXH2yYiTYDXgdv8PagaRaqdJ9W+Ouw52j5VrVLVc/G9yD8V\nkYJIVar+Yym051j7RGQAsE1VV1B77y+df58x42ZR2wy0DTlvi+9/g9rKnOS/lhJ7qro3MAxQ1XeB\nBiLSPEF78dYnmbbVidNtE5EGwCzgBVWN9II52r667KXqu1PVPcDbQPWQvpR8f9HsOdy+C4ErROQr\n4GXgZyLyXLUyaf37jAc3i9pHwOkicoqI5AO/Bv5drcy/8fuyiUh3YLeqfpsqeyLSSkTE//kCfC4x\nkeY2nMDJttWJk23zP+cZYJ2qPhalmGPti8Wew+1rIb4dKRCRRsClwIpqxZxsX532nGyfqv5BVduq\n6qnAb4D3VLW6z2ha/z7jwW3Ot0FU9bCI3ArMxTcx+oyqrhdf1AGq+pSqviMil4kvnGo/MCSV9vBt\npXKziBwGDuD7whNCRF4GegEtRGQjMBbfqqvjbYvFHg62DfgJcD2wWkQCL98fgJMD9hxuX532cLZ9\nJwLPikg9fB2D51V1Qar+NmOxh7Ptq45CMOInJX+fTmLOt4ZheAo3Dz8NwzDixkTNMAxPYaJmGIan\nMFEzDMNTmKgZhuEpTNQMw/AUJmqGYXgKEzXDMDyFiZpRAxGJFHzu5POniMiFqbRh5C4makYkUh1m\n8mNgWYptGDmKiZoRFRG5U0RK/cdtIddHi2/b88Ui8pKI3BXHMzsCn6rF5xkpwrUB7UZmEZEuwGB8\nOwLXA5aLyEJ8QfBXAWfj24W1BN8OJ9Ge81tgO9BeVf8C9AfereW+YSSF9dSMaPQE3lDVg6q6H3gD\nX56DC4F/qmqFfyPG2UTZSFBEbgS2qOobQB//5T7AnFruG0ZSWE/NqA2J83N1rgUu82+Z8wP/XmDN\nVHVrpPtOVNgwrKdmRGMxcKWINBKRxsCVwCLgf4DLRaShf/vsnxNhYcH/M6KqVcDlwL+Ai4H3a7lv\nGEljPTUjEqqqK0TkH8CH/mvTVHUVgIj8G18W7m+BUmBPhGf8GDgoIgOBjsAE4DFgZpT7D/qf/WN8\n826zgF8C7wCdgKVAd2AbsM7/+VxgHHCO/2feBDqp6ouO/BaMrMREzaiBqh7r//dRfFmFqvMXVb1f\nRI4GFgIfRyhzIfCwPzXfvwBEpAdwW7T7fr7CJ5Z78eVhrQIOAfcBf8W3y+ptwF3Aqf57W4CvgUp8\nu74aOYwNP41EmOrfNvtj4PVAzs1qnIov6W4QVT1fVY9Eu++nB7AA3xbdC/Ctsu71fwY4Gp/YdQUa\n4uvl/RiYD5wPHEy8WYYXsO28DcPwFNZTMwzDU5ioGYbhKUzUDMPwFCZqhmF4ChM1wzA8hYmaYRie\nwkTNMAxPYaJmGIanMFEzDMNT/B/lDzwAxcGrNAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e1df510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "emis1 = 0.66/11.26 * 10**emiscoeff[2] * n_Harray**2\n",
    "emis2 = 0.66/11.26 * 10**emiscoeff[1] * n_Harray**2\n",
    "emis3 = 0.66/11.26 * 10**emiscoeff[0] * n_Harray**2\n",
    "emis4 = 0.66/11.26 * 10**emiscoeff[3] * n_Harray**2\n",
    "\n",
    "SB_plot1 = emis1 * vol_plot[2] / (4.*np.pi*DL**2) * wavlen * (1+z)/(hp*c) * 3.1e24 / solidangle_plot[2]\n",
    "SB_plot2 = emis2 * vol_plot[1] / (4.*np.pi*DL**2) * wavlen * (1+z)/(hp*c) * 3.1e24 / solidangle_plot[1]\n",
    "SB_plot3 = emis3 * vol_plot[0] / (4.*np.pi*DL**2) * wavlen * (1+z)/(hp*c) * 3.1e24 / solidangle_plot[0]\n",
    "SB_plot4 = emis4 * vol_plot[3] / (4.*np.pi*DL**2) * wavlen * (1+z)/(hp*c) * 3.1e24 / solidangle_plot[3]\n",
    "\n",
    "# Plot the SB as a function of the density n_H\n",
    "fig, ax1 = plt.subplots(figsize=(4,4))\n",
    "ax1.set_xlabel(r'log $\\rho$/$\\rho_{mean}$')\n",
    "ax1.set_ylabel('log SB [ph cm$^{-2}$ s$^{-1}$ sr$^{-1}$]')\n",
    "#ax1.legend(loc=4)\n",
    "#ax1.get_legend()\n",
    "\n",
    "ax1.plot(np.arange(0,4.1,0.1),np.log10(SB_plot1),'k-',linewidth=2)                          # logT ~ 3\n",
    "ax1.plot(np.arange(0,4.1,0.1),np.log10(SB_plot2),'k--',linewidth=2)        # logT ~ 3\n",
    "#ax1.plot(np.arange(0,4.1,0.1),np.log10(SB_plot3),'k:',linewidth=2)                          # logT ~ 6\n",
    "ax1.plot(np.arange(0,4.1,0.1),np.log10(SB_plot4),'k:',linewidth=2)                          # logT ~ 6\n",
    "\n",
    "#plot vertical and horizontal lines to indicate actual SB and density\n",
    "SB_plot = emis_plot * vol_plot / (4.*np.pi*DL**2) * wavlen * (1+z)/(hp*c) * 3.1e24 / solidangle_plot\n",
    "print(SB_plot)\n",
    "plt.plot([0.,4.0],np.log10([SB_plot[2],SB_plot[2]]),'k-',linewidth=1,color='gray')\n",
    "plt.plot([0.,4.0],np.log10([SB_plot[1],SB_plot[1]]),'k--',linewidth=1,color='gray')\n",
    "#plt.plot([0.,4.0],np.log10([SB_plot[0],SB_plot[0]]),'k:',linewidth=1,color='gray')\n",
    "plt.plot([0.,4.0],np.log10([SB_plot[3],SB_plot[3]]),'k:',linewidth=1,color='gray')\n",
    "plt.plot([densityfactor[2],densityfactor[2]],[-6,6],'k-',linewidth=1,color='gray')\n",
    "plt.plot([densityfactor[1],densityfactor[1]],[-6,6],'k--',linewidth=1,color='gray')\n",
    "#plt.plot([densityfactor[0],densityfactor[0]],[-6,6],'k:',linewidth=1,color='gray')\n",
    "plt.plot([densityfactor[3],densityfactor[3]],[-6,6],'k:',linewidth=1,color='gray')\n",
    "\n",
    "ax3 = ax1.twinx()\n",
    "ax3.set_yticks([-6,-4,-2,0,2,4,6])\n",
    "ax3.set_yticklabels([-6,-4,-2,0,2,4,6])\n",
    "ax3.set_ylabel('log SB [ph cm$^{-2}$ s$^{-1}$ sr$^{-1}$]')\n",
    "#convert to mag/arsec^2\n",
    "\n",
    "\n",
    "plt.savefig('FigHalphaSBBoxes.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}