Devin Silvia avatar Devin Silvia committed 11fa65b

updated carbon files to follow the new roman numeral convention.

Comments (0)

Files changed (3)

dengo/carbon_cooling.py

+"""
+Author: Devin Silvia <devin.silvia@gmail.com>
+Affiliation: UC Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2012 Matthew Turk.  All Rights Reserved.
+
+  This file is part of the dengo package.
+
+  This file is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
 from reaction_classes import Species, ion_cooling_rate, species_registry
+import docutils.utils.roman as roman
 
-for i in range(7):
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'C'
+atomicNumber = 6
+atomicWeight = 12
+nIons = atomicNumber + 1
+
+for i in range(nIons):
     ion_state = i + 1
-    speciesName = "c_%s" % ion_state
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
     # Check if the species already exists
     # in the species registry, if it does
     # we don't want to create it again
     if (speciesName in species_registry) == False:
-        s = Species(speciesName, 12, i)
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
     else:
         s = species_registry[speciesName]
     ion_cooling_rate(s)

dengo/carbon_rates.py

 """
 
 from reaction_classes import Species, chianti_rate, species_registry
+import docutils.utils.roman as roman
 
-for i in range(7):
+# Note: the atomic/species properties have to be hard-coded for this
+# we may want to come up with a better solution here...
+atomicSymbol = 'C'
+atomicNumber = 6
+atomicWeight = 12
+nIons = atomicNumber + 1
+
+for i in range(nIons):
     ion_state = i + 1
-    speciesName = "c_%s" % ion_state
+    speciesName = "%s%s" %(atomicSymbol, roman.toRoman(ion_state))
     # Check if the species already exists
     # in the species registry, if it does
     # we don't want to create it again
     if (speciesName in species_registry) == False:
-        s = Species(speciesName, 12, i)
+        s = Species(speciesName, atomicNumber, atomicWeight, i)
     else:
         s = species_registry[speciesName]
 
-    if ion_state != 7:
+    if ion_state != nIons:
         # we need to do this to make sure the 'ion_state + 1' species
         # exists when chianti_rate is called
-        speciesNamePlusOne = "c_%s" % (ion_state+1)
+        speciesNamePlusOne = "%s%s" % (atomicSymbol, roman.toRoman(ion_state+1))
         if (speciesNamePlusOne in species_registry) == False:
-            splusone = Species(speciesNamePlusOne, 12, i+1)
+            splusone = Species(speciesNamePlusOne, atomicNumber, atomicWeight, i+1)
+
     chianti_rate(s)

examples/run_carbon_network.py

 
 # If only a subset of species are wanted put them here
 # and change the commented lines below
-want = ('c_5', 'c_6', 'c_7', 'de', 'ge')
+want = ('CIII', 'CIV', 'CV', 'de', 'ge')
 
 carbon = ChemicalNetwork()
 carbon.add_energy_term()
 
-for ca in cooling_registry.values():
-    # The following line can be used to specify a subset of species
-    #if all(sp.name in want for sp in ca.species):
-    if ca.name.startswith("c_"):
-       carbon.add_cooling(ca)
+# for ca in cooling_registry.values():
+#     # The following line can be used to specify a subset of species
+#     #if all(sp.name in want for sp in ca.species):
+#     if ca.name.startswith("C"):
+#        carbon.add_cooling(ca)
 
 for s in reaction_registry.values():
     # The following line can be used to specify a subset of species
     #if all(sp.name in want for sp in s.species):
-    if s.name.startswith("c_"):
+    if s.name.startswith("C"):
         carbon.add_reaction(s)
 
 # This defines the temperature range for the rate tables
 
 # Generate initial conditions (switch to False to disable this)
 generate_initial_conditions = True
+start_neutral = False
+initial_state = 'CVII' # if using start_neutral = True, this should be the _1 state
+                      # if starting in a close to CIE state, this should be the
+                      # ion species that will be closest for one for the given Temp.
+                      # note, this will get hairy when we're not testing a single T
 
 if generate_initial_conditions:
-    import numpy as na
+    import numpy as np
     NCELLS = 1
-    density = 1.0
-    init_array = na.ones(NCELLS) 
-    X = 1e-6
-
+    density = 1
+    init_array = np.ones(NCELLS)
     init_values = dict()
     init_values['density'] = density * init_array
-    initial_state = 'c_1'
     init_values[initial_state] = init_array.copy() # use conservation to set this below
-    # populate initial fractional values for the other species
-    for s in sorted(carbon.required_species):
-        if s.name != 'ge' and s.name != initial_state:
-            if s.name == 'de':
-                continue
-            else:
-                init_values[s.name] = X * init_array
-            init_values[initial_state] -= init_values[s.name]
-    init_values['de'] = init_array * 0.0
-    for s in sorted(carbon.required_species):
-        if s.name in ("ge", "de"): continue
-        init_values['de'] += init_values[s.name] * s.free_electrons
-        print "Adding %0.5e to electrons from %s" % (
-            (init_values[s.name] * s.free_electrons)[0], s.name)
-    print "Total de: %0.5e" % (init_values['de'][0])
+           
+    # set up initial temperatures values used to define ge
+    temperature = np.logspace(4, 6.7, NCELLS)
+    temperature[:] = 1e7; # need to remove this line for the above one to matter
+    init_values['T'] = temperature
+
+    if start_neutral:
+        X = 1e-6
+ 
+        # populate initial fractional values for the other species
+        for s in sorted(carbon.required_species):
+            if s.name != 'ge' and s.name != initial_state:
+                if s.name == 'de':
+                    continue
+                else:
+                    init_values[s.name] = X * init_array
+                init_values[initial_state] -= init_values[s.name]
+        init_values['de'] = init_array * 0.0
+        for s in sorted(carbon.required_species):
+            if s.name in ("ge", "de"): continue
+            init_values['de'] += init_values[s.name] * s.free_electrons
+            print "Adding %0.5e to electrons from %s" % (
+                (init_values[s.name] * s.free_electrons)[0], s.name)
+        print "Total de: %0.5e" % (init_values['de'][0])
+
+    else:
+        # Unless neutral is desired, we'll start in a perturbed CIE solution
+        import chianti.core as ch
+
+        # populate initial fractional values for the species
+        for s in sorted(carbon.required_species):
+            if s.name != 'ge' and s.name != initial_state:
+                if s.name == 'de':
+                    continue
+                else:
+                    ion = ch.ion(s.name, temperature=init_values['T'])
+                    ion.ioneqOne()
+                    ion_frac = ion.IoneqOne
+                    if ion_frac == 0.0:
+                        init_values[s.name] = tiny * init_array
+                    else:
+                        ion_frac[ion_frac < tiny] = tiny
+                        init_values[s.name] = ion_frac * init_array
+                
+                # add some random noise
+                init_values[s.name] += 0.5 * init_values[s.name] * np.random.randn(NCELLS)
+                # in case something went negative:
+                init_values[s.name][init_values[s.name] < tiny] = tiny
+                
+                # conservation...
+                init_values[initial_state] -= init_values[s.name]
+        
+        init_values[initial_state][init_values[initial_state] < tiny] = tiny
+        init_values['de'] = init_array * 0.0
+        for s in sorted(carbon.required_species):
+            if s.name in ("ge", "de"): continue
+            init_values['de'] += init_values[s.name] * s.free_electrons
+            print "Adding %0.5e to electrons from %s" % (
+                (init_values[s.name] * s.free_electrons)[0], s.name)
+        print "Total de: %0.5e" % (init_values['de'][0])
 
     # convert to masses to multiplying by the density factor and the species weight
     for s in carbon.required_species:
             init_values[s.name] *= (density) # this is still just number density
         else:
             init_values[s.name] *= (density * s.weight)
-    
+
     #compute new total density and number density
     density = 0.0
     number_density = 0.0
         if s.name == 'de': continue
         density += init_values[s.name][0]
 
-    # set up initial temperatures values used to define ge
-    temperature = na.logspace(4, 6.7, NCELLS)
-    temperature[:] = 1e7; # need to remove this line for the above one to matter
-    init_values['T'] = temperature
-
     # calculate ge (very crudely)
     gamma = 5.e0/3.e0
     init_values['ge'] = ((temperature * number_density * kboltz)
                          / (density * mh * (gamma - 1)))
-    
+
     # Write the initial conditions file
     create_initial_conditions(init_values, 'carbon')
 
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.