{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib notebook\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "hide_input": true, "scrolled": true }, "outputs": [ { "data": { "text/latex": [ "from Roskar 2014\n", "\$$\n", "P_J = (4 \\Delta x_{min})^2 \\frac{G}{\\pi \\gamma}\\rho^2\n", "\$$\n", "where $\\Delta x_{min} = l_{box} \\,/ \\,2^{lmax}$, $\\gamma = 5/3$ and also the equilibrium temperature is defined as:\n", " \n", "\$$\n", "T_{eq} = \\frac{5000}{\\sqrt n_H}\n", "\$$\n", "here number density of hydrogen is expressed in (H/cm^3).Now, it can be proved that the n_star of a hydro run is given by:\n", "\$$\n", "n_{star} = \\frac{2k_b T_{eq}}\n", " {G(4\\pi \\Delta x_{min})}\n", "\$$\n", "Most of the heavy lifting is inside the calculation of $n_H$ that is as follows\n", "\$$\n", "n_H = \\frac{2k_b M_{\\%}^{-1} (3n_H^*)^{-1}} {G(4 \\Delta x_{min})^2}\n", "\$$\n", "\n", "thi is done until $|n_H-n_H^*|/n_H$>0.0001 for convergence. notice that here the is a molecular gas so $m_{\\%} = (0.76m_p+0.24m_{He})$." ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%latex\n", "from Roskar 2014\n", "\$$\n", "P_J = (4 \\Delta x_{min})^2 \\frac{G}{\\pi \\gamma}\\rho^2\n", "\$$\n", "where $\\Delta x_{min} = l_{box} \\,/ \\,2^{lmax}$, $\\gamma = 5/3$ and also the equilibrium temperature is defined as:\n", " \n", "\$$\n", "T_{eq} = \\frac{5000}{\\sqrt n_H}\n", "\$$\n", "here number density of hydrogen is expressed in (H/cm^3).Now, it can be proved that the n_star of a hydro run is given by:\n", "\$$\n", "n_{star} = \\frac{2k_b T_{eq}}\n", " {G(4\\pi \\Delta x_{min})}\n", "\$$\n", "Most of the heavy lifting is inside the calculation of $n_H$ that is as follows\n", "\$$\n", "n_H = \\frac{2k_b M_{\\%}^{-1} (3n_H^*)^{-1}} {G(4 \\Delta x_{min})^2}\n", "\$$\n", "\n", "thi is done until $|n_H-n_H^*|/n_H$>0.0001 for convergence. notice that here the is a molecular gas so $m_{\\%} = (0.76m_p+0.24m_{He})$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# then \n", "# constants\n", "k_b = 1.38e-23 # J K^-1 \n", "m_He = 6.646468e-27 # kg\n", "m_p = 1.672e-27 # kg\n", "G = 6.67e-11 # m^3 kg^-1 s^-2 \n", "m_mg = (0.76 * m_p + 0.24*m_He)\n", "pctocm = 3.08567758e18" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Reproducing Roskar 2013 value of $p_J$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "With Delta_x = 35.285949707 as in Roskar 2013 we reproduce the value for p_j = n_h = 18.3080473755\n", "and Teq = 1168.55454942\n" ] } ], "source": [ "# reproducing Roskar 2013\n", "levelmax = 20\n", "Delta_x = (37e6 / 2.**levelmax)\n", "\n", "n_h, aux=5., 20.\n", "i=0\n", "while (np.abs(n_h-aux)/n_h ) > 1e-4:\n", " aux = np.copy(n_h)\n", " n_h = (2*np.pi*k_b*1e4*np.sqrt(0.3))/\\\n", " (G* m_mg*np.sqrt(aux)*(4*Delta_x*pctocm/100)**2)\n", " n_h /= (m_p*1e6)#\n", "\n", "Teq = 5000/np.sqrt(n_h)\n", " \n", "print \"With Delta_x = {0} as in Roskar 2013 we reproduce the value for p_j = n_h = {1}\".format(Delta_x,n_h)\n", "print \"and Teq = {0}\".format(Teq)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# My value\n", "for levelmax =17 and boxlength = 25 Mpc " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "the resulting values for my sim are: \n", "p_j = 4.863 m_h / cc\n", "n_star = 8.336 m_H / cc\n" ] } ], "source": [ "# my box\n", "levelmax = 18# max reached by Zoom DMO\n", "box_len = 25e6 # parsec\n", "Delta_x = (box_len / 2.**levelmax)\n", "n_h, aux=5., 20.\n", "i=0\n", "while (np.abs(n_h-aux)/n_h ) > 1e-4:\n", " aux = np.copy(n_h)\n", " n_h = (2*np.pi*k_b*1e4*np.sqrt(0.3))/\\\n", " (G* m_mg*np.sqrt(aux)*(4*Delta_x*pctocm/100)**2) # kg per cubic meter\n", " n_h /= (m_p*1e6) # H per cubic centimeter\n", " \n", "# so my n_star is\n", "n_star = (2*np.pi*k_b*1e4*np.sqrt(0.3/n_h))/\\\n", " (G*m_p*(4.*Delta_x*pctocm/100)**2)\n", " \n", " \n", "n_star /= (m_p*1e6)\n", "print \"the resulting values for my sim are: \" \n", "print \"p_j = {0:.3f} m_h / cc\".format(n_h)\n", "print \"n_star = {0:.3f} m_H / cc\".format(n_star)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "mstar=n_star*(1./(2.**levelmax))**3. /n_h" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "7.612191598560314e-16" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(n_star/n_h)*(1./(2.**levelmax))**3." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "7.612191598560314e-16" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mstar" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "190.73486328125" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Delta_x\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.12" }, "latex_envs": { "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 } }, "nbformat": 4, "nbformat_minor": 1 }