Depletion changes

src/synthesizer_grids/cloudy/README.md

@@ -59,7 +59,7 @@ python run_cloudy.py \
 
 ### Create synthesizer grid
 
-Finally, we need to create a new `synthesizer` grid object containing the cloudy computed lines and spectra. This is acheived by running `create_synthesizer_grid.py`.
+Finally, we need to create a new `synthesizer` grid object containing the cloudy computed lines and spectra. This is acheived by running `create_synthesizer_grid.py`. By default spectra are not extracted; if this is desired it necessary to add the `include_spectra` argument as below.
 
 ```
 python create_synthesizer_grid.py \
@@ -68,5 +68,6 @@ python create_synthesizer_grid.py \
   --cloudy-output-dir=/path/to/cloudy/outputs \
   --cloudy-paramfile=c23.01-sps \
   --cloudy-paramfile-extra=test_suite/reference_ionisation_parameter
+  --include_spectra \
 ```
 

src/synthesizer_grids/cloudy/create_cloudy_input_grid.py

@@ -11,9 +11,7 @@
 
 import numpy as np
 import yaml
-from synthesizer.abundances import (
-    Abundances,
-)
+from synthesizer.abundances import Abundances, depletion_models
 from synthesizer.exceptions import InconsistentParameter
 from synthesizer.grid import Grid
 from synthesizer.photoionisation import cloudy17, cloudy23
@@ -60,17 +58,47 @@ def create_cloudy_input(
         if len(k.split(".")) > 1:
             if k.split(".")[0] == "abundance_scalings":
                 kk = k.split(".")[1]
-                # convert to synthesizer standard
-                parameters["abundance_scalings"][kk] = v
+
+                if v is not None:
+                    # The value (v) here could be a bona fide string or number,
+                    # try to convert to a
+                    try:
+                        value = float(v)
+                    except ValueError:
+                        value = v
+
+                    # convert to synthesizer standard
+                    parameters["abundance_scalings"][kk] = value
+
+    # Initialise depletion model
+    if "depletion_model" in parameters.keys():
+        if parameters["depletion_model"] is None:
+            depletion_model = None
+        else:
+            if "depletion_scale" not in parameters.keys():
+                depletion_scale = 1.0
+            else:
+                depletion_scale = parameters["depletion_scale"]
+
+            depletion_model = getattr(
+                depletion_models, parameters["depletion_model"]
+            )(scale=depletion_scale)
+
+    else:
+        depletion_model = None
+
+    print(parameters["metallicities"])
+    print(parameters["reference_abundance"])
+    print(parameters["alpha_enhancement"])
+    print(parameters["abundance_scalings"])
 
     # Create synthesizer.Abundance object
     abundances = Abundances(
         metallicity=float(parameters["metallicities"]),
         reference=parameters["reference_abundance"],
         alpha=parameters["alpha_enhancement"],
         abundances=parameters["abundance_scalings"],
-        depletion_model=parameters["depletion_model"],
-        depletion_scale=parameters["depletion_scale"],
+        depletion_model=depletion_model,
     )
 
     # Define the ionisation parameter and geometry
@@ -112,7 +140,7 @@ def create_cloudy_input(
     elif parameters["ionisation_parameter_model"] == "fixed":
         ionisation_parameter = parameters["ionisation_parameter"]
 
-    # If the model is not regnoised raise an exception
+    # If the model is not recognised raise an exception
     else:
         raise InconsistentParameter(
             f"""ERROR: do not understand U model choice:

src/synthesizer_grids/cloudy/create_synthesizer_grid.py

@@ -66,6 +66,12 @@
     "stop_column_density": "column_densities",
     "alpha_enhancement": "alpha_enhancements",
     "reference_ionisation_parameter": "reference_ionisation_parameters",
+    "abundance_scalings.nitrogen_to_oxygen": (
+        "abundance_scalings.nitrogen_to_oxygen"
+    ),
+    "abundance_scalings.carbon_to_oxygen": (
+        "abundance_scalings.carbon_to_oxygen"
+    ),
 }
 
 
@@ -187,8 +193,6 @@ def create_empty_grid(
         data=weight_var,
     )
 
-    print(incident_grid._extract_axes)
-
     # Now write out the grid axes, we first do the incident grid axes so we
     # can extract their metadata and then any extras
     for axis, extract_axis in zip(
@@ -404,7 +408,7 @@ def check_cloudy_runs(
                     )
 
                 # Print the indices and the last line of the cloudy out file
-                print(incident_index, photoionisation_index, last_line)
+                print("    ", incident_index, photoionisation_index, last_line)
 
     number_of_failed_models = len(failed_list)
 
@@ -662,14 +666,14 @@ def add_lines(
     lines["luminosity"] = np.empty((*new_shape, nlines))
 
     # ... but only save continuum values if spectra are provided.
-    if calculate_continuum:
-        continuum_quantities = [
-            "transmitted",
-            "incident",
-            "nebular_continuum",
-            "total_continuum",
-        ]
+    continuum_quantities = [
+        "transmitted",
+        "incident",
+        "nebular_continuum",
+        "total_continuum",
+    ]
 
+    if calculate_continuum:
         spectra_ = {}
 
         # Calculate incideted_contiuum
@@ -693,6 +697,13 @@ def add_lines(
         for continuum_quantity in continuum_quantities:
             lines[continuum_quantity] = np.empty((*new_shape, nlines))
 
+    # If we're not using spectra, we still need to create output arrays for
+    # the continuum quantities, but this time they are zeros.
+    else:
+        # Define output arrays
+        for continuum_quantity in continuum_quantities:
+            lines[continuum_quantity] = np.zeros((*new_shape, nlines))
+
     # Array for recording cloudy failures
     failures = np.zeros(new_shape)
 
@@ -779,9 +790,8 @@ def add_lines(
     lines["luminosity"] *= erg / s
 
     # If continuum values are calculated add units
-    if calculate_continuum:
-        for continuum_quantity in continuum_quantities:
-            lines[continuum_quantity] *= erg / s / Hz
+    for continuum_quantity in continuum_quantities:
+        lines[continuum_quantity] *= erg / s / Hz
 
     # Write the lines out
     new_grid.write_lines(
@@ -825,7 +835,7 @@ def add_lines(
         "--include-spectra",
         action="store_true",
         help="Should the spectra be included in the grid?",
-        default=True,
+        default=False,
         required=False,
     )
 
@@ -855,9 +865,9 @@ def add_lines(
 
     print(" " * 80)
     print("-" * 80)
-    print(incident_grid_name)
-    print(cloudy_param_file)
-    print(extra_cloudy_param_file)
+    print("incident_grid_name:", incident_grid_name)
+    print("cloudy_param_file:", cloudy_param_file)
+    print("extra_cloudy_param_file:", extra_cloudy_param_file)
 
     # Check for extensions
     # Create _file and _name versions (with and w/o extensions)
@@ -895,8 +905,6 @@ def add_lines(
         new_grid_name += "-" + extra_cloudy_param_name
 
     new_grid_file = new_grid_name + ".hdf5"
-    print(cloudy_param_name, cloudy_param_file)
-    print(new_grid_name, new_grid_file)
 
     # Open the incident grid using synthesizer
     incident_grid = Grid(
@@ -935,8 +943,12 @@ def add_lines(
         fixed_photoionisation_params |= photoionisation_fixed_params_
         variable_photoionisation_params |= photoionisation_variable_params_
 
-    print(fixed_photoionisation_params)
-    print(variable_photoionisation_params)
+    print("fixed photoionisation parameters:")
+    for k, v in fixed_photoionisation_params.items():
+        print(f"    {k}: {v}")
+    print("variable photoionisation parameters:")
+    for k, v in variable_photoionisation_params.items():
+        print(f"    {k}: {v}")
 
     # If we have photoionisation parameters that vary we need to calculate the
     # model list
@@ -977,7 +989,6 @@ def add_lines(
         new_axes += photoionisation_axes
         new_axes_values |= variable_photoionisation_params
 
-    print(new_axes)
     for key, value in new_axes_values.items():
         print(key, value)
 
@@ -1006,7 +1017,7 @@ def add_lines(
     # recording failures with zeros in the spectra and line quantities. The
     # code will also save an array of failure flags.
 
-    # Now add spectra
+    # Now add spectra.
     if include_spectra:
         lam, spectra = add_spectra(
             new_grid,

src/synthesizer_grids/cloudy/params/c23.01-blr-test.yaml

@@ -26,7 +26,7 @@ ionisation_parameter: [0.01, 0.1]  # value of ionisation
 hydrogen_density: [1.0e+9, 1.0e+10]  
 
 # STOPPING CRITERIA
-stop_column_density: [22, 23]  # stopping column density
+stop_column_density: [23, 25]  # stopping column density
 
 
 # MISC COMMANDS

src/synthesizer_grids/cloudy/params/c23.01-nlr-test.yaml

@@ -22,7 +22,7 @@ geometry: planeparallel
 
 # IONISATION PARAMETER 
 ionisation_parameter_model: fixed  # which ionisation parameter model to use. `ref` assumes a varying ionisation parameter at a fixed reference age and metallicity
-ionisation_parameter: [0.001, 0.01, 0.1]  # value of ionisation 
+ionisation_parameter: [0.001, 0.01]  # value of ionisation 
 # DENSITY
 hydrogen_density: [1.0e+3, 1.0e+4]  
 

src/synthesizer_grids/cloudy/params/c23.01-nlr.yaml

@@ -22,7 +22,7 @@ geometry: planeparallel
 
 # IONISATION PARAMETER 
 ionisation_parameter_model: fixed  # which ionisation parameter model to use. `ref` assumes a varying ionisation parameter at a fixed reference age and metallicity
-ionisation_parameter: [0.001, 0.01, 0.1]  # value of ionisation 
+ionisation_parameter: [0.0001, 0.001, 0.01, 0.1]  # value of ionisation 
 # DENSITY
 hydrogen_density: [1.0e+2, 1.0e+3, 1.0e+4, 1.0e+5]  
 

src/synthesizer_grids/cloudy/params/c23.01-sps-UnH.yaml

@@ -0,0 +1,45 @@
+# Vary just U and n_H
+
+# CLOUDY VERSION
+cloudy_version: c23.01
+
+# ABUNDANCE PATTERN AND DEPLETION
+reference_abundance: GalacticConcordance # the reference abundance pattern to assume
+alpha_enhancement: 0.0  # alpha element enhancement
+abundance_scalings:
+  nitrogen: GalacticConcordance
+  carbon: GalacticConcordance
+depletion_model: Jenkins2009_Gunasekera2021 # the depletion model.
+depletion_scale: 0.5 # the depletion scale factor. For linearly scaled depletion models (e.g. CloudyClassic, Gutkin2016) this is by default 1.0. For Jenkins2009 this should be 0.5 (fstar).
+
+# GRAINS
+grains: Orion # include ISM grains
+
+# GEOMETRY
+geometry: spherical
+
+# IONISATION PARAMETER 
+ionisation_parameter_model: fixed  # which ionisation parameter model to use. `ref` assumes a varying ionisation parameter at a fixed reference age and metallicity
+ionisation_parameter: [0.0001, 0.001, 0.01, 0.1]  
+
+# DENSITY
+hydrogen_density: [1.0e+1, 1.0e+2, 1.0e+3, 1.0e+4]
+
+# STOPPING CRITERIA
+stop_T: 500  # stopping temperature
+stop_efrac: -2  # limiting ratio of electron to H densities
+
+# MISC COMMANDS
+CMB:  # include CMB heating
+T_floor: 100  # lower gas temperature floor
+cosmic_rays: true  # flag for inclusion of cosmic ray heating
+# covering_factor: 1.0  #
+radius: 0.01  # log of the radius of the cloud, in parsecs
+turbulence: 100 # turbulence
+z: 0.0  # redshift
+
+# OUTPUT COMMANDS
+resolution: 1.0 # energy resolution relative to the default
+output_cont: true
+output_linelist: linelist-standard.dat
+iterate_to_convergence: true

src/synthesizer_grids/cloudy/params/c23.01-sps-grid.yaml

@@ -40,6 +40,6 @@ z: 0.0  # redshift
 
 # OUTPUT COMMANDS
 resolution: 1.0 # energy resolution relative to the default
-output_cont: false
+output_cont: true
 output_linelist: linelist-standard.dat
 iterate_to_convergence: true

src/synthesizer_grids/cloudy/params/c23.01-sps.yaml

@@ -25,8 +25,8 @@ reference_age: 1.0e+6  # log10(age/year) at which to set reference ionisation pa
 reference_metallicity: 0.01  # metallicity at which to set reference ionisation parameter, for U_model='ref'
 
 # DENSITY
-hydrogen_density: 1.0e+3
-constant_density: True
+hydrogen_density: 3.16e+2 # hydrogen density in cm^-3 (==10^2.5)
+constant_density: True  # flag for constant density
 
 # STOPPING CRITERIA
 stop_T: 500  # stopping temperature

src/synthesizer_grids/cloudy/params/test_suite/CO.yaml

@@ -0,0 +1,5 @@
+abundance_scalings:
+  carbon:
+  carbon_to_oxygen: [-1.5, -1.0, -0.34, 0.0] # scaling for Carbon, either float relative to Solar or string defining the in-built function to use
+
+

src/synthesizer_grids/cloudy/params/test_suite/NO.yaml

@@ -0,0 +1,3 @@
+abundance_scalings:
+  nitrogen: 
+  nitrogen_to_oxygen: [-2.0, -1.5, -1.03, -0.5] 

src/synthesizer_grids/cloudy/params/test_suite/fixed_ionisation_parameter.yaml

@@ -1,4 +1,4 @@
 # IONISATION PARAMETER 
-geometry: planeparallel
+geometry: sphere
 ionisation_parameter_model: fixed  # which ionisation parameter model to use. `ref` assumes a varying ionisation parameter at a fixed reference age and metallicity
 ionisation_parameter: [0.0001, 0.001, 0.01, 0.1] # value of ionisation parameter at reference value, for U_model='ref'

src/synthesizer_grids/cloudy/run_cloudy.py

@@ -110,13 +110,13 @@
     # the models we want to run
 
     # If a list is provided open the file containing the list
-    if args.list_file is not False:
+    if args.list_file is not None:
         incident_indices, photoionisation_indices = np.loadtxt(
             args.list_file, dtype=int
         )
 
         # If an index is also provided just choose this model
-        if args.list_index is not False:
+        if args.list_index is not None:
             incident_indices = [incident_indices[args.list_index]]
             photoionisation_indices = [
                 photoionisation_indices[args.list_index]
@@ -154,6 +154,9 @@
         incident_indices = incident_indices.flatten()
         photoionisation_indices = photoionisation_indices.flatten()
 
+    print(incident_indices)
+    print(photoionisation_indices)
+
     # Loop over the list of indices and run each one
     for incident_index, photoionisation_index in zip(
         incident_indices, photoionisation_indices