Commits

John Wise committed e85e964

Adding fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas.py

Comments (0)

Files changed (4)

html/fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas-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">import</span> <span class="nn">os</span><span class="o">,</span> <span class="nn">sys</span>
+<span class="kn">from</span> <span class="nn">yt.mods</span> <span class="kn">import</span> <span class="o">*</span>
+<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">na</span>
+
+<span class="k">if</span> <span class="nb">len</span><span class="p">(</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">2</span><span class="p">:</span>
+    <span class="k">print</span> <span class="s">&quot;usage: python </span><span class="si">%s</span><span class="s"> amr_file&quot;</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="mi">0</span><span class="p">]</span>
+    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">()</span>
+
+<span class="n">fname</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="k">if</span> <span class="ow">not</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">fname</span><span class="p">):</span>
+    <span class="k">print</span> <span class="s">&quot;</span><span class="si">%s</span><span class="s"> not found&quot;</span> <span class="o">%</span> <span class="n">fname</span>
+    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">()</span>
+
+<span class="n">pf</span> <span class="o">=</span> <span class="n">load</span><span class="p">(</span><span class="n">fname</span><span class="p">)</span>
+
+<span class="c"># Read &quot;yt hop&quot; output</span>
+<span class="n">hop_output</span> <span class="o">=</span> <span class="s">&quot;</span><span class="si">%s</span><span class="s">.hop&quot;</span> <span class="o">%</span> <span class="n">pf</span><span class="o">.</span><span class="n">basename</span>
+<span class="n">lines</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">hop_output</span><span class="p">)</span><span class="o">.</span><span class="n">readlines</span><span class="p">()</span>
+<span class="n">nhalos</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">lines</span><span class="p">)</span><span class="o">-</span><span class="mi">2</span>  <span class="c"># -2 header lines</span>
+<span class="n">centers</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">nhalos</span><span class="p">,</span><span class="mi">3</span><span class="p">))</span>
+<span class="n">dm_mass</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nhalos</span><span class="p">)</span>
+<span class="n">npart</span> <span class="o">=</span> <span class="n">na</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">nhalos</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="s">&#39;int&#39;</span><span class="p">)</span>
+<span class="n">i</span> <span class="o">=</span> <span class="mi">0</span>
+<span class="k">for</span> <span class="n">l</span> <span class="ow">in</span> <span class="n">lines</span><span class="p">:</span>
+    <span class="k">if</span> <span class="n">l</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s">&#39;#&#39;</span><span class="p">):</span> <span class="k">continue</span>
+    <span class="n">values</span> <span class="o">=</span> <span class="nb">map</span><span class="p">(</span><span class="nb">float</span><span class="p">,</span> <span class="n">l</span><span class="o">.</span><span class="n">split</span><span class="p">())</span>
+    <span class="n">centers</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">values</span><span class="p">[</span><span class="mi">4</span><span class="p">:</span><span class="mi">7</span><span class="p">]</span>
+    <span class="n">dm_mass</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">values</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
+    <span class="c">#radii[i] = values[13]</span>
+    <span class="n">npart</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">values</span><span class="p">[</span><span class="mi">2</span><span class="p">])</span>
+    <span class="n">i</span> <span class="o">+=</span> <span class="mi">1</span>
+<span class="k">del</span> <span class="n">lines</span>
+
+<span class="c"># Calculate rvir [cm-&gt;code] directly from mass</span>
+<span class="n">radii</span> <span class="o">=</span> <span class="p">(</span> <span class="mf">6.673e-8</span> <span class="o">*</span> <span class="p">(</span><span class="n">dm_mass</span> <span class="o">*</span> <span class="mf">1.989e33</span><span class="p">)</span> <span class="o">/</span> <span class="mf">100.0</span> <span class="o">/</span> <span class="n">pf</span><span class="o">.</span><span class="n">omega_matter</span> <span class="o">/</span> \
+          <span class="p">(</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">pf</span><span class="o">.</span><span class="n">current_redshift</span><span class="p">)</span><span class="o">**</span><span class="mi">3</span> <span class="o">/</span> \
+          <span class="p">(</span><span class="n">pf</span><span class="o">.</span><span class="n">hubble_constant</span> <span class="o">*</span> <span class="mi">100</span> <span class="o">/</span> <span class="mf">3.086e19</span><span class="p">)</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.</span><span class="o">/</span><span class="mi">3</span><span class="p">)</span> <span class="o">/</span> <span class="n">pf</span><span class="p">[</span><span class="s">&#39;cm&#39;</span><span class="p">]</span>
+
+<span class="c"># Look in these halos and record some data</span>
+<span class="c"># 1. gas fraction</span>
+<span class="c"># 2. cold, dense gas fraction (&gt;100 cm^-3, T &lt; T_vir / 2)</span>
+<span class="c"># 3. H2 mass fraction</span>
+
+<span class="n">output_name</span><span class="o">=</span> <span class="s">&quot;</span><span class="si">%s</span><span class="s">.fgas&quot;</span> <span class="o">%</span> <span class="n">pf</span><span class="o">.</span><span class="n">basename</span>
+<span class="n">output</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">output_name</span><span class="p">,</span> <span class="s">&quot;w&quot;</span><span class="p">)</span>
+<span class="n">output</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s">&quot;# Halo#   TotalMass   GasMass   GasFraction  DenseMass   AvgH2Frac   H2Mass</span><span class="se">\n</span><span class="s">&quot;</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="n">nhalos</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">centers</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">radii</span><span class="p">[</span><span class="n">i</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">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">fh2</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;H2I_Fraction&quot;</span><span class="p">,</span> <span class="s">&quot;CellMassCode&quot;</span><span class="p">)</span>
+    <span class="n">cd_cut</span> <span class="o">=</span> <span class="n">sp</span><span class="o">.</span><span class="n">cut_region</span><span class="p">([</span><span class="s">&quot;grid[&#39;Density&#39;] &gt; 1.673e-24&quot;</span><span class="p">,</span> <span class="s">&quot;grid[&#39;Temperature&#39;] &lt; 500&quot;</span><span class="p">])</span>
+    <span class="n">mdense</span> <span class="o">=</span> <span class="n">cd_cut</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">line</span> <span class="o">=</span> <span class="s">&quot;</span><span class="si">%d</span><span class="s"> </span><span class="si">%g</span><span class="s"> </span><span class="si">%g</span><span class="s"> </span><span class="si">%g</span><span class="s"> </span><span class="si">%g</span><span class="s"> </span><span class="si">%g</span><span class="s"> </span><span class="si">%g</span><span class="se">\n</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">mtot</span><span class="p">,</span> <span class="n">mgas</span><span class="p">,</span> <span class="n">mgas</span><span class="o">/</span><span class="n">mtot</span><span class="p">,</span> <span class="n">mdense</span><span class="p">,</span> <span class="n">fh2</span><span class="p">,</span> <span class="n">fh2</span><span class="o">*</span><span class="n">mgas</span><span class="p">)</span>
+    <span class="n">output</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">line</span><span class="p">)</span>
+    <span class="k">print</span> <span class="n">line</span><span class="p">[:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+    <span class="k">del</span> <span class="n">sp</span>
+<span class="n">output</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
+</pre></div>
+</body>
+</html>
 var inventory_data = [
  {
   "username": "John Wise <jwise@astro.princeton.edu>", 
-  "modtime": 1300885454, 
+  "modtime": 1308145368, 
+  "name": "fgas.py", 
+  "descr": "Routine to calculate halo properties", 
+  "htmlname": "html/fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas-py.html", 
+  "modified": "Wed Jun 15 09:42:48 2011", 
+  "lastmod_hash": "tip", 
+  "fullname": "posts/fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas.py"
+ }, 
+ {
+  "username": "John Wise <jwise@astro.princeton.edu>", 
+  "modtime": 1300899854.0, 
   "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", 
+  "modified": "Wed Mar 23 13:04:14 2011", 
+  "lastmod_hash": "ef7f9a83b8b09f4450c2f7ba56d3c75096f28afa", 
   "fullname": "posts/cc35ab50-a9a0-4c2b-b000-c5c2f47cc6b3-analysis2.py"
  }
 ];

posts/fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas.py

+import os, sys
+from yt.mods import *
+import numpy as na
+
+if len(sys.argv) != 2:
+    print "usage: python %s amr_file" % sys.argv[0]
+    sys.exit()
+
+fname = sys.argv[-1]
+if not os.path.exists(fname):
+    print "%s not found" % fname
+    sys.exit()
+
+pf = load(fname)
+
+# Read "yt hop" output
+hop_output = "%s.hop" % pf.basename
+lines = open(hop_output).readlines()
+nhalos = len(lines)-2  # -2 header lines
+centers = na.zeros((nhalos,3))
+dm_mass = na.zeros(nhalos)
+npart = na.zeros(nhalos, dtype='int')
+i = 0
+for l in lines:
+    if l.startswith('#'): continue
+    values = map(float, l.split())
+    centers[i,:] = values[4:7]
+    dm_mass[i] = values[1]
+    #radii[i] = values[13]
+    npart[i] = int(values[2])
+    i += 1
+del lines
+
+# Calculate rvir [cm->code] directly from mass
+radii = ( 6.673e-8 * (dm_mass * 1.989e33) / 100.0 / pf.omega_matter / \
+          (1.0 + pf.current_redshift)**3 / \
+          (pf.hubble_constant * 100 / 3.086e19)**2 )**(1./3) / pf['cm']
+
+# Look in these halos and record some data
+# 1. gas fraction
+# 2. cold, dense gas fraction (>100 cm^-3, T < T_vir / 2)
+# 3. H2 mass fraction
+
+output_name= "%s.fgas" % pf.basename
+output = open(output_name, "w")
+output.write("# Halo#   TotalMass   GasMass   GasFraction  DenseMass   AvgH2Frac   H2Mass\n")
+for i in range(nhalos):
+    sp = pf.h.sphere(centers[i,:], radii[i])
+    mtot = sp.quantities["TotalMass"]()
+    mgas = sp.quantities["TotalQuantity"]("CellMassMsun")[0]
+    fh2 = sp.quantities["WeightedAverageQuantity"]("H2I_Fraction", "CellMassCode")
+    cd_cut = sp.cut_region(["grid['Density'] > 1.673e-24", "grid['Temperature'] < 500"])
+    mdense = cd_cut.quantities["TotalQuantity"]("CellMassMsun")[0]
+    line = "%d %g %g %g %g %g %g\n" % (i, mtot, mgas, mgas/mtot, mdense, fh2, fh2*mgas)
+    output.write(line)
+    print line[:-1]
+    del sp
+output.close()

posts/fedc5a46-bf59-446a-821e-acf7633e5fc7-fgas.py.desc

+Routine to calculate halo properties