Commits

John Wise committed ef7f9a8

Adding cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2.py

Comments (0)

Files changed (4)

html/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2-py.html

+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
+   "http://www.w3.org/TR/html4/strict.dtd">
+
+<html>
+<head>
+  <title></title>
+  <meta http-equiv="content-type" content="text/html; charset=latin1">
+  <style type="text/css">
+td.linenos { background-color: #f0f0f0; padding-right: 10px; }
+span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
+pre { line-height: 125%; }
+body .hll { background-color: #ffffcc }
+body  { background: #f8f8f8; }
+body .c { color: #408080; font-style: italic } /* Comment */
+body .err { border: 1px solid #FF0000 } /* Error */
+body .k { color: #008000; font-weight: bold } /* Keyword */
+body .o { color: #666666 } /* Operator */
+body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
+body .cp { color: #BC7A00 } /* Comment.Preproc */
+body .c1 { color: #408080; font-style: italic } /* Comment.Single */
+body .cs { color: #408080; font-style: italic } /* Comment.Special */
+body .gd { color: #A00000 } /* Generic.Deleted */
+body .ge { font-style: italic } /* Generic.Emph */
+body .gr { color: #FF0000 } /* Generic.Error */
+body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
+body .gi { color: #00A000 } /* Generic.Inserted */
+body .go { color: #808080 } /* Generic.Output */
+body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
+body .gs { font-weight: bold } /* Generic.Strong */
+body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
+body .gt { color: #0040D0 } /* Generic.Traceback */
+body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
+body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
+body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
+body .kp { color: #008000 } /* Keyword.Pseudo */
+body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
+body .kt { color: #B00040 } /* Keyword.Type */
+body .m { color: #666666 } /* Literal.Number */
+body .s { color: #BA2121 } /* Literal.String */
+body .na { color: #7D9029 } /* Name.Attribute */
+body .nb { color: #008000 } /* Name.Builtin */
+body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
+body .no { color: #880000 } /* Name.Constant */
+body .nd { color: #AA22FF } /* Name.Decorator */
+body .ni { color: #999999; font-weight: bold } /* Name.Entity */
+body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
+body .nf { color: #0000FF } /* Name.Function */
+body .nl { color: #A0A000 } /* Name.Label */
+body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
+body .nt { color: #008000; font-weight: bold } /* Name.Tag */
+body .nv { color: #19177C } /* Name.Variable */
+body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
+body .w { color: #bbbbbb } /* Text.Whitespace */
+body .mf { color: #666666 } /* Literal.Number.Float */
+body .mh { color: #666666 } /* Literal.Number.Hex */
+body .mi { color: #666666 } /* Literal.Number.Integer */
+body .mo { color: #666666 } /* Literal.Number.Oct */
+body .sb { color: #BA2121 } /* Literal.String.Backtick */
+body .sc { color: #BA2121 } /* Literal.String.Char */
+body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
+body .s2 { color: #BA2121 } /* Literal.String.Double */
+body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
+body .sh { color: #BA2121 } /* Literal.String.Heredoc */
+body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
+body .sx { color: #008000 } /* Literal.String.Other */
+body .sr { color: #BB6688 } /* Literal.String.Regex */
+body .s1 { color: #BA2121 } /* Literal.String.Single */
+body .ss { color: #19177C } /* Literal.String.Symbol */
+body .bp { color: #008000 } /* Name.Builtin.Pseudo */
+body .vc { color: #19177C } /* Name.Variable.Class */
+body .vg { color: #19177C } /* Name.Variable.Global */
+body .vi { color: #19177C } /* Name.Variable.Instance */
+body .il { color: #666666 } /* Literal.Number.Integer.Long */
+
+  </style>
+</head>
+<body>
+<h2></h2>
+
+<div class="highlight"><pre><span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">use</span><span class="p">;</span> <span class="n">use</span><span class="p">(</span><span class="s">&#39;Agg&#39;</span><span class="p">)</span>
+<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">rc</span>
+<span class="n">rc</span><span class="p">(</span><span class="s">&#39;font&#39;</span><span class="p">,</span> <span class="n">size</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
+<span class="kn">from</span> <span class="nn">matplotlib.ticker</span> <span class="kn">import</span> <span class="n">NullFormatter</span><span class="p">,</span> <span class="n">LinearLocator</span><span class="p">,</span> <span class="n">LogLocator</span>
+
+<span class="kn">import</span> <span class="nn">os</span><span class="o">,</span> <span class="nn">sys</span>
+<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">na</span>
+<span class="kn">import</span> <span class="nn">h5py</span> <span class="kn">as</span> <span class="nn">h5</span>
+<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>
+<span class="k">try</span><span class="p">:</span>
+    <span class="kn">from</span> <span class="nn">mpi4py</span> <span class="kn">import</span> <span class="n">MPI</span>
+    <span class="n">parallel</span> <span class="o">=</span> <span class="bp">True</span>
+<span class="k">except</span><span class="p">:</span>
+    <span class="n">parallel</span> <span class="o">=</span> <span class="bp">False</span>
+<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">log10</span>
+<span class="kn">from</span> <span class="nn">yt.mods</span> <span class="kn">import</span> <span class="o">*</span>
+<span class="kn">from</span> <span class="nn">myanyl</span> <span class="kn">import</span> <span class="o">*</span>
+<span class="kn">from</span> <span class="nn">yt.utilities.cosmology</span> <span class="kn">import</span> <span class="o">*</span>
+
+<span class="n">fn</span> <span class="o">=</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+<span class="n">output_file</span> <span class="o">=</span> <span class="s">&quot;anyl-files/&quot;</span> <span class="o">+</span> <span class="n">fn</span><span class="o">.</span><span class="n">split</span><span class="p">(</span><span class="s">&quot;/&quot;</span><span class="p">)[</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="s">&quot;_anyl.h5&quot;</span>
+<span class="c">#output_file = fn + &quot;_anyl.h5&quot;</span>
+<span class="n">pf</span> <span class="o">=</span> <span class="n">EnzoStaticOutput</span><span class="p">(</span><span class="n">fn</span><span class="p">,</span> <span class="n">data_style</span><span class="o">=</span><span class="s">&quot;enzo_packed_3d&quot;</span><span class="p">)</span>
+<span class="n">nhalos</span> <span class="o">=</span> <span class="mi">10</span>
+<span class="n">projwidth</span> <span class="o">=</span> <span class="mf">5.0</span>  <span class="c"># kpc</span>
+<span class="n">append_output</span> <span class="o">=</span> <span class="bp">True</span>
+<span class="n">zlim</span> <span class="o">=</span> <span class="p">{</span><span class="s">&#39;Density&#39;</span><span class="p">:</span> <span class="p">(</span><span class="mf">3e-27</span><span class="p">,</span> <span class="mf">1e-22</span><span class="p">)}</span>
+<span class="n">zticks</span> <span class="o">=</span> <span class="p">[</span><span class="mi">30</span><span class="p">,</span><span class="mi">25</span><span class="p">,</span><span class="mi">20</span><span class="p">,</span><span class="mi">15</span><span class="p">,</span><span class="mi">12</span><span class="p">,</span><span class="mi">10</span><span class="p">,</span><span class="mi">9</span><span class="p">,</span><span class="mi">8</span><span class="p">,</span><span class="mi">7</span><span class="p">,</span><span class="mi">6</span><span class="p">]</span>
+
+<span class="k">if</span> <span class="n">parallel</span><span class="p">:</span>
+    <span class="n">comm</span> <span class="o">=</span> <span class="n">MPI</span><span class="o">.</span><span class="n">COMM_WORLD</span>
+    <span class="n">rank</span> <span class="o">=</span> <span class="n">comm</span><span class="o">.</span><span class="n">Get_rank</span><span class="p">()</span>
+    <span class="n">nproc</span> <span class="o">=</span> <span class="n">comm</span><span class="o">.</span><span class="n">Get_size</span><span class="p">()</span>
+<span class="k">else</span><span class="p">:</span>
+    <span class="n">rank</span> <span class="o">=</span> <span class="mi">0</span>
+    <span class="n">nproc</span> <span class="o">=</span> <span class="mi">1</span>
+
+<span class="c">########################################################################</span>
+<span class="k">def</span> <span class="nf">filter_particle_type</span><span class="p">(</span><span class="n">source</span><span class="p">,</span> <span class="n">fields</span><span class="p">,</span> <span class="n">_type</span><span class="p">):</span>
+    <span class="c">#fields = ensure_list(fields)</span>
+    <span class="n">result</span> <span class="o">=</span> <span class="p">{}</span>
+    <span class="k">for</span> <span class="n">f</span> <span class="ow">in</span> <span class="n">fields</span><span class="p">:</span> <span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">array</span><span class="p">([])</span>
+    <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">grid</span> <span class="ow">in</span> <span class="nb">enumerate</span><span class="p">(</span><span class="n">source</span><span class="o">.</span><span class="n">_grids</span><span class="p">):</span>
+        <span class="k">if</span> <span class="n">i</span> <span class="o">%</span> <span class="n">nproc</span> <span class="o">==</span> <span class="n">rank</span><span class="p">:</span>
+            <span class="n">sel</span> <span class="o">=</span> <span class="n">grid</span><span class="p">[</span><span class="s">&#39;particle_type&#39;</span><span class="p">]</span> <span class="o">==</span> <span class="n">_type</span>
+            <span class="k">if</span> <span class="n">sel</span><span class="o">.</span><span class="n">any</span><span class="p">():</span>
+                <span class="k">for</span> <span class="n">f</span> <span class="ow">in</span> <span class="n">fields</span><span class="p">:</span>
+                    <span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">concatenate</span><span class="p">((</span><span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">],</span> <span class="n">grid</span><span class="p">[</span><span class="n">f</span><span class="p">][</span><span class="n">sel</span><span class="p">]))</span>
+            <span class="n">grid</span><span class="o">.</span><span class="n">clear_data</span><span class="p">()</span>
+    <span class="k">if</span> <span class="n">parallel</span><span class="p">:</span>
+        <span class="k">for</span> <span class="n">f</span> <span class="ow">in</span> <span class="n">fields</span><span class="p">:</span>
+            <span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">comm</span><span class="o">.</span><span class="n">gather</span><span class="p">(</span><span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">])</span>
+            <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+                <span class="n">temp</span> <span class="o">=</span> <span class="p">[]</span>
+                <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">]:</span>
+                    <span class="k">if</span> <span class="n">i</span> <span class="o">!=</span> <span class="p">[]</span> <span class="ow">and</span> <span class="n">i</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">:</span>
+                        <span class="n">temp</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">concatenate</span><span class="p">((</span><span class="n">temp</span><span class="p">,</span> <span class="n">i</span><span class="p">))</span>
+                <span class="n">result</span><span class="p">[</span><span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">temp</span>
+    <span class="k">return</span> <span class="n">result</span>
+<span class="c">########################################################################</span>
+<span class="k">def</span> <span class="nf">sphere_filter</span><span class="p">(</span><span class="n">particles</span><span class="p">,</span> <span class="n">sphere</span><span class="p">):</span>
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="k">return</span> <span class="bp">None</span>
+    <span class="n">rr</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">sphere</span><span class="o">.</span><span class="n">center</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">particles</span><span class="p">[</span><span class="s">&quot;particle_position_x&quot;</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span>
+                 <span class="p">(</span><span class="n">sphere</span><span class="o">.</span><span class="n">center</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">particles</span><span class="p">[</span><span class="s">&quot;particle_position_y&quot;</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span>
+                 <span class="p">(</span><span class="n">sphere</span><span class="o">.</span><span class="n">center</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="n">particles</span><span class="p">[</span><span class="s">&quot;particle_position_z&quot;</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
+    <span class="n">inside</span> <span class="o">=</span> <span class="n">rr</span> <span class="o">&lt;</span> <span class="n">sphere</span><span class="o">.</span><span class="n">radius</span>
+    <span class="k">if</span> <span class="n">inside</span><span class="o">.</span><span class="n">any</span><span class="p">():</span>
+        <span class="k">for</span> <span class="n">f</span> <span class="ow">in</span> <span class="n">particles</span><span class="o">.</span><span class="n">keys</span><span class="p">():</span>
+            <span class="n">particles</span><span class="p">[</span><span class="n">f</span><span class="p">]</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">f</span><span class="p">][</span><span class="n">inside</span><span class="p">]</span>
+    <span class="k">else</span><span class="p">:</span>
+        <span class="n">particles</span> <span class="o">=</span> <span class="bp">None</span>
+    <span class="k">return</span> <span class="n">particles</span>
+<span class="c">########################################################################</span>
+<span class="k">def</span> <span class="nf">one_halo</span><span class="p">(</span><span class="n">pf</span><span class="p">,</span> <span class="n">num</span><span class="p">,</span> <span class="n">center</span><span class="p">,</span> <span class="n">mass</span><span class="p">,</span> <span class="n">output</span><span class="p">):</span>
+
+    <span class="n">haloname</span> <span class="o">=</span> <span class="s">&quot;Halo</span><span class="si">%6.6d</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">num</span>
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">fp</span> <span class="o">=</span> <span class="n">h5</span><span class="o">.</span><span class="n">File</span><span class="p">(</span><span class="n">output</span><span class="p">)</span>
+        <span class="n">exists</span> <span class="o">=</span> <span class="n">haloname</span> <span class="ow">in</span> <span class="n">fp</span><span class="o">.</span><span class="n">listnames</span><span class="p">()</span>
+    <span class="k">else</span><span class="p">:</span>
+        <span class="n">exists</span> <span class="o">=</span> <span class="bp">None</span>
+    <span class="k">if</span> <span class="n">parallel</span><span class="p">:</span>
+        <span class="n">exists</span> <span class="o">=</span> <span class="n">comm</span><span class="o">.</span><span class="n">bcast</span><span class="p">([</span><span class="n">exists</span><span class="p">])</span>
+    <span class="k">if</span> <span class="n">exists</span><span class="p">[</span><span class="mi">0</span><span class="p">]:</span>
+        <span class="k">if</span> <span class="n">append_output</span><span class="p">:</span>
+            <span class="k">return</span>
+        <span class="k">else</span><span class="p">:</span>
+            <span class="k">del</span> <span class="n">pf</span><span class="p">[</span><span class="n">haloname</span><span class="p">]</span>
+
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">gid</span> <span class="o">=</span> <span class="n">fp</span><span class="o">.</span><span class="n">create_group</span><span class="p">(</span><span class="n">haloname</span><span class="p">)</span>
+    
+    <span class="c"># Calculate virial radius from mass, then find r200</span>
+    <span class="n">omegam</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">pf</span><span class="o">.</span><span class="n">get_parameter</span><span class="p">(</span><span class="s">&quot;CosmologyOmegaMatterNow&quot;</span><span class="p">))</span>
+    <span class="n">omegal</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">pf</span><span class="o">.</span><span class="n">get_parameter</span><span class="p">(</span><span class="s">&quot;CosmologyOmegaLambdaNow&quot;</span><span class="p">))</span>
+    <span class="n">zz</span> <span class="o">=</span> <span class="n">pf</span><span class="p">[</span><span class="s">&quot;CosmologyCurrentRedshift&quot;</span><span class="p">]</span>
+    <span class="n">hh</span> <span class="o">=</span> <span class="n">pf</span><span class="p">[</span><span class="s">&quot;CosmologyHubbleConstantNow&quot;</span><span class="p">]</span>
+    <span class="n">cosmo</span> <span class="o">=</span> <span class="n">Cosmology</span><span class="p">()</span>
+    <span class="n">ez</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">omegal</span> <span class="o">+</span> <span class="n">omegam</span> <span class="o">*</span> <span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">zz</span><span class="p">)</span><span class="o">**</span><span class="mi">3</span><span class="p">)</span>
+    <span class="n">omegamz</span> <span class="o">=</span> <span class="n">omegam</span> <span class="o">*</span> <span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">zz</span><span class="p">)</span><span class="o">**</span><span class="mi">3</span> <span class="o">/</span> <span class="n">ez</span><span class="o">**</span><span class="mi">2</span>
+    <span class="n">hz</span> <span class="o">=</span> <span class="p">(</span><span class="n">hh</span> <span class="o">*</span> <span class="mf">100.0</span> <span class="o">/</span> <span class="mf">3.086e19</span><span class="p">)</span> <span class="o">*</span> <span class="n">ez</span>
+    <span class="n">rvir</span> <span class="o">=</span> <span class="p">(</span><span class="mf">6.673e-8</span> <span class="o">*</span> <span class="n">mass</span> <span class="o">*</span> <span class="mf">1.989e33</span> <span class="o">/</span>
+            <span class="p">(</span><span class="mf">100.0</span> <span class="o">*</span> <span class="n">omegamz</span> <span class="o">*</span> <span class="n">hz</span><span class="o">**</span><span class="mi">2</span><span class="p">))</span><span class="o">**</span><span class="p">(</span><span class="mf">1.0</span><span class="o">/</span><span class="mi">3</span><span class="p">)</span>
+    <span class="n">sp</span> <span class="o">=</span> <span class="n">pf</span><span class="o">.</span><span class="n">h</span><span class="o">.</span><span class="n">sphere</span><span class="p">(</span><span class="n">center</span><span class="p">,</span> <span class="mf">1.5</span><span class="o">*</span><span class="n">rvir</span><span class="o">/</span><span class="n">pf</span><span class="p">[</span><span class="s">&quot;cm&quot;</span><span class="p">])</span>
+    <span class="n">prof</span> <span class="o">=</span> <span class="n">BinnedProfile1D</span><span class="p">(</span><span class="n">sp</span><span class="p">,</span> <span class="mi">100</span><span class="p">,</span> <span class="s">&quot;Radius&quot;</span><span class="p">,</span> <span class="mf">0.1</span><span class="o">*</span><span class="n">rvir</span><span class="p">,</span> <span class="mf">1.5</span><span class="o">*</span><span class="n">rvir</span><span class="p">)</span>
+    <span class="n">prof</span><span class="o">.</span><span class="n">add_fields</span><span class="p">([</span><span class="s">&quot;TotalMassMsun&quot;</span><span class="p">,</span> <span class="s">&quot;CellVolume&quot;</span><span class="p">])</span>
+    <span class="n">rhoc</span> <span class="o">=</span> <span class="n">cosmo</span><span class="o">.</span><span class="n">CriticalDensity</span><span class="p">(</span><span class="n">pf</span><span class="p">[</span><span class="s">&quot;CosmologyCurrentRedshift&quot;</span><span class="p">])</span>
+    <span class="n">avg_overden</span> <span class="o">=</span> <span class="mf">1.989e33</span><span class="o">*</span><span class="n">prof</span><span class="p">[</span><span class="s">&quot;TotalMassMsun&quot;</span><span class="p">]</span> <span class="o">/</span> <span class="n">prof</span><span class="p">[</span><span class="s">&quot;CellVolume&quot;</span><span class="p">]</span> <span class="o">/</span> \
+                  <span class="p">(</span><span class="n">rhoc</span><span class="o">*</span><span class="n">omegam</span><span class="p">)</span>
+    <span class="n">avg_overden</span> <span class="o">=</span> <span class="n">avg_overden</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="mi">0</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+    <span class="n">rr</span> <span class="o">=</span> <span class="n">prof</span><span class="p">[</span><span class="s">&#39;Radius&#39;</span><span class="p">][</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="mi">0</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+    <span class="n">r200</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">interp</span><span class="p">(</span><span class="mf">200.0</span><span class="p">,</span> <span class="n">avg_overden</span><span class="p">,</span> <span class="n">rr</span><span class="p">)</span>
+    <span class="k">del</span> <span class="n">prof</span>
+    <span class="k">del</span> <span class="n">sp</span>
+
+    <span class="c"># Redefine sphere with r200 instead of rvir</span>
+    <span class="n">sp</span> <span class="o">=</span> <span class="n">pf</span><span class="o">.</span><span class="n">h</span><span class="o">.</span><span class="n">sphere</span><span class="p">(</span><span class="n">center</span><span class="p">,</span> <span class="n">r200</span><span class="o">/</span><span class="n">pf</span><span class="p">[</span><span class="s">&quot;cm&quot;</span><span class="p">])</span>
+
+    <span class="c"># Global halo properties</span>
+    <span class="n">mgas</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;TotalQuantity&quot;</span><span class="p">](</span><span class="s">&quot;CellMassMsun&quot;</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
+    <span class="n">z2</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;WeightedAverageQuantity&quot;</span><span class="p">](</span>
+        <span class="s">&quot;Metallicity&quot;</span><span class="p">,</span> <span class="s">&quot;CellMass&quot;</span><span class="p">)</span>
+    <span class="n">z3</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;WeightedAverageQuantity&quot;</span><span class="p">](</span>
+        <span class="s">&quot;Metallicity3&quot;</span><span class="p">,</span> <span class="s">&quot;CellMass&quot;</span><span class="p">)</span>
+    <span class="n">mtot</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;TotalMass&quot;</span><span class="p">]()</span>
+    <span class="n">gas_cm</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;CenterOfMass&quot;</span><span class="p">]()</span>
+    <span class="n">gas_spin</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;BaryonSpinParameter&quot;</span><span class="p">]()</span>
+    <span class="n">dm_spin</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">quantities</span><span class="p">[</span><span class="s">&quot;ParticleSpinParameter&quot;</span><span class="p">]()</span>
+    <span class="n">gas_frac</span> <span class="o">=</span> <span class="n">mgas</span> <span class="o">/</span> <span class="n">mtot</span>
+    
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;GasMass&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">mgas</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;TotalMass&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">mtot</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;GasFraction&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gas_frac</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;Metallicity2&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">log10</span><span class="p">(</span><span class="n">z2</span><span class="p">)</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;Metallicity3&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">log10</span><span class="p">(</span><span class="n">z3</span><span class="p">)</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;GasCoM&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gas_cm</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;DarkMatterCoM&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">center</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;r200kpc&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">r200</span> <span class="o">/</span> <span class="mf">3.086e21</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;SpinGas&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gas_spin</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;SpinDM&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">dm_spin</span>
+
+    <span class="c"># Stellar properties</span>
+    <span class="n">star_type</span> <span class="o">=</span> <span class="mi">7</span>
+    <span class="n">star_fields</span> <span class="o">=</span> <span class="p">[</span><span class="s">&quot;ParticleMassMsun&quot;</span><span class="p">,</span> <span class="s">&quot;metallicity_fraction&quot;</span><span class="p">,</span>
+                   <span class="s">&quot;creation_time&quot;</span><span class="p">,</span> <span class="s">&quot;particle_index&quot;</span><span class="p">,</span> <span class="s">&quot;particle_position_x&quot;</span><span class="p">,</span>
+                   <span class="s">&quot;particle_position_y&quot;</span><span class="p">,</span> <span class="s">&quot;particle_position_z&quot;</span><span class="p">]</span>
+    <span class="n">Stars2</span> <span class="o">=</span> <span class="n">filter_particle_type</span><span class="p">(</span><span class="n">sp</span><span class="p">,</span> <span class="n">star_fields</span><span class="p">,</span> <span class="n">star_type</span><span class="p">)</span>
+    <span class="n">Stars2</span> <span class="o">=</span> <span class="n">sphere_filter</span><span class="p">(</span><span class="n">Stars2</span><span class="p">,</span> <span class="n">sp</span><span class="p">)</span>
+
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">nstars2</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">Stars2</span><span class="p">[</span><span class="n">star_fields</span><span class="p">[</span><span class="mi">0</span><span class="p">]])</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;NumberOfPop2Stars&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">nstars2</span>
+        <span class="k">print</span> <span class="s">&quot;</span><span class="si">%d</span><span class="s"> Pop II stars&quot;</span> <span class="o">%</span> <span class="n">nstars2</span>
+        <span class="k">if</span> <span class="n">nstars2</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+            <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">]</span> <span class="o">/=</span> <span class="mf">0.0204</span>
+            <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">])</span>
+            <span class="n">star_age</span> <span class="o">=</span> <span class="p">(</span><span class="n">pf</span><span class="p">[</span><span class="s">&quot;InitialTime&quot;</span><span class="p">]</span> <span class="o">-</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;creation_time&quot;</span><span class="p">])</span> <span class="o">*</span> \
+                       <span class="n">pf</span><span class="p">[</span><span class="s">&quot;years&quot;</span><span class="p">]</span> <span class="o">/</span> <span class="mf">1e6</span>
+            <span class="n">hz</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">histogram</span><span class="p">(</span><span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">],</span>
+                              <span class="n">weights</span><span class="o">=</span><span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;ParticleMassMsun&quot;</span><span class="p">],</span>
+                              <span class="nb">range</span> <span class="o">=</span> <span class="p">(</span><span class="o">-</span><span class="mi">6</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="n">bins</span><span class="o">=</span><span class="mi">48</span><span class="p">)</span>
+            <span class="n">ht</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">histogram</span><span class="p">(</span><span class="n">star_age</span><span class="p">,</span>
+                              <span class="n">weights</span><span class="o">=</span><span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;ParticleMassMsun&quot;</span><span class="p">],</span>
+                              <span class="n">bins</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_pdf&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">hz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_bin&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">hz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_age_pdf&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_age_bin&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_total_mass&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;ParticleMassMsun&quot;</span><span class="p">]</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_metallicity&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">star_age</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;ParticleMassMsun&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_position_x&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;particle_position_x&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_position_y&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;particle_position_y&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_position_z&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;particle_position_z&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop2_IDs&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars2</span><span class="p">[</span><span class="s">&quot;particle_index&quot;</span><span class="p">]</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="s">&quot;int&quot;</span><span class="p">)</span>
+
+    <span class="n">star_type</span> <span class="o">=</span> <span class="mi">5</span>
+    <span class="n">star_fields</span> <span class="o">=</span> <span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">,</span> <span class="s">&quot;creation_time&quot;</span><span class="p">,</span>
+                   <span class="s">&quot;particle_index&quot;</span><span class="p">,</span> <span class="s">&quot;particle_position_x&quot;</span><span class="p">,</span>
+                   <span class="s">&quot;particle_position_y&quot;</span><span class="p">,</span> <span class="s">&quot;particle_position_z&quot;</span><span class="p">]</span>
+    <span class="n">Stars3</span> <span class="o">=</span> <span class="n">filter_particle_type</span><span class="p">(</span><span class="n">sp</span><span class="p">,</span> <span class="n">star_fields</span><span class="p">,</span> <span class="n">star_type</span><span class="p">)</span>
+    <span class="n">Stars3</span> <span class="o">=</span> <span class="n">sphere_filter</span><span class="p">(</span><span class="n">Stars3</span><span class="p">,</span> <span class="n">sp</span><span class="p">)</span>
+
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">nstars3</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">Stars3</span><span class="p">[</span><span class="n">star_fields</span><span class="p">[</span><span class="mi">0</span><span class="p">]])</span>
+        <span class="n">gid</span><span class="p">[</span><span class="s">&quot;NumberOfPop3Stars&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">nstars3</span>
+        <span class="k">if</span> <span class="n">nstars3</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+            <span class="n">star_age</span> <span class="o">=</span> <span class="p">(</span><span class="n">pf</span><span class="p">[</span><span class="s">&quot;InitialTime&quot;</span><span class="p">]</span> <span class="o">-</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;creation_time&quot;</span><span class="p">])</span> <span class="o">*</span> \
+                       <span class="n">pf</span><span class="p">[</span><span class="s">&quot;years&quot;</span><span class="p">]</span> <span class="o">/</span> <span class="mf">1e6</span>
+            <span class="n">ht</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">histogram</span><span class="p">(</span><span class="n">star_age</span><span class="p">,</span> <span class="n">bins</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_age_pdf&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_age_bin&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_age&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">star_age</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_IDs&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;particle_index&quot;</span><span class="p">]</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="s">&#39;int&#39;</span><span class="p">)</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_metallicity&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;metallicity_fraction&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_position_x&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;particle_position_x&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_position_y&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;particle_position_y&quot;</span><span class="p">]</span>
+            <span class="n">gid</span><span class="p">[</span><span class="s">&quot;pop3_position_z&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Stars3</span><span class="p">[</span><span class="s">&quot;particle_position_z&quot;</span><span class="p">]</span>
+
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">fp</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
+        
+<span class="c">########################################################################</span>
+<span class="k">def</span> <span class="nf">plot_halo</span><span class="p">(</span><span class="n">_pf</span><span class="p">,</span> <span class="n">num</span><span class="p">,</span> <span class="n">output</span><span class="p">,</span> <span class="n">proj_fields</span><span class="o">=</span><span class="p">[</span><span class="s">&quot;Density&quot;</span><span class="p">],</span>
+              <span class="n">proj_dims</span><span class="o">=</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]):</span>
+    <span class="n">area_normalize</span> <span class="o">=</span> <span class="mf">1e2</span>  <span class="c"># Scatter area normalization (=1 pix in area)</span>
+    <span class="k">if</span> <span class="n">parallel</span><span class="p">:</span>
+        <span class="n">comm</span><span class="o">.</span><span class="n">barrier</span><span class="p">()</span>
+    <span class="n">fp</span> <span class="o">=</span> <span class="n">h5</span><span class="o">.</span><span class="n">File</span><span class="p">(</span><span class="n">output</span><span class="p">,</span> <span class="s">&quot;r&quot;</span><span class="p">)</span>
+    <span class="n">gid</span> <span class="o">=</span> <span class="n">fp</span><span class="p">[</span><span class="s">&quot;Halo</span><span class="si">%6.6d</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">num</span><span class="p">]</span>
+    <span class="n">data</span> <span class="o">=</span> <span class="p">{}</span>
+    <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="n">gid</span><span class="o">.</span><span class="n">listnames</span><span class="p">():</span>
+        <span class="n">data</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">gid</span><span class="p">[</span><span class="n">k</span><span class="p">]</span><span class="o">.</span><span class="n">value</span>
+    <span class="n">fp</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
+    <span class="c"># Convert stellar age to age of the universe</span>
+    <span class="n">tnow</span> <span class="o">=</span> <span class="mf">538.0</span> <span class="o">*</span> <span class="p">((</span><span class="n">_pf</span><span class="p">[</span><span class="s">&quot;CosmologyCurrentRedshift&quot;</span><span class="p">]</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="mf">10.0</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mf">1.5</span><span class="p">)</span>
+    <span class="k">if</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;NumberOfPop2Stars&quot;</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">tnow</span> <span class="o">-</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span>
+        <span class="n">massive</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mi">500</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span>
+        <span class="c"># Ignore SFR histogram because it uses 20 bins.  Instead use a</span>
+        <span class="c"># constant width.</span>
+        <span class="n">bin_width</span> <span class="o">=</span> <span class="mf">5.0</span>  <span class="c"># Myr</span>
+        <span class="n">nbins</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">())</span> <span class="o">/</span> <span class="n">bin_width</span><span class="p">)</span>
+        <span class="k">if</span> <span class="n">nbins</span> <span class="o">&lt;=</span> <span class="mi">1</span><span class="p">:</span>
+            <span class="n">nbins</span> <span class="o">=</span> <span class="p">[</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">],</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">bin_width</span><span class="p">,</span>
+                     <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="mi">2</span><span class="o">*</span><span class="n">bin_width</span><span class="p">]</span>
+        <span class="n">ht</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">histogram</span><span class="p">(</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">],</span> <span class="n">weights</span><span class="o">=</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">],</span>
+                              <span class="n">bins</span><span class="o">=</span><span class="n">nbins</span><span class="p">)</span>
+            
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_pdf&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_bin&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">ht</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
+    
+    <span class="c">#pc = PlotCollection(_pf, center=data[&#39;GasCoM&#39;])</span>
+    <span class="c">#colors = [&#39;b&#39;, &#39;c&#39;, &#39;w&#39;, &#39;y&#39;, &#39;g&#39;, &#39;r&#39;]  # star marker colors</span>
+    <span class="c">#ages = [20.0, 50.0, 100.0, 200.0, 500.0, 1000.0] # Myr</span>
+    <span class="c">#title1 = r&#39;(z = %0.2f) Halo %d&#39; % (_pf[&quot;CosmologyCurrentRedshift&quot;], num)</span>
+    <span class="c">#if data[&quot;NumberOfPop2Stars&quot;] &gt; 0:</span>
+    <span class="c">#    title2 = r&#39;log M = %0.2f, $f_{gas}$ = %0.3g, log M$_\star$ = %0.3g&#39; % \</span>
+    <span class="c">#             (log10(data[&quot;TotalMass&quot;]), data[&quot;GasFraction&quot;],</span>
+    <span class="c">#              log10(data[&quot;pop2_total_mass&quot;]))</span>
+    <span class="c">#else:</span>
+    <span class="c">#    title2 = r&#39;log M = %0.2f, $f_{gas}$ = %0.3g&#39; % \</span>
+    <span class="c">#             (log10(data[&quot;TotalMass&quot;]), data[&quot;GasFraction&quot;])</span>
+    <span class="c">#targ1 = {&#39;ha&#39;: &#39;left&#39;, &#39;color&#39;:&#39;w&#39;, &#39;va&#39;:&#39;top&#39;, &#39;size&#39;:&#39;medium&#39;,</span>
+    <span class="c">#         &#39;weight&#39;:&#39;bold&#39;}</span>
+    <span class="c">#targ2 = {&#39;ha&#39;: &#39;left&#39;, &#39;color&#39;:&#39;w&#39;, &#39;va&#39;:&#39;bottom&#39;, &#39;size&#39;:&#39;medium&#39;,</span>
+    <span class="c">#         &#39;weight&#39;:&#39;bold&#39;}</span>
+    <span class="c">#carg = {&#39;color&#39;:&#39;w&#39;}</span>
+    <span class="c">#for f in proj_fields:</span>
+    <span class="c">#    for dim in proj_dims:</span>
+    <span class="c">#        p=pc.add_projection(f, dim, &quot;Density&quot;)</span>
+    <span class="c">#        p.set_zlim(zlim[f][0], zlim[f][1])</span>
+    <span class="c">#        width = max(projwidth, 2*data[&quot;r200kpc&quot;])</span>
+    <span class="c">#        pc.set_width(width, &quot;kpc&quot;)</span>
+    <span class="c">#        boundL = data[&quot;GasCoM&quot;] - 0.5*width/pf[&quot;kpc&quot;]</span>
+    <span class="c">#        boundR = data[&quot;GasCoM&quot;] + 0.5*width/pf[&quot;kpc&quot;]</span>
+    <span class="c">#        p.modify[&quot;sphere&quot;](data[&quot;GasCoM&quot;], data[&quot;r200kpc&quot;]/pf[&quot;kpc&quot;],</span>
+    <span class="c">#                           circle_args=carg)</span>
+    <span class="c">#        tpos1 = [boundL[x_dict[p.data.axis]], boundR[y_dict[p.data.axis]]]</span>
+    <span class="c">#        tpos2 = [boundL[x_dict[p.data.axis]], boundL[y_dict[p.data.axis]]]</span>
+    <span class="c">#        p.modify[&quot;point&quot;](tpos1, title1, text_args=targ1)</span>
+    <span class="c">#        p.modify[&quot;point&quot;](tpos2, title2, text_args=targ2)</span>
+    <span class="c">#        for i in range(data[&quot;NumberOfPop2Stars&quot;]):</span>
+    <span class="c">#            pos = na.array(</span>
+    <span class="c">#                [data[&quot;pop2_position_x&quot;][i], data[&quot;pop2_position_y&quot;][i],</span>
+    <span class="c">#                 data[&quot;pop2_position_z&quot;][i]])</span>
+    <span class="c">#            age_bin = na.searchsorted(ages, data[&quot;pop2_age&quot;][i])</span>
+    <span class="c">#            ar = {&#39;color&#39;:colors[age_bin]}</span>
+    <span class="c">#            inside = na.logical_and(pos&gt;boundL, pos&lt;boundR).all()</span>
+    <span class="c">#            if inside:</span>
+    <span class="c">#                p.modify[&quot;marker&quot;](pos, marker=&quot;*&quot;, plot_args=ar)</span>
+    <span class="c">#pc.save(&quot;%s-halo%4.4d&quot; % (pf,num))</span>
+    <span class="c">#del pc</span>
+    <span class="k">if</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">and</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;NumberOfPop2Stars&quot;</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">nullfmt</span> <span class="o">=</span> <span class="n">NullFormatter</span><span class="p">()</span>
+        <span class="c"># Definitions for axes</span>
+        <span class="n">left</span><span class="p">,</span> <span class="n">width</span> <span class="o">=</span> <span class="mf">0.175</span><span class="p">,</span> <span class="mf">0.6</span>
+        <span class="n">bottom</span><span class="p">,</span> <span class="n">height</span> <span class="o">=</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mf">0.625</span>
+        <span class="n">left_h</span> <span class="o">=</span> <span class="n">left</span><span class="o">+</span><span class="n">width</span>
+        <span class="c">#bottom_h = left+width-0.025</span>
+        <span class="n">bottom_h</span> <span class="o">=</span> <span class="n">bottom</span><span class="o">+</span><span class="n">height</span>
+
+        <span class="n">rect_scatter</span> <span class="o">=</span> <span class="p">[</span><span class="n">left</span><span class="p">,</span> <span class="n">bottom</span><span class="p">,</span> <span class="n">width</span><span class="p">,</span> <span class="n">height</span><span class="p">]</span>
+        <span class="n">rect_histx</span> <span class="o">=</span> <span class="p">[</span><span class="n">left</span><span class="p">,</span> <span class="n">bottom_h</span><span class="p">,</span> <span class="n">width</span><span class="p">,</span> <span class="mf">0.175</span><span class="p">]</span>
+        <span class="n">rect_histy</span> <span class="o">=</span> <span class="p">[</span><span class="n">left_h</span><span class="p">,</span> <span class="n">bottom</span><span class="p">,</span> <span class="mf">0.2</span><span class="p">,</span> <span class="n">height</span><span class="p">]</span>
+
+        <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">figsize</span><span class="o">=</span><span class="p">(</span><span class="mi">6</span><span class="p">,</span><span class="mi">6</span><span class="p">))</span>
+
+        <span class="n">axScatter</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">axes</span><span class="p">(</span><span class="n">rect_scatter</span><span class="p">)</span>
+        <span class="n">axHistx</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">axes</span><span class="p">(</span><span class="n">rect_histx</span><span class="p">)</span>
+        <span class="n">axHisty</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">axes</span><span class="p">(</span><span class="n">rect_histy</span><span class="p">)</span>
+
+        <span class="c"># No labels</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">xaxis</span><span class="o">.</span><span class="n">set_major_formatter</span><span class="p">(</span><span class="n">nullfmt</span><span class="p">)</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">yaxis</span><span class="o">.</span><span class="n">set_major_formatter</span><span class="p">(</span><span class="n">nullfmt</span><span class="p">)</span>
+
+        <span class="c"># Scatter plot (only include massive star clusters)</span>
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">][</span><span class="n">massive</span><span class="p">]</span>
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity&quot;</span><span class="p">][</span><span class="n">massive</span><span class="p">]</span>
+        <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">]</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">][</span><span class="n">massive</span><span class="p">]</span>
+        <span class="n">isort</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">argsort</span><span class="p">(</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">])[::</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+        <span class="n">x</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age&quot;</span><span class="p">][</span><span class="n">isort</span><span class="p">]</span>
+        <span class="n">y</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity&quot;</span><span class="p">][</span><span class="n">isort</span><span class="p">]</span>
+        <span class="n">area</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_mass&quot;</span><span class="p">][</span><span class="n">isort</span><span class="p">]</span> <span class="o">/</span> <span class="n">area_normalize</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span>
+                          <span class="n">color</span><span class="o">=</span><span class="s">&#39;#888888&#39;</span><span class="p">,</span> <span class="n">s</span><span class="o">=</span><span class="n">area</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=</span><span class="mf">0.75</span><span class="p">,</span> <span class="n">edgecolors</span><span class="o">=</span><span class="s">&#39;k&#39;</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span><span class="o">-</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">)</span>
+        <span class="n">xl</span> <span class="o">=</span> <span class="n">axScatter</span><span class="o">.</span><span class="n">get_xlim</span><span class="p">()</span>
+        <span class="n">yl</span> <span class="o">=</span> <span class="n">axScatter</span><span class="o">.</span><span class="n">get_ylim</span><span class="p">()</span>
+        <span class="n">xl</span> <span class="o">=</span> <span class="p">(</span><span class="mf">150.0</span><span class="p">,</span> <span class="mf">1.01</span><span class="o">*</span><span class="n">tnow</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">xl</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">xl</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">r&#39;Time [Myr]&#39;</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">r&#39;[Z/H]&#39;</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">yaxis</span><span class="o">.</span><span class="n">set_label_coords</span><span class="p">(</span><span class="o">-</span><span class="mf">0.19</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="bp">True</span><span class="p">)</span>
+        <span class="n">margin</span> <span class="o">=</span> <span class="p">[</span><span class="mf">0.02</span> <span class="o">*</span> <span class="p">(</span><span class="n">xl</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">xl</span><span class="p">[</span><span class="mi">0</span><span class="p">]),</span> <span class="mf">0.02</span> <span class="o">*</span> <span class="p">(</span><span class="n">yl</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">yl</span><span class="p">[</span><span class="mi">0</span><span class="p">])]</span>
+        <span class="c">#txt = &#39;Halo %d&#39; % num</span>
+        <span class="c">#axScatter.text(xl[1]-margin[0], yl[0]+margin[1], txt, ha=&quot;right&quot;,</span>
+        <span class="c">#               size=&#39;small&#39;)</span>
+
+        <span class="c"># Draw log-linear &quot;fits&quot; to merged stellar populations for halo 0</span>
+        <span class="k">if</span> <span class="n">num</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+            <span class="n">axScatter</span><span class="o">.</span><span class="n">plot</span><span class="p">([</span><span class="mi">450</span><span class="p">,</span><span class="mi">750</span><span class="p">],</span> <span class="p">[</span><span class="o">-</span><span class="mf">3.5</span><span class="p">,</span> <span class="o">-</span><span class="mf">2.5</span><span class="p">],</span> <span class="s">&quot;k--&quot;</span><span class="p">)</span>
+            <span class="n">axScatter</span><span class="o">.</span><span class="n">plot</span><span class="p">([</span><span class="mi">450</span><span class="p">,</span><span class="mi">750</span><span class="p">],</span> <span class="p">[</span><span class="o">-</span><span class="mf">3.2</span><span class="p">,</span> <span class="o">-</span><span class="mf">2.1</span><span class="p">],</span> <span class="s">&quot;k--&quot;</span><span class="p">)</span>
+
+        <span class="c"># Plot sample masses</span>
+        <span class="n">props</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span><span class="n">edgecolors</span><span class="o">=</span><span class="s">&#39;k&#39;</span><span class="p">)</span>
+        <span class="n">masses</span> <span class="o">=</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span>
+        <span class="n">rbig</span> <span class="o">=</span> <span class="p">(</span><span class="mf">10.0</span><span class="o">**</span><span class="n">masses</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">/</span><span class="n">area_normalize</span><span class="o">/</span><span class="n">na</span><span class="o">.</span><span class="n">pi</span><span class="p">)</span><span class="o">**</span><span class="mf">0.5</span>
+        <span class="k">for</span> <span class="n">m</span> <span class="ow">in</span> <span class="n">masses</span><span class="p">:</span>
+            <span class="n">label</span> <span class="o">=</span> <span class="s">r&#39;$10^{</span><span class="si">%d</span><span class="s">}$&#39;</span> <span class="o">%</span> <span class="n">m</span>
+            <span class="n">rthis</span> <span class="o">=</span> <span class="p">(</span><span class="mf">10.0</span><span class="o">**</span><span class="n">m</span><span class="o">/</span><span class="n">area_normalize</span><span class="o">/</span><span class="n">na</span><span class="o">.</span><span class="n">pi</span><span class="p">)</span><span class="o">**</span><span class="mf">0.5</span>
+            <span class="n">axScatter</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">xl</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">+</span><span class="mi">3</span><span class="o">*</span><span class="n">margin</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span>
+                              <span class="n">yl</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mi">3</span><span class="o">*</span><span class="n">margin</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="p">(</span><span class="n">rbig</span><span class="o">-</span><span class="n">rthis</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">yl</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">yl</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span><span class="o">/</span><span class="mi">400</span><span class="p">,</span>
+                              <span class="n">c</span><span class="o">=</span><span class="s">&quot;#FFFFFF&quot;</span><span class="p">,</span>
+                              <span class="n">s</span><span class="o">=</span><span class="mf">10.0</span><span class="o">**</span><span class="n">m</span><span class="o">/</span><span class="n">area_normalize</span><span class="p">,</span>
+                              <span class="n">label</span><span class="o">=</span><span class="n">label</span><span class="p">,</span> <span class="o">**</span><span class="n">props</span><span class="p">)</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">xl</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">xl</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
+        <span class="n">axScatter</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span><span class="n">yl</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">yl</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
+
+        <span class="c"># Histograms</span>
+        <span class="n">nmajor</span> <span class="o">=</span> <span class="mi">5</span>
+        <span class="n">width</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_bin&quot;</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_bin&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span>
+        <span class="nb">bin</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_bin&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+        <span class="n">height</span> <span class="o">=</span> <span class="mf">1e-6</span><span class="o">*</span><span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_age_pdf&quot;</span><span class="p">]</span> <span class="o">/</span> <span class="n">width</span>
+        <span class="n">nearest_power10</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">na</span><span class="o">.</span><span class="n">log10</span><span class="p">(</span><span class="n">height</span><span class="o">.</span><span class="n">max</span><span class="p">()))</span><span class="o">-</span><span class="mi">1</span>
+        <span class="n">height</span> <span class="o">/=</span> <span class="mf">10.0</span><span class="o">**</span><span class="n">nearest_power10</span>
+<span class="c">#        bin = 0.5 * (data[&quot;pop2_age_bin&quot;][0:-1] + data[&quot;pop2_age_bin&quot;][1:])</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">bar</span><span class="p">(</span><span class="nb">bin</span><span class="p">,</span> <span class="n">height</span><span class="p">,</span> <span class="n">width</span><span class="p">,</span> <span class="n">color</span><span class="o">=</span><span class="s">&quot;#AAAAAA&quot;</span><span class="p">)</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">r&#39;SFR [$10^{</span><span class="si">%d</span><span class="s">} M_\odot$/yr]&#39;</span> <span class="o">%</span> <span class="n">nearest_power10</span><span class="p">)</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">yaxis</span><span class="o">.</span><span class="n">set_label_coords</span><span class="p">(</span><span class="o">-</span><span class="mf">0.19</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">)</span>
+        <span class="n">x0</span><span class="p">,</span><span class="n">x1</span> <span class="o">=</span> <span class="n">axScatter</span><span class="o">.</span><span class="n">get_xlim</span><span class="p">()</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span> <span class="n">x0</span><span class="p">,</span><span class="n">x1</span> <span class="p">)</span>
+        
+        <span class="c">#axHistx.yaxis.set_major_locator(LinearLocator(numticks=nmajor))</span>
+        <span class="c">#ticks = axHistx.yaxis.get_major_ticks()</span>
+        <span class="c">#ticks[0].label1On = False</span>
+
+        <span class="n">width</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_bin&quot;</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_bin&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span>
+        <span class="nb">bin</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_bin&quot;</span><span class="p">][</span><span class="mi">0</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+<span class="c">#        bin = 0.5 * (data[&quot;pop2_metallicity_bin&quot;][0:-1] +</span>
+<span class="c">#                     data[&quot;pop2_metallicity_bin&quot;][1:])</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">barh</span><span class="p">(</span><span class="nb">bin</span><span class="p">,</span> <span class="n">data</span><span class="p">[</span><span class="s">&quot;pop2_metallicity_pdf&quot;</span><span class="p">],</span> <span class="n">height</span><span class="o">=</span><span class="n">width</span><span class="p">,</span>
+                     <span class="n">color</span><span class="o">=</span><span class="s">&quot;#AAAAAA&quot;</span><span class="p">)</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">r&#39;M$_\star$ [$M_\odot$]&#39;</span><span class="p">)</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span> <span class="n">axScatter</span><span class="o">.</span><span class="n">get_ylim</span><span class="p">()</span> <span class="p">)</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">xaxis</span><span class="o">.</span><span class="n">set_label_position</span><span class="p">(</span><span class="s">&#39;top&#39;</span><span class="p">)</span>
+        <span class="n">axHisty</span><span class="o">.</span><span class="n">xaxis</span><span class="o">.</span><span class="n">set_major_locator</span><span class="p">(</span><span class="n">LinearLocator</span><span class="p">(</span><span class="n">numticks</span><span class="o">=</span><span class="n">nmajor</span><span class="p">))</span>
+        <span class="n">ticks</span> <span class="o">=</span> <span class="n">axHisty</span><span class="o">.</span><span class="n">xaxis</span><span class="o">.</span><span class="n">get_major_ticks</span><span class="p">()</span>
+        <span class="k">for</span> <span class="n">t</span> <span class="ow">in</span> <span class="n">ticks</span><span class="p">:</span>
+            <span class="n">t</span><span class="o">.</span><span class="n">label1On</span> <span class="o">=</span> <span class="bp">False</span>
+            <span class="n">t</span><span class="o">.</span><span class="n">label2On</span> <span class="o">=</span> <span class="bp">True</span>
+            <span class="n">t</span><span class="o">.</span><span class="n">label2</span><span class="o">.</span><span class="n">set_rotation</span><span class="p">(</span><span class="mi">270</span><span class="p">)</span>
+        <span class="n">ticks</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">label2On</span> <span class="o">=</span> <span class="bp">False</span>        
+
+        <span class="n">x0</span><span class="p">,</span><span class="n">x1</span> <span class="o">=</span> <span class="n">axHistx</span><span class="o">.</span><span class="n">get_xlim</span><span class="p">()</span>
+        <span class="n">y0</span><span class="p">,</span><span class="n">y1</span> <span class="o">=</span> <span class="n">axHistx</span><span class="o">.</span><span class="n">get_ylim</span><span class="p">()</span>
+        <span class="n">p1</span> <span class="o">=</span> <span class="n">axHistx</span><span class="o">.</span><span class="n">twiny</span><span class="p">()</span>
+        <span class="n">p1</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;Redshift&quot;</span><span class="p">)</span>
+        <span class="n">tnow</span> <span class="o">=</span> <span class="mf">538.0</span> <span class="o">*</span> <span class="p">((</span><span class="n">_pf</span><span class="p">[</span><span class="s">&quot;CosmologyCurrentRedshift&quot;</span><span class="p">]</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="mf">10.0</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mf">1.5</span><span class="p">)</span>
+        <span class="n">tticks</span> <span class="o">=</span> <span class="p">[]</span>
+        <span class="k">for</span> <span class="n">tick</span> <span class="ow">in</span> <span class="n">zticks</span><span class="p">:</span>
+            <span class="n">thisz</span> <span class="o">=</span> <span class="mf">538.0</span> <span class="o">*</span> <span class="p">((</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">tick</span><span class="p">)</span><span class="o">/</span><span class="mi">10</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mf">1.5</span><span class="p">)</span>
+            <span class="n">tticks</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">thisz</span><span class="p">)</span>
+        <span class="c">#p1.set_xscale(&#39;log&#39;)</span>
+        <span class="n">p1</span><span class="o">.</span><span class="n">set_xticks</span><span class="p">(</span><span class="n">tticks</span><span class="p">)</span>
+        <span class="n">p1</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">x0</span><span class="p">,</span><span class="n">x1</span><span class="p">)</span>
+        <span class="n">p1</span><span class="o">.</span><span class="n">set_xticklabels</span><span class="p">([</span><span class="s">&quot;</span><span class="si">%d</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">e</span> <span class="k">for</span> <span class="n">e</span> <span class="ow">in</span> <span class="n">zticks</span><span class="p">])</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">x0</span><span class="p">,</span><span class="n">x1</span><span class="p">)</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span><span class="n">y0</span><span class="p">,</span><span class="n">y1</span><span class="p">)</span>
+        <span class="n">axHistx</span><span class="o">.</span><span class="n">yaxis</span><span class="o">.</span><span class="n">set_major_locator</span><span class="p">(</span><span class="n">LinearLocator</span><span class="p">(</span><span class="n">numticks</span><span class="o">=</span><span class="n">nmajor</span><span class="p">))</span>
+        <span class="n">ticks</span> <span class="o">=</span> <span class="n">axHistx</span><span class="o">.</span><span class="n">yaxis</span><span class="o">.</span><span class="n">get_major_ticks</span><span class="p">()</span>
+        <span class="n">ticks</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">label1On</span> <span class="o">=</span> <span class="bp">False</span>
+        <span class="n">name</span> <span class="o">=</span> <span class="s">&quot;Halo </span><span class="si">%d</span><span class="s">&quot;</span> <span class="o">%</span> <span class="n">num</span>
+        <span class="k">if</span> <span class="n">name</span> <span class="o">!=</span> <span class="bp">None</span><span class="p">:</span>
+            <span class="n">axHistx</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="n">x0</span><span class="o">+</span><span class="mf">0.02</span><span class="o">*</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x0</span><span class="p">),</span><span class="n">y1</span><span class="o">-</span><span class="mf">0.05</span><span class="o">*</span><span class="p">(</span><span class="n">y1</span><span class="o">-</span><span class="n">y0</span><span class="p">),</span><span class="n">name</span><span class="p">,</span>
+                         <span class="n">ha</span><span class="o">=</span><span class="s">&quot;left&quot;</span><span class="p">,</span><span class="n">va</span><span class="o">=</span><span class="s">&quot;top&quot;</span><span class="p">)</span>
+
+        <span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s">&quot;</span><span class="si">%s</span><span class="s">-halo</span><span class="si">%4.4d</span><span class="s">-stars.eps&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">_pf</span><span class="p">,</span><span class="n">num</span><span class="p">))</span>
+        <span class="n">plt</span><span class="o">.</span><span class="n">clf</span><span class="p">()</span>
+
+    <span class="k">return</span>
+    
+<span class="c">########################################################################</span>
+
+<span class="k">if</span> <span class="n">append_output</span> <span class="ow">is</span> <span class="bp">False</span> <span class="ow">and</span> <span class="n">rank</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
+    <span class="k">if</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">output_file</span><span class="p">):</span>
+        <span class="n">os</span><span class="o">.</span><span class="n">remove</span><span class="p">(</span><span class="n">output_file</span><span class="p">)</span>
+
+<span class="c"># Make projection of fields</span>
+<span class="n">fields</span> <span class="o">=</span> <span class="p">[</span><span class="s">&quot;Density&quot;</span><span class="p">]</span>
+<span class="c">#pc = PlotCollection(pf, center=[0.5]*3)</span>
+<span class="c">#for f in fields:</span>
+<span class="c">#    for dim in range(3):</span>
+<span class="c">#        pc.add_projection(f, dim, &quot;Density&quot;)</span>
+<span class="c">#del pc</span>
+
+<span class="c"># Load halos</span>
+<span class="n">cycle</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">pf</span><span class="o">.</span><span class="n">get_parameter</span><span class="p">(</span><span class="s">&quot;InitialCycleNumber&quot;</span><span class="p">))</span>
+<span class="n">halos</span> <span class="o">=</span> <span class="n">readhalos</span><span class="p">(</span><span class="n">cycle</span><span class="o">=</span><span class="n">cycle</span><span class="p">,</span> <span class="n">maxhalos</span><span class="o">=</span><span class="n">nhalos</span><span class="p">)</span>
+<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">halos</span><span class="p">[</span><span class="s">&quot;mass&quot;</span><span class="p">])):</span>
+    <span class="n">one_halo</span><span class="p">(</span><span class="n">pf</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">halos</span><span class="p">[</span><span class="s">&quot;pos&quot;</span><span class="p">][</span><span class="n">i</span><span class="p">],</span> <span class="n">halos</span><span class="p">[</span><span class="s">&quot;mass&quot;</span><span class="p">][</span><span class="n">i</span><span class="p">],</span> <span class="n">output_file</span><span class="p">)</span>
+    <span class="n">plot_halo</span><span class="p">(</span><span class="n">pf</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">output_file</span><span class="p">)</span>
+
+<span class="n">MPI</span><span class="o">.</span><span class="n">Finalize</span><span class="p">()</span>
+</pre></div>
+</body>
+</html>
-var inventory_data = [];
+var inventory_data = [
+ {
+  "username": "John Wise <jwise@astro.princeton.edu>", 
+  "modtime": 1300885454, 
+  "name": "analysis2.py", 
+  "descr": "Test paste.  Dwarf galaxy analysis script.", 
+  "htmlname": "html/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2-py.html", 
+  "modified": "Wed Mar 23 09:04:14 2011", 
+  "lastmod_hash": "tip", 
+  "fullname": "posts/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2.py"
+ }
+];

posts/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2.py

+from matplotlib import use; use('Agg')
+from matplotlib import rc
+rc('font', size=20)
+from matplotlib.ticker import NullFormatter, LinearLocator, LogLocator
+
+import os, sys
+import numpy as na
+import h5py as h5
+import matplotlib.pyplot as plt
+try:
+    from mpi4py import MPI
+    parallel = True
+except:
+    parallel = False
+from math import log10
+from yt.mods import *
+from myanyl import *
+from yt.utilities.cosmology import *
+
+fn = sys.argv[-1]
+output_file = "anyl-files/" + fn.split("/")[1] + "_anyl.h5"
+#output_file = fn + "_anyl.h5"
+pf = EnzoStaticOutput(fn, data_style="enzo_packed_3d")
+nhalos = 10
+projwidth = 5.0  # kpc
+append_output = True
+zlim = {'Density': (3e-27, 1e-22)}
+zticks = [30,25,20,15,12,10,9,8,7,6]
+
+if parallel:
+    comm = MPI.COMM_WORLD
+    rank = comm.Get_rank()
+    nproc = comm.Get_size()
+else:
+    rank = 0
+    nproc = 1
+
+########################################################################
+def filter_particle_type(source, fields, _type):
+    #fields = ensure_list(fields)
+    result = {}
+    for f in fields: result[f] = na.array([])
+    for i,grid in enumerate(source._grids):
+        if i % nproc == rank:
+            sel = grid['particle_type'] == _type
+            if sel.any():
+                for f in fields:
+                    result[f] = na.concatenate((result[f], grid[f][sel]))
+            grid.clear_data()
+    if parallel:
+        for f in fields:
+            result[f] = comm.gather(result[f])
+            if rank == 0:
+                temp = []
+                for i in result[f]:
+                    if i != [] and i != None:
+                        temp = na.concatenate((temp, i))
+                result[f] = temp
+    return result
+########################################################################
+def sphere_filter(particles, sphere):
+    if rank > 0:
+        return None
+    rr = na.sqrt((sphere.center[0] - particles["particle_position_x"])**2 +
+                 (sphere.center[1] - particles["particle_position_y"])**2 +
+                 (sphere.center[2] - particles["particle_position_z"])**2)
+    inside = rr < sphere.radius
+    if inside.any():
+        for f in particles.keys():
+            particles[f] = particles[f][inside]
+    else:
+        particles = None
+    return particles
+########################################################################
+def one_halo(pf, num, center, mass, output):
+
+    haloname = "Halo%6.6d" % num
+    if rank == 0:
+        fp = h5.File(output)
+        exists = haloname in fp.listnames()
+    else:
+        exists = None
+    if parallel:
+        exists = comm.bcast([exists])
+    if exists[0]:
+        if append_output:
+            return
+        else:
+            del pf[haloname]
+
+    if rank == 0:
+        gid = fp.create_group(haloname)
+    
+    # Calculate virial radius from mass, then find r200
+    omegam = float(pf.get_parameter("CosmologyOmegaMatterNow"))
+    omegal = float(pf.get_parameter("CosmologyOmegaLambdaNow"))
+    zz = pf["CosmologyCurrentRedshift"]
+    hh = pf["CosmologyHubbleConstantNow"]
+    cosmo = Cosmology()
+    ez = na.sqrt(omegal + omegam * (1+zz)**3)
+    omegamz = omegam * (1+zz)**3 / ez**2
+    hz = (hh * 100.0 / 3.086e19) * ez
+    rvir = (6.673e-8 * mass * 1.989e33 /
+            (100.0 * omegamz * hz**2))**(1.0/3)
+    sp = pf.h.sphere(center, 1.5*rvir/pf["cm"])
+    prof = BinnedProfile1D(sp, 100, "Radius", 0.1*rvir, 1.5*rvir)
+    prof.add_fields(["TotalMassMsun", "CellVolume"])
+    rhoc = cosmo.CriticalDensity(pf["CosmologyCurrentRedshift"])
+    avg_overden = 1.989e33*prof["TotalMassMsun"] / prof["CellVolume"] / \
+                  (rhoc*omegam)
+    avg_overden = avg_overden[-1:0:-1]
+    rr = prof['Radius'][-1:0:-1]
+    r200 = na.interp(200.0, avg_overden, rr)
+    del prof
+    del sp
+
+    # Redefine sphere with r200 instead of rvir
+    sp = pf.h.sphere(center, r200/pf["cm"])
+
+    # Global halo properties
+    mgas = sp.quantities["TotalQuantity"]("CellMassMsun")[0]
+    z2 = sp.quantities["WeightedAverageQuantity"](
+        "Metallicity", "CellMass")
+    z3 = sp.quantities["WeightedAverageQuantity"](
+        "Metallicity3", "CellMass")
+    mtot = sp.quantities["TotalMass"]()
+    gas_cm = sp.quantities["CenterOfMass"]()
+    gas_spin = sp.quantities["BaryonSpinParameter"]()
+    dm_spin = sp.quantities["ParticleSpinParameter"]()
+    gas_frac = mgas / mtot
+    
+    if rank == 0:
+        gid["GasMass"] = mgas
+        gid["TotalMass"] = mtot
+        gid["GasFraction"] = gas_frac
+        gid["Metallicity2"] = log10(z2)
+        gid["Metallicity3"] = log10(z3)
+        gid["GasCoM"] = gas_cm
+        gid["DarkMatterCoM"] = center
+        gid["r200kpc"] = r200 / 3.086e21
+        gid["SpinGas"] = gas_spin
+        gid["SpinDM"] = dm_spin
+
+    # Stellar properties
+    star_type = 7
+    star_fields = ["ParticleMassMsun", "metallicity_fraction",
+                   "creation_time", "particle_index", "particle_position_x",
+                   "particle_position_y", "particle_position_z"]
+    Stars2 = filter_particle_type(sp, star_fields, star_type)
+    Stars2 = sphere_filter(Stars2, sp)
+
+    if rank == 0:
+        nstars2 = len(Stars2[star_fields[0]])
+        gid["NumberOfPop2Stars"] = nstars2
+        print "%d Pop II stars" % nstars2
+        if nstars2 > 0:
+            Stars2["metallicity_fraction"] /= 0.0204
+            Stars2["metallicity_fraction"] = na.log10(Stars2["metallicity_fraction"])
+            star_age = (pf["InitialTime"] - Stars2["creation_time"]) * \
+                       pf["years"] / 1e6
+            hz = na.histogram(Stars2["metallicity_fraction"],
+                              weights=Stars2["ParticleMassMsun"],
+                              range = (-6, 0), bins=48)
+            ht = na.histogram(star_age,
+                              weights=Stars2["ParticleMassMsun"],
+                              bins=20)
+            gid["pop2_metallicity_pdf"] = hz[0]
+            gid["pop2_metallicity_bin"] = hz[1]
+            gid["pop2_age_pdf"] = ht[0]
+            gid["pop2_age_bin"] = ht[1]
+            gid["pop2_total_mass"] = Stars2["ParticleMassMsun"].sum()
+            gid["pop2_metallicity"] = Stars2["metallicity_fraction"]
+            gid["pop2_age"] = star_age
+            gid["pop2_mass"] = Stars2["ParticleMassMsun"]
+            gid["pop2_position_x"] = Stars2["particle_position_x"]
+            gid["pop2_position_y"] = Stars2["particle_position_y"]
+            gid["pop2_position_z"] = Stars2["particle_position_z"]
+            gid["pop2_IDs"] = Stars2["particle_index"].astype("int")
+
+    star_type = 5
+    star_fields = ["metallicity_fraction", "creation_time",
+                   "particle_index", "particle_position_x",
+                   "particle_position_y", "particle_position_z"]
+    Stars3 = filter_particle_type(sp, star_fields, star_type)
+    Stars3 = sphere_filter(Stars3, sp)
+
+    if rank == 0:
+        nstars3 = len(Stars3[star_fields[0]])
+        gid["NumberOfPop3Stars"] = nstars3
+        if nstars3 > 0:
+            star_age = (pf["InitialTime"] - Stars3["creation_time"]) * \
+                       pf["years"] / 1e6
+            ht = na.histogram(star_age, bins=20)
+            gid["pop3_age_pdf"] = ht[0]
+            gid["pop3_age_bin"] = ht[1]
+            gid["pop3_age"] = star_age
+            gid["pop3_IDs"] = Stars3["particle_index"].astype('int')
+            gid["pop3_metallicity"] = Stars3["metallicity_fraction"]
+            gid["pop3_position_x"] = Stars3["particle_position_x"]
+            gid["pop3_position_y"] = Stars3["particle_position_y"]
+            gid["pop3_position_z"] = Stars3["particle_position_z"]
+
+    if rank == 0:
+        fp.close()
+        
+########################################################################
+def plot_halo(_pf, num, output, proj_fields=["Density"],
+              proj_dims=[0,1,2]):
+    area_normalize = 1e2  # Scatter area normalization (=1 pix in area)
+    if parallel:
+        comm.barrier()
+    fp = h5.File(output, "r")
+    gid = fp["Halo%6.6d" % num]
+    data = {}
+    for k in gid.listnames():
+        data[k] = gid[k].value
+    fp.close()
+    # Convert stellar age to age of the universe
+    tnow = 538.0 * ((_pf["CosmologyCurrentRedshift"]+1)/10.0)**(-1.5)
+    if data["NumberOfPop2Stars"] > 0:
+        data["pop2_age"] = tnow - data["pop2_age"]
+        massive = na.where(data["pop2_mass"] > 500)[0]
+        # Ignore SFR histogram because it uses 20 bins.  Instead use a
+        # constant width.
+        bin_width = 5.0  # Myr
+        nbins = int((data["pop2_age"].max() - data["pop2_age"].min()) / bin_width)
+        if nbins <= 1:
+            nbins = [data["pop2_age"][0], data["pop2_age"][0] + bin_width,
+                     data["pop2_age"][0] + 2*bin_width]
+        ht = na.histogram(data["pop2_age"], weights=data["pop2_mass"],
+                              bins=nbins)
+            
+        data["pop2_age_pdf"] = ht[0]
+        data["pop2_age_bin"] = ht[1]
+    
+    #pc = PlotCollection(_pf, center=data['GasCoM'])
+    #colors = ['b', 'c', 'w', 'y', 'g', 'r']  # star marker colors
+    #ages = [20.0, 50.0, 100.0, 200.0, 500.0, 1000.0] # Myr
+    #title1 = r'(z = %0.2f) Halo %d' % (_pf["CosmologyCurrentRedshift"], num)
+    #if data["NumberOfPop2Stars"] > 0:
+    #    title2 = r'log M = %0.2f, $f_{gas}$ = %0.3g, log M$_\star$ = %0.3g' % \
+    #             (log10(data["TotalMass"]), data["GasFraction"],
+    #              log10(data["pop2_total_mass"]))
+    #else:
+    #    title2 = r'log M = %0.2f, $f_{gas}$ = %0.3g' % \
+    #             (log10(data["TotalMass"]), data["GasFraction"])
+    #targ1 = {'ha': 'left', 'color':'w', 'va':'top', 'size':'medium',
+    #         'weight':'bold'}
+    #targ2 = {'ha': 'left', 'color':'w', 'va':'bottom', 'size':'medium',
+    #         'weight':'bold'}
+    #carg = {'color':'w'}
+    #for f in proj_fields:
+    #    for dim in proj_dims:
+    #        p=pc.add_projection(f, dim, "Density")
+    #        p.set_zlim(zlim[f][0], zlim[f][1])
+    #        width = max(projwidth, 2*data["r200kpc"])
+    #        pc.set_width(width, "kpc")
+    #        boundL = data["GasCoM"] - 0.5*width/pf["kpc"]
+    #        boundR = data["GasCoM"] + 0.5*width/pf["kpc"]
+    #        p.modify["sphere"](data["GasCoM"], data["r200kpc"]/pf["kpc"],
+    #                           circle_args=carg)
+    #        tpos1 = [boundL[x_dict[p.data.axis]], boundR[y_dict[p.data.axis]]]
+    #        tpos2 = [boundL[x_dict[p.data.axis]], boundL[y_dict[p.data.axis]]]
+    #        p.modify["point"](tpos1, title1, text_args=targ1)
+    #        p.modify["point"](tpos2, title2, text_args=targ2)
+    #        for i in range(data["NumberOfPop2Stars"]):
+    #            pos = na.array(
+    #                [data["pop2_position_x"][i], data["pop2_position_y"][i],
+    #                 data["pop2_position_z"][i]])
+    #            age_bin = na.searchsorted(ages, data["pop2_age"][i])
+    #            ar = {'color':colors[age_bin]}
+    #            inside = na.logical_and(pos>boundL, pos<boundR).all()
+    #            if inside:
+    #                p.modify["marker"](pos, marker="*", plot_args=ar)
+    #pc.save("%s-halo%4.4d" % (pf,num))
+    #del pc
+    if rank == 0 and data["NumberOfPop2Stars"] > 0:
+        nullfmt = NullFormatter()
+        # Definitions for axes
+        left, width = 0.175, 0.6
+        bottom, height = 0.1, 0.625
+        left_h = left+width
+        #bottom_h = left+width-0.025
+        bottom_h = bottom+height
+
+        rect_scatter = [left, bottom, width, height]
+        rect_histx = [left, bottom_h, width, 0.175]
+        rect_histy = [left_h, bottom, 0.2, height]
+
+        plt.figure(1, figsize=(6,6))
+
+        axScatter = plt.axes(rect_scatter)
+        axHistx = plt.axes(rect_histx)
+        axHisty = plt.axes(rect_histy)
+
+        # No labels
+        axHistx.xaxis.set_major_formatter(nullfmt)
+        axHisty.yaxis.set_major_formatter(nullfmt)
+
+        # Scatter plot (only include massive star clusters)
+        data["pop2_age"] = data["pop2_age"][massive]
+        data["pop2_metallicity"] = data["pop2_metallicity"][massive]
+        data["pop2_mass"] = data["pop2_mass"][massive]
+        isort = na.argsort(data["pop2_mass"])[::-1]
+        x = data["pop2_age"][isort]
+        y = data["pop2_metallicity"][isort]
+        area = data["pop2_mass"][isort] / area_normalize
+        axScatter.scatter(x,y,
+                          color='#888888', s=area, alpha=0.75, edgecolors='k')
+        axScatter.set_ylim(-4.0, 0.0)
+        xl = axScatter.get_xlim()
+        yl = axScatter.get_ylim()
+        xl = (150.0, 1.01*tnow)
+        axScatter.set_xlim(xl[0], xl[1])
+        axScatter.set_xlabel(r'Time [Myr]')
+        axScatter.set_ylabel(r'[Z/H]')
+        axScatter.yaxis.set_label_coords(-0.19, 0.5)
+        axScatter.grid(True)
+        margin = [0.02 * (xl[1]-xl[0]), 0.02 * (yl[1]-yl[0])]
+        #txt = 'Halo %d' % num
+        #axScatter.text(xl[1]-margin[0], yl[0]+margin[1], txt, ha="right",
+        #               size='small')
+
+        # Draw log-linear "fits" to merged stellar populations for halo 0
+        if num == 0:
+            axScatter.plot([450,750], [-3.5, -2.5], "k--")
+            axScatter.plot([450,750], [-3.2, -2.1], "k--")
+
+        # Plot sample masses
+        props = dict(edgecolors='k')
+        masses = [4,3]
+        rbig = (10.0**masses[0]/area_normalize/na.pi)**0.5
+        for m in masses:
+            label = r'$10^{%d}$' % m
+            rthis = (10.0**m/area_normalize/na.pi)**0.5
+            axScatter.scatter(xl[0]+3*margin[0],
+                              yl[1]-3*margin[1] + (rbig-rthis)*(yl[1]-yl[0])/400,
+                              c="#FFFFFF",
+                              s=10.0**m/area_normalize,
+                              label=label, **props)
+        axScatter.set_xlim(xl[0], xl[1])
+        axScatter.set_ylim(yl[0], yl[1])
+
+        # Histograms
+        nmajor = 5
+        width = data["pop2_age_bin"][1] - data["pop2_age_bin"][0]
+        bin = data["pop2_age_bin"][0:-1]
+        height = 1e-6*data["pop2_age_pdf"] / width
+        nearest_power10 = int(na.log10(height.max()))-1
+        height /= 10.0**nearest_power10
+#        bin = 0.5 * (data["pop2_age_bin"][0:-1] + data["pop2_age_bin"][1:])
+        axHistx.bar(bin, height, width, color="#AAAAAA")
+        axHistx.set_ylabel(r'SFR [$10^{%d} M_\odot$/yr]' % nearest_power10)
+        axHistx.yaxis.set_label_coords(-0.19, 0.5)
+        x0,x1 = axScatter.get_xlim()
+        axHistx.set_xlim( x0,x1 )
+        
+        #axHistx.yaxis.set_major_locator(LinearLocator(numticks=nmajor))
+        #ticks = axHistx.yaxis.get_major_ticks()
+        #ticks[0].label1On = False
+
+        width = data["pop2_metallicity_bin"][1] - data["pop2_metallicity_bin"][0]
+        bin = data["pop2_metallicity_bin"][0:-1]
+#        bin = 0.5 * (data["pop2_metallicity_bin"][0:-1] +
+#                     data["pop2_metallicity_bin"][1:])
+        axHisty.barh(bin, data["pop2_metallicity_pdf"], height=width,
+                     color="#AAAAAA")
+        axHisty.set_xlabel(r'M$_\star$ [$M_\odot$]')
+        axHisty.set_ylim( axScatter.get_ylim() )
+        axHisty.xaxis.set_label_position('top')
+        axHisty.xaxis.set_major_locator(LinearLocator(numticks=nmajor))
+        ticks = axHisty.xaxis.get_major_ticks()
+        for t in ticks:
+            t.label1On = False
+            t.label2On = True
+            t.label2.set_rotation(270)
+        ticks[0].label2On = False        
+
+        x0,x1 = axHistx.get_xlim()
+        y0,y1 = axHistx.get_ylim()
+        p1 = axHistx.twiny()
+        p1.set_xlabel("Redshift")
+        tnow = 538.0 * ((_pf["CosmologyCurrentRedshift"]+1)/10.0)**(-1.5)
+        tticks = []
+        for tick in zticks:
+            thisz = 538.0 * ((1.0 + tick)/10)**(-1.5)
+            tticks.append(thisz)
+        #p1.set_xscale('log')
+        p1.set_xticks(tticks)
+        p1.set_xlim(x0,x1)
+        p1.set_xticklabels(["%d" % e for e in zticks])
+        axHistx.set_xlim(x0,x1)
+        axHistx.set_ylim(y0,y1)
+        axHistx.yaxis.set_major_locator(LinearLocator(numticks=nmajor))
+        ticks = axHistx.yaxis.get_major_ticks()
+        ticks[0].label1On = False
+        name = "Halo %d" % num
+        if name != None:
+            axHistx.text(x0+0.02*(x1-x0),y1-0.05*(y1-y0),name,
+                         ha="left",va="top")
+
+        plt.savefig("%s-halo%4.4d-stars.eps" % (_pf,num))
+        plt.clf()
+
+    return
+    
+########################################################################
+
+if append_output is False and rank == 0:
+    if os.path.exists(output_file):
+        os.remove(output_file)
+
+# Make projection of fields
+fields = ["Density"]
+#pc = PlotCollection(pf, center=[0.5]*3)
+#for f in fields:
+#    for dim in range(3):
+#        pc.add_projection(f, dim, "Density")
+#del pc
+
+# Load halos
+cycle = int(pf.get_parameter("InitialCycleNumber"))
+halos = readhalos(cycle=cycle, maxhalos=nhalos)
+for i in range(len(halos["mass"])):
+    one_halo(pf, i, halos["pos"][i], halos["mass"][i], output_file)
+    plot_halo(pf, i, output_file)
+
+MPI.Finalize()
+

posts/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2.py.desc

+Test paste.  Dwarf galaxy analysis script.