Merge branch 'Correct_additional_EAGLE_cooling_abundances' into 'master'

Correct additional eagle cooling abundances

Closes #509

See merge request !708

doc/RTD/source/SubgridModels/Basic/index.rst

+.. Basic sub-grid model
+   Matthieu Schaller, 20th December 2018
+
+
+Basic model
+===========
+
+
+Cooling: Analytic models
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Currently, we have 3 different simple cooling models (const-lambda, const-du
+and Compton). These are all based on analytic formulas and can be used
+to quickly understand how the cooling interacts with the rest of the
+code before moving to more complex models.
+
+Equations
+---------
+
+The first table compares the different analytical cooling while the next ones
+are specific to a given cooling.  The quantities are the internal energy (\\( u
+\\)), the density \\( rho \\), the element mass fraction (\\( X_i \\)), the
+cooling function (\\(\\Lambda\\), the proton mass (\\( m_H \\)) and the time
+step condition (\\( t\_\\text{step}\\)).  If not specified otherwise, all
+cooling contains a temperature floor avoiding negative temperature.
+
+.. csv-table:: Analytical Cooling
+   :header: "Variable", "Const-Lambda", "Const-du"
+
+   "\\( \\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \\)", "\\( -\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \\)", "const"
+   "\\( \\Delta t\_\\text{max} \\)", "\\( t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "\\( t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)"
+
+TODO: Add description of the parameters and units.
+
+TODO: Add Compton cooling model
+
+How to Implement a New Cooling
+------------------------------
+
+The developer should provide at least one function for:
+ * writing the cooling name in HDF5
+ * cooling a particle
+ * the maximal time step possible
+ * initializing a particle
+ * computing the total energy radiated by a particle
+ * initializing the cooling parameters
+ * printing the cooling type
+
+For implementation details, see ``src/cooling/none/cooling.h``
+
+See :ref:`new_option` for the full list of changes required.

doc/RTD/source/SubgridModels/EAGLE/index.rst

+.. EAGLE sub-grid model
+   Matthieu Schaller, 20th December 2018
+
+
+EAGLE model
+===========
+
+Chemical tracers
+~~~~~~~~~~~~~~~~
+
+The gas particles in the EAGLE model carry metal abundance information
+in the form of metal mass fractions. We follow the following 9
+elements: `H`, `He`, `C`, `N`, `O`, `Ne`, `Mg`, `Si` and `Fe`. We
+additionally follow the total metal mass fraction `Z`. This is
+typically larger than the sum of the 7 metals that are individually
+traced since this will also contain the contribution of all the
+elements that are not individually followed.
+
+As part of the diagnostics, we additionally trace the elements coming
+from the different stellar evolution channels. We store for each
+particle the total mass coming from all the SNIa that enriched that
+particle and the metal mass fraction from SNIa. This is the fraction
+of the *total* gas mass that is in the form of metals originating from
+SNIa stars. By construction this fraction will be smaller than the
+total metal mass fraction. The same tracers exist for the SNII and AGB
+channels. Finally, we also compute the iron gas fraction from
+SNIa. This it the fraction of the *total* gas mass that is made of
+iron originating from SNIa explosions. 
+
+We finally also compute the smoothed versions of the individual
+element mass fractions, of the total metal mass fractions, and of the
+iron gas fraction from SNIa.
+
+The chemistry module in ``src/chemistry/EAGLE`` includes all the arrays
+that are added to the particles and the functions used to compute the
+smoothed elements.
+
+In the snapshots, we output:
+
++------------------------------+--------------------------------+-----------+-----------------------------+
+| Name                         | Description                    | Units     | Comments                    | 
++==============================+================================+===========+=============================+
+| ``ElementAbundance``         | | Fraction of the gas mass     | [-]       | | Array of length           |
+|                              | | in the different elements    |           | | 9 for each particle       |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``SmoothedElementAbundance`` | | Fraction of the gas mass     | [-]       | | Array of length           |
+|                              | | in the different elements    |           | | 9 for each particle       |
+|                              | | smoothed over SPH neighbours |           |                             |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``Metallicity``              | | Fraction of the gas mass     | [-]       |                             |
+|                              | | in *all* metals              |           |                             |
++------------------------------+--------------------------------+-----------+-----------------------------+
+| ``SmoothedMetallicity``      | | Fraction of the gas mass     | [-]       |                             |
+|                              | | in *all* metals              |           |                             |
+|                              | | smoothed over SPH neighbours |           |                             |
++------------------------------+--------------------------------+-----------+-----------------------------+
+
+Cooling: Wiersma+2008a
+~~~~~~~~~~~~~~~~~~~~~~
+
+Particle tracers
+~~~~~~~~~~~~~~~~
+
+Star formation: Schaye+2008
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Stellar enrichment: Wiersma+2008b
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Supernova feedback: Schaye+2012
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Black-hole creation
+~~~~~~~~~~~~~~~~~~~
+
+Black-hole accretion
+~~~~~~~~~~~~~~~~~~~~
+
+AGN feedback
+~~~~~~~~~~~~

doc/RTD/source/Cooling/index.rst → doc/RTD/source/SubgridModels/GEAR/index.rst

-.. Equation of State
-   Loic Hausammann, 7th April 2018
+.. GEAR sub-grid model
+   Matthieu Schaller, 20th December 2018
 
-.. _cooling:
 
-Cooling
-=======
+GEAR model
+===========
 
-Currently, we have 5 different cooling (EAGLE, Grackle, const-lambda, const-du
-and none).  Three of them are easily solved analytically (const-lambda,
-const-du and none) while the two last requires complex chemical networks.
 
-
-Equations
----------
-
-The first table compares the different analytical cooling while the next ones
-are specific to a given cooling.  The quantities are the internal energy (\\( u
-\\)), the density \\( rho \\), the element mass fraction (\\( X_i \\)), the
-cooling function (\\(\\Lambda\\), the proton mass (\\( m_H \\)) and the time
-step condition (\\( t\_\\text{step}\\)).  If not specified otherwise, all
-cooling contains a temperature floor avoiding negative temperature.
-
-.. csv-table:: Analytical Cooling
-   :header: "Variable", "Const-Lambda", "Const-du", "None"
-
-   "\\( \\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \\)", "\\( -\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \\)", "const", "0"
-   "\\( \\Delta t\_\\text{max} \\)", "\\( t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "\\( t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \\)", "None"
-
-
-Grackle
-~~~~~~~
+Cooling: Grackle
+~~~~~~~~~~~~~~~~
    
 Grackle is a chemistry and cooling library presented in `B. Smith et al. 2016 <https://arxiv.org/abs/1610.09591>`_ 
 (do not forget to cite if used).  Four different modes are available:
@@ -57,24 +35,3 @@ Finish with ``make machine-linux-gnu; make && make install``.
 If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
 
 You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
-
-Eagle
-~~~~~
-
-TODO
-
-How to Implement a New Cooling
-------------------------------
-
-The developer should provide at least one function for:
- * writing the cooling name in HDF5
- * cooling a particle
- * the maximal time step possible
- * initializing a particle
- * computing the total energy radiated by a particle
- * initializing the cooling parameters
- * printing the cooling type
-
-For implementation details, see ``src/cooling/none/cooling.h``
-
-See :ref:`new_option` for the full list of changes required.

doc/RTD/source/SubgridModels/index.rst

+.. Subgrid Models
+   Matthieu Schaller, 20th December 2018
+
+.. _subgrid:
+   
+Galaxy Formation Subgrid Models
+===============================
+
+Multiple models are available in SWIFT. The 'Basic' model can
+be use as an empty canvas to be copied to create additional models.
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Available models:
+	      
+   Basic/index	      
+   EAGLE/index
+   GEAR/index

doc/RTD/source/index.rst

@@ -19,7 +19,7 @@ difference is the parameter file that will need to be adapted for SWIFT.
    ParameterFiles/index
    InitialConditions/index
    HydroSchemes/index
-   Cooling/index
+   SubgridModels/index
    EquationOfState/index
    ExternalPotentials/index
    NewOption/index

examples/CoolingBox/coolingBox.yml

@@ -70,8 +70,6 @@ EAGLEChemistry:
   InitAbundance_Magnesium: 5.907064e-4
   InitAbundance_Silicon:   6.825874e-4
   InitAbundance_Iron:      1.1032152e-3
-  CalciumOverSilicon:      0.0941736
-  SulphurOverSilicon:      0.6054160
 
 GearChemistry:
   InitialMetallicity: 0.01295

examples/CoolingRates/cooling_rates.yml

@@ -26,9 +26,6 @@ EAGLEChemistry:
   InitAbundance_Magnesium: 5.907064e-4
   InitAbundance_Silicon:   6.825874e-4
   InitAbundance_Iron:      1.1032152e-3
-  CalciumOverSilicon:      0.0941736
-  SulphurOverSilicon:      0.6054160
-
 
 EagleCooling:
   filename:                ./coolingtables/

examples/EAGLE_12/eagle_12.yml

@@ -68,8 +68,6 @@ EAGLEChemistry: 	    # Solar abundances
   InitAbundance_Magnesium:  5.907064e-4
   InitAbundance_Silicon:    6.825874e-4
   InitAbundance_Iron:       1.1032152e-3
-  CalciumOverSilicon:       0.0941736
-  SulphurOverSilicon:       0.6054160
 
 EagleCooling:
   filename:                ./coolingtables/

examples/EAGLE_25/eagle_25.yml

@@ -77,8 +77,6 @@ EAGLEChemistry: 	    # Solar abundances
   InitAbundance_Magnesium:  5.907064e-4
   InitAbundance_Silicon:    6.825874e-4
   InitAbundance_Iron:       1.1032152e-3
-  CalciumOverSilicon:       0.0941736
-  SulphurOverSilicon:       0.6054160
 
 EagleCooling:
   filename:                ./coolingtables/

examples/EAGLE_50/eagle_50.yml

@@ -70,8 +70,6 @@ EAGLEChemistry: 	    # Solar abundances
   InitAbundance_Magnesium:  5.907064e-4
   InitAbundance_Silicon:    6.825874e-4
   InitAbundance_Iron:       1.1032152e-3
-  CalciumOverSilicon:       0.0941736
-  SulphurOverSilicon:       0.6054160
 
 EagleCooling:
   filename:                ./coolingtables/

examples/EAGLE_6/eagle_6.yml

@@ -81,8 +81,6 @@ EAGLEChemistry: 	    # Solar abundances
   InitAbundance_Magnesium:  5.907064e-4
   InitAbundance_Silicon:    6.825874e-4
   InitAbundance_Iron:       1.1032152e-3
-  CalciumOverSilicon:       0.0941736
-  SulphurOverSilicon:       0.6054160
 
 EagleCooling:
   filename:                ./coolingtables/

examples/SantaBarbara/santa_barbara.yml

@@ -70,8 +70,6 @@ EAGLEChemistry:
   InitAbundance_Magnesium: 0.0
   InitAbundance_Silicon:   0.0
   InitAbundance_Iron:      0.0
-  CalciumOverSilicon:      0.0
-  SulphurOverSilicon:      0.0
 
 EagleCooling:
   filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/

examples/SmallCosmoVolume_cooling/small_cosmo_volume.yml

@@ -82,5 +82,3 @@ EAGLEChemistry:
   InitAbundance_Magnesium:  0.
   InitAbundance_Silicon: 0.
   InitAbundance_Iron: 0.
-  CalciumOverSilicon: 0.
-  SulphurOverSilicon: 0.

examples/parameter_example.yml

@@ -248,12 +248,15 @@ LambdaCooling:
   lambda_nH2_cgs:              1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
   cooling_tstep_mult:          1.0   # (Optional) Dimensionless pre-factor for the time-step condition.
 
+# Parameters of the EAGLE cooling model (Wiersma+08 cooling tables).
 EagleCooling:
-  filename:                ./coolingtables/  # Location of the Wiersma+08 cooling tables
-  reionisation_redshift:   11.5              # Redshift of Hydrogen re-ionization
-  He_reion_z_centre:       3.5               # Redshift of the centre of the Helium re-ionization Gaussian
-  He_reion_z_sigma:        0.5               # Spread in redshift of the  Helium re-ionization Gaussian
-  He_reion_ev_pH:          2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+  filename:                  ./coolingtables/  # Location of the Wiersma+08 cooling tables
+  reionisation_redshift:     8.5               # Redshift of Hydrogen re-ionization
+  He_reion_z_centre:         3.5               # Redshift of the centre of the Helium re-ionization Gaussian
+  He_reion_z_sigma:          0.5               # Spread in redshift of the  Helium re-ionization Gaussian
+  He_reion_ev_pH:            2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+  CalciumOverSiliconInSolar: 1.                # (Optional) Ratio of Ca/Si to use in units of solar. If set to 1, the code uses [Ca/Si] = 0, i.e. Ca/Si = 0.0941736.
+  SulphurOverSiliconInSolar: 1.                # (Optional) Ratio of S/Si to use in units of solar. If set to 1, the code uses [S/Si] = 0, i.e. S/Si = 0.6054160.
   
 # Cooling with Grackle 3.0
 GrackleCooling:
@@ -282,8 +285,6 @@ EAGLEChemistry:
   InitAbundance_Magnesium: 0.000        # Inital fraction of particle mass in Magnesium
   InitAbundance_Silicon:   0.000        # Inital fraction of particle mass in Silicon
   InitAbundance_Iron:      0.000        # Inital fraction of particle mass in Iron
-  CalciumOverSilicon:      0.0941736    # Constant ratio of Calcium over Silicon abundance
-  SulphurOverSilicon:      0.6054160    # Constant ratio of Sulphur over Silicon abundance
 
 # Structure finding options (requires velociraptor)
 StructureFinding:

src/chemistry/EAGLE/chemistry.h

@@ -188,12 +188,6 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
       data->initial_metal_mass_fraction[elem] =
           parser_get_param_float(parameter_file, buffer);
     }
-
-    /* Read the constant ratios */
-    data->calcium_over_silicon_ratio = parser_get_param_float(
-        parameter_file, "EAGLEChemistry:CalciumOverSilicon");
-    data->sulphur_over_silicon_ratio = parser_get_param_float(
-        parameter_file, "EAGLEChemistry:SulphurOverSilicon");
   }
 }
 

src/chemistry/EAGLE/chemistry_struct.h

@@ -45,12 +45,6 @@ struct chemistry_global_data {
 
   /*! Fraction of the particle mass in *all* metals at the start of the run */
   float initial_metal_mass_fraction_total;
-
-  /*! Constant ratio of Calcium over Silicium */
-  float calcium_over_silicon_ratio;
-
-  /*! Constant ratio of Calcium over Silicium */
-  float sulphur_over_silicon_ratio;
 };
 
 /**

src/cooling/EAGLE/cooling.c

@@ -466,7 +466,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
      (See cosmology theory document for the derivation) */
   const double delta_redshift = -dt * cosmo->H * cosmo->a_inv;
 
-  /* Get this particle's abundance ratios */
+  /* Get this particle's abundance ratios compared to solar
+   * Note that we need to add S and Ca that are in the tables but not tracked
+   * by the particles themselves.
+   * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
   float abundance_ratio[chemistry_element_count + 2];
   abundance_ratio_to_solar(p, cooling, abundance_ratio);
 
@@ -739,10 +742,10 @@ void cooling_init_backend(struct swift_params *parameter_file,
                           cooling->cooling_table_path);
   cooling->H_reion_z = parser_get_param_float(
       parameter_file, "EagleCooling:reionisation_redshift");
-  cooling->calcium_over_silicon_ratio = parser_get_param_float(
-      parameter_file, "EAGLEChemistry:CalciumOverSilicon");
-  cooling->sulphur_over_silicon_ratio = parser_get_param_float(
-      parameter_file, "EAGLEChemistry:SulphurOverSilicon");
+  cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float(
+      parameter_file, "EAGLECooling:CalciumOverSiliconInSolar", 1.f);
+  cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float(
+      parameter_file, "EAGLECooling:SulphurOverSiliconInSolar", 1.f);
   cooling->He_reion_z_centre =
       parser_get_param_float(parameter_file, "EagleCooling:He_reion_z_centre");
   cooling->He_reion_z_sigma =
@@ -866,6 +869,7 @@ void cooling_clean(struct cooling_function_data *cooling) {
   free(cooling->HeFrac);
   free(cooling->Therm);
   free(cooling->SolarAbundances);
+  free(cooling->SolarAbundances_inv);
 
   /* Free the tables */
   free(cooling->table.metal_heating);

src/cooling/EAGLE/cooling_rates.h

@@ -29,50 +29,68 @@
 #include "interpolate.h"
 
 /**
- * @brief Calculate ratio of particle element abundances
- * to solar abundance.
-
- * Multiple if statements are necessary because order of elements
- * in tables is different from chemistry_element enum.
- * Tables: H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe
- * Enum: H, He, C, N, O, Ne, Mg, Si, Fe
- * The order in ratio_solar is:
- * H, He, C, N, O, Ne, Mg, Si, Fe, S, Ca
- * Hence Fe, S, Ca need to be treated separately to be put in the
- * correct place in the output array.
- *
- * @param p Pointer to #part struct
- * @param cooling #cooling_function_data struct
- * @param ratio_solar Pointer to array or ratios
+ * @brief Compute ratio of mass fraction to solar mass fraction
+ * for each element carried by a given particle.
+ *
+ * The solar abundances are taken from the tables themselves.
+ *
+ * The EAGLE chemistry model does not track S and Ca. We assume
+ * that their abundance with respect to solar is the same as
+ * the ratio for Si.
+ * We optionally apply a correction if the user asked for a different
+ * ratio.
+ *
+ * We also re-order the elements such that they match the order of the
+ * tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
+ *
+ * The solar abundances table is arranged as
+ * [H, He, C, N, O, Ne, Mg, Si, S, C, Fe].
+ *
+ * @param p Pointer to #part struct.
+ * @param cooling #cooling_function_data struct.
+ * @param ratio_solar (return) Array of ratios to solar abundances.
  */
 __attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
     const struct part *p, const struct cooling_function_data *cooling,
     float ratio_solar[chemistry_element_count + 2]) {
 
-  /* compute ratios for all elements */
-  for (enum chemistry_element elem = chemistry_element_H;
-       elem < chemistry_element_count; elem++) {
-    if (elem == chemistry_element_Fe) {
-      /* NOTE: solar abundances have iron last with calcium and sulphur directly
-       * before, hence +2 */
-      ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
-                          cooling->SolarAbundances[elem + 2];
-    } else {
-      ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
-                          cooling->SolarAbundances[elem];
-    }
-  }
+  ratio_solar[0] = p->chemistry_data.metal_mass_fraction[chemistry_element_H] *
+                   cooling->SolarAbundances_inv[0];
+
+  ratio_solar[1] = p->chemistry_data.metal_mass_fraction[chemistry_element_He] *
+                   cooling->SolarAbundances_inv[1];
+
+  ratio_solar[2] = p->chemistry_data.metal_mass_fraction[chemistry_element_C] *
+                   cooling->SolarAbundances_inv[2];
+
+  ratio_solar[3] = p->chemistry_data.metal_mass_fraction[chemistry_element_N] *
+                   cooling->SolarAbundances_inv[3];
 
-  /* assign ratios for Ca and S, note positions of these elements occur before
-   * Fe */
-  ratio_solar[chemistry_element_count] =
-      p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
-      cooling->sulphur_over_silicon_ratio /
-      cooling->SolarAbundances[chemistry_element_count - 1];
-  ratio_solar[chemistry_element_count + 1] =
-      p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
-      cooling->calcium_over_silicon_ratio /
-      cooling->SolarAbundances[chemistry_element_count];
+  ratio_solar[4] = p->chemistry_data.metal_mass_fraction[chemistry_element_O] *
+                   cooling->SolarAbundances_inv[4];
+
+  ratio_solar[5] = p->chemistry_data.metal_mass_fraction[chemistry_element_Ne] *
+                   cooling->SolarAbundances_inv[5];
+
+  ratio_solar[6] = p->chemistry_data.metal_mass_fraction[chemistry_element_Mg] *
+                   cooling->SolarAbundances_inv[6];
+
+  ratio_solar[7] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7];
+
+  /* For S, we use the same ratio as Si */
+  ratio_solar[8] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7] *
+                   cooling->S_over_Si_ratio_in_solar;
+
+  /* For Ca, we use the same ratio as Si */
+  ratio_solar[9] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7] *
+                   cooling->Ca_over_Si_ratio_in_solar;
+
+  ratio_solar[10] =
+      p->chemistry_data.metal_mass_fraction[chemistry_element_Fe] *
+      cooling->SolarAbundances_inv[10];
 }
 
 /**
@@ -585,7 +603,7 @@ INLINE static double eagle_metal_cooling_rate(
 #endif
   }
 
-  const double abundance_ratio =
+  const double electron_abundance_ratio =
       H_plus_He_electron_abundance / solar_electron_abundance;
 
   /**********************/
@@ -596,22 +614,28 @@ INLINE static double eagle_metal_cooling_rate(
    * electron abundance to solar electron abundance then by the ratio of the
    * particle metal abundance to solar metal abundance. */
 
-  double lambda_metal[eagle_cooling_N_metal];
+  double lambda_metal[eagle_cooling_N_metal + 2] = {0.};
 
   if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) {
 
-    for (int elem = 0; elem < eagle_cooling_N_metal; elem++) {
+    /* Loop over the metals (ignore H and He) */
+    for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {
+
+      if (solar_ratio[elem] > 0.) {
 
-      lambda_metal[elem] =
-          interpolation_3d_no_x(cooling->table.metal_heating,   /* */
-                                elem, n_H_index, T_index,       /* */
-                                /*delta_elem=*/0.f, d_n_H, d_T, /* */
-                                eagle_cooling_N_metal,          /* */
-                                eagle_cooling_N_density,        /* */
-                                eagle_cooling_N_temperature);   /* */
+        /* Note that we do not interpolate along the x-axis
+         * (element dimension) */
+        lambda_metal[elem] =
+            interpolation_3d_no_x(cooling->table.metal_heating,   /* */
+                                  elem - 2, n_H_index, T_index,   /* */
+                                  /*delta_elem=*/0.f, d_n_H, d_T, /* */
+                                  eagle_cooling_N_metal,          /* */
+                                  eagle_cooling_N_density,        /* */
+                                  eagle_cooling_N_temperature);   /* */
 
-      lambda_metal[elem] *= abundance_ratio;
-      lambda_metal[elem] *= solar_ratio[elem + 2];
+        lambda_metal[elem] *= electron_abundance_ratio;
+        lambda_metal[elem] *= solar_ratio[elem];
+      }
 
 #ifdef TO_BE_DONE
       /* compute values at temperature gridpoints above and below input
@@ -636,19 +660,25 @@ INLINE static double eagle_metal_cooling_rate(
 
   } else {
 
-    for (int elem = 0; elem < eagle_cooling_N_metal; elem++) {
+    /* Loop over the metals (ignore H and He) */
+    for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) {
 
-      lambda_metal[elem] = interpolation_4d_no_x(
-          cooling->table.metal_heating,                /* */
-          elem, /*z_index=*/0, n_H_index, T_index,     /* */
-          /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
-          eagle_cooling_N_metal,                       /* */
-          eagle_cooling_N_loaded_redshifts,            /* */
-          eagle_cooling_N_density,                     /* */
-          eagle_cooling_N_temperature);                /* */
+      if (solar_ratio[elem] > 0.) {
 
-      lambda_metal[elem] *= abundance_ratio;
-      lambda_metal[elem] *= solar_ratio[elem + 2];
+        /* Note that we do not interpolate along the x-axis
+         * (element dimension) */
+        lambda_metal[elem] = interpolation_4d_no_x(
+            cooling->table.metal_heating,                /* */
+            elem - 2, /*z_index=*/0, n_H_index, T_index, /* */
+            /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */
+            eagle_cooling_N_metal,                       /* */
+            eagle_cooling_N_loaded_redshifts,            /* */
+            eagle_cooling_N_density,                     /* */
+            eagle_cooling_N_temperature);                /* */
+
+        lambda_metal[elem] *= electron_abundance_ratio;
+        lambda_metal[elem] *= solar_ratio[elem];
+      }
 
 #ifdef TO_BE_DONE
       /* compute values at temperature gridpoints above and below input
@@ -675,14 +705,14 @@ INLINE static double eagle_metal_cooling_rate(
   }
 
   if (element_lambda != NULL) {
-    for (int elem = 0; elem < eagle_cooling_N_metal; ++elem) {
-      element_lambda[elem + 2] = lambda_metal[elem];
+    for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
+      element_lambda[elem] = lambda_metal[elem];
     }
   }
 
   /* Sum up all the contributions */
   double Lambda_net = Lambda_free + Lambda_Compton;
-  for (int elem = 0; elem < eagle_cooling_N_metal; ++elem) {
+  for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) {
     Lambda_net += lambda_metal[elem];
   }
 

src/cooling/EAGLE/cooling_struct.h

@@ -65,20 +65,23 @@ struct cooling_function_data {
   /*! Internal energy bins */
   float *Therm;
 
-  /*! Solar mass fractions */
+  /*! Mass fractions of elements for solar abundances (from the tables) */