Page MenuHomec4science

GeneratingVelocities.html
No OneTemporary

File Metadata

Created
Fri, Feb 7, 19:21

GeneratingVelocities.html

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Generating velocities &mdash; pNbody v4 documentation</title>
<link rel="stylesheet" href="../_static/default.css" type="text/css" />
<link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: '../',
VERSION: '4',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="../_static/jquery.js"></script>
<script type="text/javascript" src="../_static/underscore.js"></script>
<script type="text/javascript" src="../_static/doctools.js"></script>
<link rel="top" title="pNbody v4 documentation" href="../index.html" />
<link rel="up" title="Generating initial conditions" href="InitialConditions.html" />
<link rel="next" title="Selectig particles" href="Selection.html" />
<link rel="prev" title="Generating mass profiles" href="GeneratingMassProfiles.html" />
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="../py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="../np-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="Selection.html" title="Selectig particles"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="GeneratingMassProfiles.html" title="Generating mass profiles"
accesskey="P">previous</a> |</li>
<li><a href="../index.html">pNbody v4 documentation</a> &raquo;</li>
<li><a href="InitialConditions.html" accesskey="U">Generating initial conditions</a> &raquo;</li>
</ul>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="generating-velocities">
<h1>Generating velocities<a class="headerlink" href="#generating-velocities" title="Permalink to this headline"></a></h1>
<p>Setting initial velocities may be trivial or very complex, depending
on the aim. Of course, its allways possible to directely set the
variable &#8216;&#8217;nb.vel&#8217;&#8217; to some values. However, in galactic dynamics, this
is usually nor the right way to do.</p>
<p><strong>pNbody</strong> offers different methods to set the velocities in order to ensure
the system to be in a resonably good equilibrium.</p>
<table border="1" class="docutils">
<colgroup>
<col width="60%" />
<col width="12%" />
<col width="28%" />
</colgroup>
<thead valign="bottom">
<tr><th class="head">method name</th>
<th class="head">geometry</th>
<th class="head">method</th>
</tr>
</thead>
<tbody valign="top">
<tr><td><tt class="xref py py-func docutils literal"><span class="pre">NbodyDefault.Get_Velocities_From_AdaptativeSpherical_Grid()</span></tt></td>
<td>spherical</td>
<td>jeans equation+self-adaptative grid</td>
</tr>
<tr><td><tt class="xref py py-func docutils literal"><span class="pre">NbodyDefault.Get_Velocities_From_Spherical_Grid()</span></tt></td>
<td>spherical</td>
<td>jeans equation</td>
</tr>
<tr><td><tt class="xref py py-func docutils literal"><span class="pre">NbodyDefault.Get_Velocities_From_Cylindrical_Grid()</span></tt></td>
<td>cylindrical</td>
<td>jeans equation</td>
</tr>
</tbody>
</table>
<p>Some examples showing how to set the initial velocities for spherical or axi-symetric models are provided in the
example directory. To get it, type:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">pNbody_copy</span><span class="o">-</span><span class="n">examples</span>
</pre></div>
</div>
<p>by default, this create a directory in your home <tt class="docutils literal"><span class="pre">~/pnbody_examples</span></tt>.
The scripts are in the <tt class="docutils literal"><span class="pre">~/pnbody_examples/ic</span></tt> directory.</p>
<div class="section" id="spherical-coordinate">
<h2>Spherical coordinate<a class="headerlink" href="#spherical-coordinate" title="Permalink to this headline"></a></h2>
<p>In spherical coordinates, supposing that the velocities are isotropics. the Jeans equations is reduced to one integral,
giving the value of the velocity dispertion of a component <img class="math" src="../_images/math/34857b3ba74ce5cd8607f3ebd23e9015908ada71.png" alt="i"/> :</p>
<div class="math">
<p><img src="../_images/math/a5f2e7938b43b33c24e120d4766d15391227b717.png" alt="\sigma_{i}^2(r) = \frac{1}{\rho_{i}(r)}\int_r^\infty\! dr' \,\rho_{i}(r')\, \partial_{r'} \Phi(r')." /></p>
</div></div>
<div class="section" id="example-set-a-plummer-sphere-to-the-jeans-equilibrium">
<h2>Example : Set a Plummer sphere to the jeans equilibrium<a class="headerlink" href="#example-set-a-plummer-sphere-to-the-jeans-equilibrium" title="Permalink to this headline"></a></h2>
<p>First, generate the Plummer sphere:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">pNbody</span> <span class="kn">import</span> <span class="n">ic</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">n</span> <span class="o">=</span> <span class="mi">100000</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">a</span> <span class="o">=</span> <span class="mf">1.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">rmax</span> <span class="o">=</span> <span class="mf">10.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nb</span> <span class="o">=</span> <span class="n">ic</span><span class="o">.</span><span class="n">plummer</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="n">eps</span><span class="o">=</span><span class="n">a</span><span class="p">,</span><span class="n">rmax</span><span class="o">=</span><span class="n">rmax</span><span class="p">,</span><span class="n">name</span><span class="o">=</span><span class="s">&#39;plummer.dat&#39;</span><span class="p">,</span><span class="n">ftype</span><span class="o">=</span><span class="s">&#39;gadget&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Then, we need to create a 1D grid wraping the model:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">grmin</span> <span class="o">=</span> <span class="mi">0</span> <span class="c"># grid minimal radius</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">grmax</span> <span class="o">=</span> <span class="n">rmax</span><span class="o">*</span><span class="mf">1.05</span> <span class="c"># grid maximal radius</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nr</span> <span class="o">=</span> <span class="mi">128</span> <span class="c"># number of radial bins</span>
</pre></div>
</div>
<p>Now, in order to decrease the noise in outer regions, it is usefull to
use a non linear grid. In this purpose, we need to define a function that
defines the new nodes of the grid. Here, we set it logarithmic. The inverse function
is also needed:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">rc</span> <span class="o">=</span> <span class="n">a</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">g</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">r</span><span class="p">:</span><span class="n">log</span><span class="p">(</span><span class="n">r</span><span class="o">/</span><span class="n">rc</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="c"># the function</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">gm</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">r</span><span class="p">:</span><span class="n">rc</span><span class="o">*</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="n">r</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span> <span class="c"># and its inverse</span>
</pre></div>
</div>
<p>Now it possible to invoque the magic function:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">eps</span> <span class="o">=</span> <span class="mf">0.01</span> <span class="c"># gravitational softening</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nb</span><span class="p">,</span><span class="n">phi</span><span class="p">,</span><span class="n">stats</span> <span class="o">=</span> <span class="n">nb</span><span class="o">.</span><span class="n">Get_Velocities_From_Spherical_Grid</span><span class="p">(</span><span class="n">eps</span><span class="o">=</span><span class="n">eps</span><span class="p">,</span><span class="n">nr</span><span class="o">=</span><span class="n">nr</span><span class="p">,</span><span class="n">rmax</span><span class="o">=</span><span class="n">grmax</span><span class="p">,</span><span class="n">phi</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">g</span><span class="o">=</span><span class="n">g</span><span class="p">,</span><span class="n">gm</span><span class="o">=</span><span class="n">gm</span><span class="p">,</span><span class="n">UseTree</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span><span class="n">ErrTolTheta</span><span class="o">=</span><span class="mf">0.5</span><span class="p">)</span>
</pre></div>
</div>
<p>If <tt class="docutils literal"><span class="pre">UseTree=True</span></tt> the function uses a <tt class="docutils literal"><span class="pre">treecode</span></tt> to compute the potential at each node.
A new <tt class="docutils literal"><span class="pre">nb</span></tt> object with velocities computed from the Jeans equation is returned. We also get the potential in each node
and some statistics. The potential <tt class="docutils literal"><span class="pre">phi</span></tt> is usefull if we need to run the method for another component.</p>
<p>Using the <tt class="docutils literal"><span class="pre">stats</span></tt> variable, it is possible to plot some interesting values used during the computation, like the velocity dispertion as a function of the radius:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">import</span> <span class="nn">pylab</span> <span class="kn">as</span> <span class="nn">pt</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;r&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;sigma&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<img alt="../_images/r-sigma.png" src="../_images/r-sigma.png" />
<p>or the density profile:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;r&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;rho&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">loglog</span><span class="p">()</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<img alt="../_images/r-rho.png" src="../_images/r-rho.png" />
<p>It is possible to check the velocity dispertion along the line of sight directely
by mesuring it using a grid object:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">pNbody</span> <span class="kn">import</span> <span class="n">libgrid</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">G</span> <span class="o">=</span> <span class="n">libgrid</span><span class="o">.</span><span class="n">Spherical_1d_Grid</span><span class="p">(</span><span class="n">rmin</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">rmax</span><span class="o">=</span><span class="n">rmax</span><span class="p">,</span><span class="n">nr</span><span class="o">=</span><span class="mi">128</span><span class="p">,</span><span class="n">g</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">gm</span><span class="o">=</span><span class="bp">None</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">r</span> <span class="o">=</span> <span class="n">G</span><span class="o">.</span><span class="n">get_r</span><span class="p">()</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">sigma</span> <span class="o">=</span> <span class="n">G</span><span class="o">.</span><span class="n">get_SigmaValMap</span><span class="p">(</span><span class="n">nb</span><span class="p">,</span><span class="n">nb</span><span class="o">.</span><span class="n">vz</span><span class="p">())</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">r</span><span class="p">,</span><span class="n">sigma</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s">&quot;radius&quot;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s">&quot;velocity dispertion&quot;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<img alt="../_images/r-sigma-mes.png" src="../_images/r-sigma-mes.png" />
</div>
<div class="section" id="cylindrical-coordinate">
<h2>Cylindrical coordinate<a class="headerlink" href="#cylindrical-coordinate" title="Permalink to this headline"></a></h2>
<p>As in spherical geometry, we use the Jeans equations to get an approximation of the velocity dispertions.</p>
<div class="section" id="vertical-velocity-dispersion">
<h3>Vertical velocity dispersion<a class="headerlink" href="#vertical-velocity-dispersion" title="Permalink to this headline"></a></h3>
<p>The resolution of these equations in cylindrical coordinates gives directely the dispersion along the z axis:</p>
<div class="math">
<p><img src="../_images/math/a6deaa7f455bbae6d9a3ed36b158801eed379434.png" alt="{\sigma_z}_i^2 = \frac{1}{\rho_i(z)}\int_z^\infty\! dz' \,\rho_i(z')\, \partial_{z'} \Phi(z')." /></p>
</div><p>and the mean azimuthal velocity:</p>
</div>
<div class="section" id="radial-velocity-dispersion">
<h3>Radial velocity dispersion<a class="headerlink" href="#radial-velocity-dispersion" title="Permalink to this headline"></a></h3>
<ol class="arabic simple">
<li>if <tt class="docutils literal"><span class="pre">mode_sigma_r['name']='epicyclic_approximation'</span></tt>, we use the epicyclic approximation which hold:</li>
</ol>
<div class="math">
<p><img src="../_images/math/99c1f30d1f4c875d7fb755d4282085ea3818b1da.png" alt="\sigma_R^2 = \sigma_z^2 \frac{1}{\beta^2} \frac{\nu^2}{\kappa^2}" /></p>
</div><p>the parameter <img class="math" src="../_images/math/04777e262a979166a6295fb67e24ca75bc42f3b6.png" alt="\beta^2"/> is set by the variable <tt class="docutils literal"><span class="pre">mode_sigma_r['param']</span></tt> and takes usually a value around 1.</p>
<ol class="arabic simple" start="2">
<li>if <tt class="docutils literal"><span class="pre">mode_sigma_r['name']='isothropic'</span></tt>, the value is simply:</li>
</ol>
<div class="math">
<p><img src="../_images/math/26f488dbfe0e8985861ece799eaed8f37561717c.png" alt="\sigma_R^2 = \sigma_z^2" /></p>
</div><ol class="arabic" start="3">
<li><p class="first">if <tt class="docutils literal"><span class="pre">mode_sigma_r['name']='toomre'</span></tt>, we determine the velocity dispersion uses the Safronov-Toomre parameter <img class="math" src="../_images/math/9866e3a998d628ba0941eb4fea0666ac391d149a.png" alt="Q"/>:</p>
<div class="math">
<p><img src="../_images/math/ded125c9d4c060477de06e79b24b3d76661fc17e.png" alt="\sigma_R = \frac{3.36 Q G \Sigma_i}{\kappa}" /></p>
</div><p><img class="math" src="../_images/math/9866e3a998d628ba0941eb4fea0666ac391d149a.png" alt="Q"/> is set with the variable <tt class="docutils literal"><span class="pre">mode_sigma_r['param']</span></tt>.</p>
</li>
<li><p class="first">if <tt class="docutils literal"><span class="pre">mode_sigma_r['name']='constant'</span></tt>, the value is simply constant, given by the variable <tt class="docutils literal"><span class="pre">mode_sigma_r['param']</span></tt>:</p>
</li>
</ol>
<div class="math">
<p><img src="../_images/math/b7b72642c21b9dc0c18bc48566b5a528387547b2.png" alt="\sigma_R = cte" /></p>
</div></div>
<div class="section" id="azimuthal-velocity-dispersion">
<h3>Azimuthal velocity dispersion<a class="headerlink" href="#azimuthal-velocity-dispersion" title="Permalink to this headline"></a></h3>
<ol class="arabic simple">
<li>if <tt class="docutils literal"><span class="pre">mode_sigma_p['name']='epicyclic_approximation'</span></tt>, we use the epicyclic approximation which hold:</li>
</ol>
<div class="math">
<p><img src="../_images/math/a70b008211f747612127d906de127eb49abd6f82.png" alt="\sigma_\phi^2 = \sigma_r^2 \frac{1}{4} \frac{\kappa^2}{\Omega^2}" /></p>
</div><ol class="arabic simple" start="2">
<li>if <tt class="docutils literal"><span class="pre">mode_sigma_p['name']='isothropic'</span></tt>, the value is simply:</li>
</ol>
<div class="math">
<p><img src="../_images/math/ae1c33e8b543fea77262a8f055f39575e36e5a5e.png" alt="\sigma_\phi^2 = \sigma_z^2" /></p>
</div></div>
<div class="section" id="mean-azimuthal-velocity">
<h3>Mean azimuthal velocity<a class="headerlink" href="#mean-azimuthal-velocity" title="Permalink to this headline"></a></h3>
<p>Finally, the mean azimuthal velocity is also derived directely from the Jeans equations:</p>
<div class="math">
<p><img src="../_images/math/c6c6e55195f1df3be1847c9d7ae40ae6bfa4e70b.png" alt="\langle v_{\phi}^2 \rangle = R\,\partial_{R} \Phi(R,z) + \sigma_R^2 - \sigma_\phi^2 + \frac{R}{\rho_i} \partial_R \left( \rho_i \sigma_R^2 \right)," /></p>
</div></div>
</div>
<div class="section" id="example-set-an-exponnential-disk-to-the-jeans-equilibrium">
<h2>Example : Set an exponnential disk to the jeans equilibrium<a class="headerlink" href="#example-set-an-exponnential-disk-to-the-jeans-equilibrium" title="Permalink to this headline"></a></h2>
<p>Lest try to put an exponnential disk at the Jeans equilibrium.
First, we generate the disk:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">pNbody</span> <span class="kn">import</span> <span class="n">ic</span>
<span class="gp">&gt;&gt;&gt; </span><span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="o">*</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">n</span> <span class="o">=</span> <span class="mi">100000</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Hz</span> <span class="o">=</span> <span class="mf">0.3</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Hr</span> <span class="o">=</span> <span class="mf">3.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Rmax</span> <span class="o">=</span> <span class="mf">30.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">Zmax</span> <span class="o">=</span> <span class="mf">3.</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nb</span> <span class="o">=</span> <span class="n">ic</span><span class="o">.</span><span class="n">expd</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">Hr</span><span class="p">,</span><span class="n">Hz</span><span class="p">,</span><span class="n">Rmax</span><span class="p">,</span><span class="n">Zmax</span><span class="p">,</span><span class="n">irand</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span><span class="n">name</span><span class="o">=</span><span class="s">&#39;expd.dat&#39;</span><span class="p">,</span><span class="n">ftype</span><span class="o">=</span><span class="s">&#39;gadget&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Then, we need to set the parameters for the cylindrical grid:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">grmin</span> <span class="o">=</span> <span class="mf">0.</span> <span class="c"># minimal grid radius</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">grmax</span> <span class="o">=</span> <span class="mf">1.1</span><span class="o">*</span><span class="n">Rmax</span> <span class="c"># maximal grid radius</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">gzmin</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.1</span><span class="o">*</span><span class="n">Zmax</span> <span class="c"># minimal grid z</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">gzmax</span> <span class="o">=</span> <span class="mf">1.1</span><span class="o">*</span><span class="n">Zmax</span> <span class="c"># maximal grid z</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nr</span> <span class="o">=</span> <span class="mi">32</span> <span class="c"># number of bins in r</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nt</span> <span class="o">=</span> <span class="mi">2</span> <span class="c"># number of bins in t</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">nz</span> <span class="o">=</span> <span class="mi">64</span><span class="o">+</span><span class="mi">1</span> <span class="c"># number of bins in z</span>
</pre></div>
</div>
<dl class="docutils">
<dt>Set some functions used to distort the grid along the radius::</dt>
<dd><div class="first last highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">rc</span> <span class="o">=</span> <span class="n">Hr</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">g</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">r</span><span class="p">:</span><span class="n">log</span><span class="p">(</span><span class="n">r</span><span class="o">/</span><span class="n">rc</span><span class="o">+</span><span class="mf">1.</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">gm</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">r</span><span class="p">:</span><span class="n">rc</span><span class="o">*</span><span class="p">(</span><span class="n">exp</span><span class="p">(</span><span class="n">r</span><span class="p">)</span><span class="o">-</span><span class="mf">1.</span><span class="p">)</span>
</pre></div>
</div>
</dd>
</dl>
<p>Set some options on how to compute velocity dispertions:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_z</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;jeans&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="bp">None</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_r</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;toomre&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="mf">1.0</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_p</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;epicyclic_approximation&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="bp">None</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">mode_sigma_z</span><span class="p">,</span><span class="n">mode_sigma_r</span><span class="p">,</span><span class="n">mode_sigma_p</span><span class="p">]</span>
</pre></div>
</div>
<p>Set the gravitational softening:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">eps</span><span class="o">=</span><span class="mf">0.1</span>
</pre></div>
</div>
<p>And finally, lunch the magic function:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">nb</span><span class="p">,</span><span class="n">phi</span><span class="p">,</span><span class="n">stats</span> <span class="o">=</span> <span class="n">nb</span><span class="o">.</span><span class="n">Get_Velocities_From_Cylindrical_Grid</span><span class="p">(</span><span class="n">select</span><span class="o">=</span><span class="s">&#39;0&#39;</span><span class="p">,</span><span class="n">disk</span><span class="o">=</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span><span class="n">eps</span><span class="o">=</span><span class="n">eps</span><span class="p">,</span><span class="n">nR</span><span class="o">=</span><span class="n">nr</span><span class="p">,</span><span class="n">nz</span><span class="o">=</span><span class="n">nz</span><span class="p">,</span><span class="n">nt</span><span class="o">=</span><span class="n">nt</span><span class="p">,</span><span class="n">Rmax</span><span class="o">=</span><span class="n">grmax</span><span class="p">,</span><span class="n">zmin</span><span class="o">=</span><span class="n">gzmin</span><span class="p">,</span><span class="n">zmax</span><span class="o">=</span><span class="n">gzmax</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="n">params</span><span class="p">,</span><span class="n">Phi</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">g</span><span class="o">=</span><span class="n">g</span><span class="p">,</span><span class="n">gm</span><span class="o">=</span><span class="n">gm</span><span class="p">)</span>
</pre></div>
</div>
<p>The parameter <tt class="docutils literal"><span class="pre">select='0'</span></tt> tells that we want to compute the velocities on particles of type 0, while
<tt class="docutils literal"><span class="pre">disk=0</span></tt> tells that what we considere as the disk is only the particles 0. This is usefull when dealing with
multi-component models.</p>
<p>The latter function return <tt class="docutils literal"><span class="pre">nb</span></tt> with the new velocities, a matrix <tt class="docutils literal"><span class="pre">phi</span></tt> containing the potential at each node of the grid,
and a dictrionary <tt class="docutils literal"><span class="pre">stats</span></tt> containing some physical usefull quantities. Lets plot some of them:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="kn">import</span> <span class="nn">pylab</span> <span class="kn">as</span> <span class="nn">pt</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;vc&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;vm&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;sr&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;sz&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;sp&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s">&#39;Radius&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s">&#39;velocity&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<img alt="../_images/jean_cyl1.png" src="../_images/jean_cyl1.png" />
<p>Lets try instead of fixing the Tomre parameter, to use for the radial velocity dispertion
the epicyclic approximation.
This is done by using the following parameters:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_z</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;jeans&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="bp">None</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_r</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;epicyclic_approximation&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="mi">1</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">mode_sigma_p</span> <span class="o">=</span> <span class="p">{</span><span class="s">&quot;name&quot;</span><span class="p">:</span><span class="s">&quot;epicyclic_approximation&quot;</span><span class="p">,</span><span class="s">&quot;param&quot;</span><span class="p">:</span><span class="bp">None</span><span class="p">}</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">mode_sigma_z</span><span class="p">,</span><span class="n">mode_sigma_r</span><span class="p">,</span><span class="n">mode_sigma_p</span><span class="p">]</span>
</pre></div>
</div>
<p>again, run the magic function:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">nb</span><span class="p">,</span><span class="n">phi</span><span class="p">,</span><span class="n">stats</span> <span class="o">=</span> <span class="n">nb</span><span class="o">.</span><span class="n">Get_Velocities_From_Cylindrical_Grid</span><span class="p">(</span><span class="n">select</span><span class="o">=</span><span class="s">&#39;0&#39;</span><span class="p">,</span><span class="n">disk</span><span class="o">=</span><span class="p">(</span><span class="mi">0</span><span class="p">),</span><span class="n">eps</span><span class="o">=</span><span class="n">eps</span><span class="p">,</span><span class="n">nR</span><span class="o">=</span><span class="n">nr</span><span class="p">,</span><span class="n">nz</span><span class="o">=</span><span class="n">nz</span><span class="p">,</span><span class="n">nt</span><span class="o">=</span><span class="n">nt</span><span class="p">,</span><span class="n">Rmax</span><span class="o">=</span><span class="n">grmax</span><span class="p">,</span><span class="n">zmin</span><span class="o">=</span><span class="n">gzmin</span><span class="p">,</span><span class="n">zmax</span><span class="o">=</span><span class="n">gzmax</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="n">params</span><span class="p">,</span><span class="n">Phi</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">g</span><span class="o">=</span><span class="n">g</span><span class="p">,</span><span class="n">gm</span><span class="o">=</span><span class="n">gm</span><span class="p">)</span>
</pre></div>
</div>
<p>Now lets plot the Tommre parameter:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;R&#39;</span><span class="p">],</span><span class="n">stats</span><span class="p">[</span><span class="s">&#39;Q&#39;</span><span class="p">])</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s">&#39;Radius&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s">&#39;Q&#39;</span><span class="p">)</span>
<span class="gp">&gt;&gt;&gt; </span><span class="n">pt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
<img alt="../_images/jean_cyl2.png" src="../_images/jean_cyl2.png" />
<div class="section" id="multi-component-systems">
<h3>Multi component systems<a class="headerlink" href="#multi-component-systems" title="Permalink to this headline"></a></h3>
<p>Examples using multicomponents systems are provided in the <tt class="docutils literal"><span class="pre">pnbody_examples/ic</span></tt> directory obtained
with the command:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="n">pNbody_copy</span><span class="o">-</span><span class="n">examples</span>
</pre></div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<p class="logo"><a href="../index.html">
<img class="logo" src="../_static/icon-small.jpg" alt="Logo"/>
</a></p>
<h3><a href="../index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference internal" href="#">Generating velocities</a><ul>
<li><a class="reference internal" href="#spherical-coordinate">Spherical coordinate</a></li>
<li><a class="reference internal" href="#example-set-a-plummer-sphere-to-the-jeans-equilibrium">Example : Set a Plummer sphere to the jeans equilibrium</a></li>
<li><a class="reference internal" href="#cylindrical-coordinate">Cylindrical coordinate</a><ul>
<li><a class="reference internal" href="#vertical-velocity-dispersion">Vertical velocity dispersion</a></li>
<li><a class="reference internal" href="#radial-velocity-dispersion">Radial velocity dispersion</a></li>
<li><a class="reference internal" href="#azimuthal-velocity-dispersion">Azimuthal velocity dispersion</a></li>
<li><a class="reference internal" href="#mean-azimuthal-velocity">Mean azimuthal velocity</a></li>
</ul>
</li>
<li><a class="reference internal" href="#example-set-an-exponnential-disk-to-the-jeans-equilibrium">Example : Set an exponnential disk to the jeans equilibrium</a><ul>
<li><a class="reference internal" href="#multi-component-systems">Multi component systems</a></li>
</ul>
</li>
</ul>
</li>
</ul>
<h4>Previous topic</h4>
<p class="topless"><a href="GeneratingMassProfiles.html"
title="previous chapter">Generating mass profiles</a></p>
<h4>Next topic</h4>
<p class="topless"><a href="Selection.html"
title="next chapter">Selectig particles</a></p>
<h3>This Page</h3>
<ul class="this-page-menu">
<li><a href="../_sources/rst/GeneratingVelocities.txt"
rel="nofollow">Show Source</a></li>
</ul>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="../search.html" method="get">
<input type="text" name="q" size="18" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="../genindex.html" title="General Index"
>index</a></li>
<li class="right" >
<a href="../py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="../np-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="Selection.html" title="Selectig particles"
>next</a> |</li>
<li class="right" >
<a href="GeneratingMassProfiles.html" title="Generating mass profiles"
>previous</a> |</li>
<li><a href="../index.html">pNbody v4 documentation</a> &raquo;</li>
<li><a href="InitialConditions.html" >Generating initial conditions</a> &raquo;</li>
</ul>
</div>
<div class="footer">
&copy; Copyright 2011, Yves Revaz.
Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.0.7.
</div>
</body>
</html>

Event Timeline