src/engine_strays.c

@@ -33,6 +33,13 @@
 /* Local headers. */
 #include "proxy.h"
 
+#ifdef WITH_MPI
+/* Number of particle types to wait for after launching the proxies. We have
+   parts, xparts, gparts, sparts, bparts and sinks to exchange, hence 6 types.
+ */
+#define MPI_REQUEST_NUMBER_PARTICLE_TYPES 6
+#endif
+
 /**
  * @brief Exchange straying particles with other nodes.
  *
@@ -57,18 +64,22 @@
  * @param ind_bpart The foreign #cell ID of each bpart.
  * @param Nbpart The number of stray bparts, contains the number of bparts
  *        received on return.
+ * @param offset_sinks The index in the sinks array as of which the foreign
+ *        parts reside (i.e. the current number of local #sink).
+ * @param ind_sink The foreign #cell ID of each sink.
+ * @param Nsink The number of stray sinks, contains the number of sinks
+ *        received on return.
  *
  * Note that this function does not mess-up the linkage between parts and
  * gparts, i.e. the received particles have correct linkeage.
  */
-void engine_exchange_strays(struct engine *e, const size_t offset_parts,
-                            const int *restrict ind_part, size_t *Npart,
-                            const size_t offset_gparts,
-                            const int *restrict ind_gpart, size_t *Ngpart,
-                            const size_t offset_sparts,
-                            const int *restrict ind_spart, size_t *Nspart,
-                            const size_t offset_bparts,
-                            const int *restrict ind_bpart, size_t *Nbpart) {
+void engine_exchange_strays(
+    struct engine *e, const size_t offset_parts, const int *restrict ind_part,
+    size_t *Npart, const size_t offset_gparts, const int *restrict ind_gpart,
+    size_t *Ngpart, const size_t offset_sparts, const int *restrict ind_spart,
+    size_t *Nspart, const size_t offset_bparts, const int *restrict ind_bpart,
+    size_t *Nbpart, const size_t offset_sinks, const int *restrict ind_sink,
+    size_t *Nsink) {
 
 #ifdef WITH_MPI
   struct space *s = e->s;
@@ -80,6 +91,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
     e->proxies[k].nr_gparts_out = 0;
     e->proxies[k].nr_sparts_out = 0;
     e->proxies[k].nr_bparts_out = 0;
+    e->proxies[k].nr_sinks_out = 0;
   }
 
   /* Put the parts into the corresponding proxies. */
@@ -204,6 +216,46 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
     /* Load the bpart into the proxy */
     proxy_bparts_load(&e->proxies[pid], &s->bparts[offset_bparts + k], 1);
 
+#ifdef WITH_CSDS
+    if (e->policy & engine_policy_csds) {
+      error("Not yet implemented.");
+    }
+#endif
+  }
+
+  /* Put the sinks into the corresponding proxies. */
+  for (size_t k = 0; k < *Nsink; k++) {
+    /* Ignore the particles we want to get rid of (inhibited, ...). */
+    if (ind_sink[k] == -1) continue;
+
+    /* Get the target node and proxy ID. */
+    const int node_id = e->s->cells_top[ind_sink[k]].nodeID;
+    if (node_id < 0 || node_id >= e->nr_nodes)
+      error("Bad node ID %i.", node_id);
+    const int pid = e->proxy_ind[node_id];
+    if (pid < 0) {
+      error(
+          "Do not have a proxy for the requested nodeID %i for part with "
+          "id=%lld, x=[%e,%e,%e].",
+          node_id, s->sinks[offset_sinks + k].id,
+          s->sinks[offset_sinks + k].x[0], s->sinks[offset_sinks + k].x[1],
+          s->sinks[offset_sinks + k].x[2]);
+    }
+
+    /* Re-link the associated gpart with the buffer offset of the sink. */
+    if (s->sinks[offset_sinks + k].gpart != NULL) {
+      s->sinks[offset_sinks + k].gpart->id_or_neg_offset =
+          -e->proxies[pid].nr_sinks_out;
+    }
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (s->sinks[offset_sinks + k].time_bin == time_bin_inhibited)
+      error("Attempting to exchange an inhibited particle");
+#endif
+
+    /* Load the sink into the proxy */
+    proxy_sinks_load(&e->proxies[pid], &s->sinks[offset_sinks + k], 1);
+
 #ifdef WITH_CSDS
     if (e->policy & engine_policy_csds) {
       error("Not yet implemented.");
@@ -252,8 +304,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
   }
 
   /* Launch the proxies. */
-  MPI_Request reqs_in[5 * engine_maxproxies];
-  MPI_Request reqs_out[5 * engine_maxproxies];
+  MPI_Request reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
+  MPI_Request reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
   for (int k = 0; k < e->nr_proxies; k++) {
     proxy_parts_exchange_first(&e->proxies[k]);
     reqs_in[k] = e->proxies[k].req_parts_count_in;
@@ -281,18 +333,21 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
   int count_gparts_in = 0;
   int count_sparts_in = 0;
   int count_bparts_in = 0;
+  int count_sinks_in = 0;
   for (int k = 0; k < e->nr_proxies; k++) {
     count_parts_in += e->proxies[k].nr_parts_in;
     count_gparts_in += e->proxies[k].nr_gparts_in;
     count_sparts_in += e->proxies[k].nr_sparts_in;
     count_bparts_in += e->proxies[k].nr_bparts_in;
+    count_sinks_in += e->proxies[k].nr_sinks_in;
   }
   if (e->verbose) {
     message(
-        "sent out %zu/%zu/%zu/%zu parts/gparts/sparts/bparts, got %i/%i/%i/%i "
+        "sent out %zu/%zu/%zu/%zu/%zu parts/gparts/sparts/bparts/sinks, got "
+        "%i/%i/%i/%i/%i "
         "back.",
-        *Npart, *Ngpart, *Nspart, *Nbpart, count_parts_in, count_gparts_in,
-        count_sparts_in, count_bparts_in);
+        *Npart, *Ngpart, *Nspart, *Nbpart, *Nsink, count_parts_in,
+        count_gparts_in, count_sparts_in, count_bparts_in, count_sinks_in);
   }
 
   /* Reallocate the particle arrays if necessary */
@@ -356,6 +411,24 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
     }
   }
 
+  if (offset_sinks + count_sinks_in > s->size_sinks) {
+    s->size_sinks = (offset_sinks + count_sinks_in) * engine_parts_size_grow;
+    struct sink *sinks_new = NULL;
+    if (swift_memalign("sinks", (void **)&sinks_new, sink_align,
+                       sizeof(struct sink) * s->size_sinks) != 0)
+      error("Failed to allocate new sink data.");
+    memcpy(sinks_new, s->sinks, sizeof(struct sink) * offset_sinks);
+    swift_free("sinks", s->sinks);
+    s->sinks = sinks_new;
+
+    /* Reset the links */
+    for (size_t k = 0; k < offset_sinks; k++) {
+      if (s->sinks[k].gpart != NULL) {
+        s->sinks[k].gpart->id_or_neg_offset = -k;
+      }
+    }
+  }
+
   if (offset_gparts + count_gparts_in > s->size_gparts) {
     s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
     struct gpart *gparts_new = NULL;
@@ -374,6 +447,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
         s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
       } else if (s->gparts[k].type == swift_type_black_hole) {
         s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+      } else if (s->gparts[k].type == swift_type_sink) {
+        s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
       }
     }
   }
@@ -382,82 +457,113 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
   int nr_in = 0, nr_out = 0;
   for (int k = 0; k < e->nr_proxies; k++) {
     if (e->proxies[k].nr_parts_in > 0) {
-      reqs_in[5 * k] = e->proxies[k].req_parts_in;
-      reqs_in[5 * k + 1] = e->proxies[k].req_xparts_in;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+          e->proxies[k].req_parts_in;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+          e->proxies[k].req_xparts_in;
       nr_in += 2;
     } else {
-      reqs_in[5 * k] = reqs_in[5 * k + 1] = MPI_REQUEST_NULL;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+          reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_gparts_in > 0) {
-      reqs_in[5 * k + 2] = e->proxies[k].req_gparts_in;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
+          e->proxies[k].req_gparts_in;
       nr_in += 1;
     } else {
-      reqs_in[5 * k + 2] = MPI_REQUEST_NULL;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_sparts_in > 0) {
-      reqs_in[5 * k + 3] = e->proxies[k].req_sparts_in;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
+          e->proxies[k].req_sparts_in;
       nr_in += 1;
     } else {
-      reqs_in[5 * k + 3] = MPI_REQUEST_NULL;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_bparts_in > 0) {
-      reqs_in[5 * k + 4] = e->proxies[k].req_bparts_in;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
+          e->proxies[k].req_bparts_in;
+      nr_in += 1;
+    } else {
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
+    }
+    if (e->proxies[k].nr_sinks_in > 0) {
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
+          e->proxies[k].req_sinks_in;
       nr_in += 1;
     } else {
-      reqs_in[5 * k + 4] = MPI_REQUEST_NULL;
+      reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
     }
 
     if (e->proxies[k].nr_parts_out > 0) {
-      reqs_out[5 * k] = e->proxies[k].req_parts_out;
-      reqs_out[5 * k + 1] = e->proxies[k].req_xparts_out;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+          e->proxies[k].req_parts_out;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+          e->proxies[k].req_xparts_out;
       nr_out += 2;
     } else {
-      reqs_out[5 * k] = reqs_out[5 * k + 1] = MPI_REQUEST_NULL;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
+          reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
+              MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_gparts_out > 0) {
-      reqs_out[5 * k + 2] = e->proxies[k].req_gparts_out;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
+          e->proxies[k].req_gparts_out;
       nr_out += 1;
     } else {
-      reqs_out[5 * k + 2] = MPI_REQUEST_NULL;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_sparts_out > 0) {
-      reqs_out[5 * k + 3] = e->proxies[k].req_sparts_out;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
+          e->proxies[k].req_sparts_out;
       nr_out += 1;
     } else {
-      reqs_out[5 * k + 3] = MPI_REQUEST_NULL;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
     }
     if (e->proxies[k].nr_bparts_out > 0) {
-      reqs_out[5 * k + 4] = e->proxies[k].req_bparts_out;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
+          e->proxies[k].req_bparts_out;
       nr_out += 1;
     } else {
-      reqs_out[5 * k + 4] = MPI_REQUEST_NULL;
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
+    }
+    if (e->proxies[k].nr_sinks_out > 0) {
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
+          e->proxies[k].req_sinks_out;
+      nr_out += 1;
+    } else {
+      reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
     }
   }
 
   /* Wait for each part array to come in and collect the new
      parts from the proxies. */
-  int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0;
+  int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0,
+      count_sinks = 0;
   for (int k = 0; k < nr_in; k++) {
     int err, pid;
-    if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid,
-                           MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
+    if ((err = MPI_Waitany(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies,
+                           reqs_in, &pid, MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
       char buff[MPI_MAX_ERROR_STRING];
       int res;
       MPI_Error_string(err, buff, &res);
       error("MPI_Waitany failed (%s).", buff);
     }
     if (pid == MPI_UNDEFINED) break;
-    // message( "request from proxy %i has arrived." , pid / 5 );
-    pid = 5 * (pid / 5);
+    // message( "request from proxy %i has arrived." , pid /
+    // MPI_REQUEST_NUMBER_PARTICLE_TYPES );
+    pid = MPI_REQUEST_NUMBER_PARTICLE_TYPES *
+          (pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES);
 
     /* If all the requests for a given proxy have arrived... */
     if (reqs_in[pid + 0] == MPI_REQUEST_NULL &&
         reqs_in[pid + 1] == MPI_REQUEST_NULL &&
         reqs_in[pid + 2] == MPI_REQUEST_NULL &&
         reqs_in[pid + 3] == MPI_REQUEST_NULL &&
-        reqs_in[pid + 4] == MPI_REQUEST_NULL) {
+        reqs_in[pid + 4] == MPI_REQUEST_NULL &&
+        reqs_in[pid + 5] == MPI_REQUEST_NULL) {
       /* Copy the particle data to the part/xpart/gpart arrays. */
-      struct proxy *prox = &e->proxies[pid / 5];
+      struct proxy *prox = &e->proxies[pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES];
       memcpy(&s->parts[offset_parts + count_parts], prox->parts_in,
              sizeof(struct part) * prox->nr_parts_in);
       memcpy(&s->xparts[offset_parts + count_parts], prox->xparts_in,
@@ -468,6 +574,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
              sizeof(struct spart) * prox->nr_sparts_in);
       memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in,
              sizeof(struct bpart) * prox->nr_bparts_in);
+      memcpy(&s->sinks[offset_sinks + count_sinks], prox->sinks_in,
+             sizeof(struct sink) * prox->nr_sinks_in);
 
 #ifdef WITH_CSDS
       if (e->policy & engine_policy_csds) {
@@ -495,6 +603,12 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
         if (prox->nr_bparts_in > 0) {
           error("TODO");
         }
+
+        /* Log the sinks */
+        if (prox->nr_sinks_in > 0) {
+          /* Not implemented yet */
+          error("TODO");
+        }
       }
 #endif
       /* for (int k = offset; k < offset + count; k++)
@@ -522,6 +636,11 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
               &s->bparts[offset_bparts + count_bparts - gp->id_or_neg_offset];
           gp->id_or_neg_offset = s->bparts - bp;
           bp->gpart = gp;
+        } else if (gp->type == swift_type_sink) {
+          struct sink *sink =
+              &s->sinks[offset_sinks + count_sinks - gp->id_or_neg_offset];
+          gp->id_or_neg_offset = s->sinks - sink;
+          sink->gpart = gp;
         }
       }
 
@@ -530,13 +649,14 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
       count_gparts += prox->nr_gparts_in;
       count_sparts += prox->nr_sparts_in;
       count_bparts += prox->nr_bparts_in;
+      count_sinks += prox->nr_sinks_in;
     }
   }
 
   /* Wait for all the sends to have finished too. */
   if (nr_out > 0)
-    if (MPI_Waitall(5 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) !=
-        MPI_SUCCESS)
+    if (MPI_Waitall(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies, reqs_out,
+                    MPI_STATUSES_IGNORE) != MPI_SUCCESS)
       error("MPI_Waitall on sends failed.");
 
   /* Free the proxy memory */
@@ -553,6 +673,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
   *Ngpart = count_gparts;
   *Nspart = count_sparts;
   *Nbpart = count_bparts;
+  *Nsink = count_sinks;
 
 #else
   error("SWIFT was not compiled with MPI support.");

src/engine_unskip.c

@@ -160,12 +160,6 @@ static void engine_do_unskip_black_holes(struct cell *c, struct engine *e) {
   /* Early abort (are we below the level where tasks are)? */
   if (!cell_get_flag(c, cell_flag_has_tasks)) return;
 
-  /* Ignore empty cells. */
-  if (c->black_holes.count == 0) return;
-
-  /* Skip inactive cells. */
-  if (!cell_is_active_black_holes(c, e)) return;
-
   /* Recurse */
   if (c->split) {
     for (int k = 0; k < 8; k++) {
@@ -428,12 +422,6 @@ void engine_unskip(struct engine *e) {
         memswap(&local_cells[k], &local_cells[num_active_cells], sizeof(int));
       num_active_cells += 1;
     }
-
-    /* Activate the top-level timestep exchange */
-#ifdef WITH_MPI
-    scheduler_activate_all_subtype(&e->sched, c->mpi.send, task_subtype_tend);
-    scheduler_activate_all_subtype(&e->sched, c->mpi.recv, task_subtype_tend);
-#endif
   }
 
   /* What kind of tasks do we have? */

src/entropy_floor/EAGLE/entropy_floor.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser 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 Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include <config.h>
+
+#ifdef HAVE_HDF5
+#include <hdf5.h>
+#endif
+
+/* Local includes */
+#include "cosmology.h"
+#include "entropy_floor.h"
+#include "hydro.h"
+#include "parser.h"
+
+/**
+ * @brief Compute the pressure from the entropy floor at a given density
+ *
+ * @param rho_phys The physical density (internal units).
+ * @param rho_com The comoving density (internal units).
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+float entropy_floor_gas_pressure(const float rho_phys, const float rho_com,
+                                 const struct cosmology *cosmo,
+                                 const struct entropy_floor_properties *props) {
+
+  /* Mean baryon density in co-moving internal units for over-density condition
+   * (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
+   * making the over-density condition a no-op) */
+  const float rho_crit_0 = cosmo->critical_density_0;
+  const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
+
+  /* Physical pressure */
+  float pressure = 0.f;
+
+  /* Are we in the regime of the Jeans equation of state? */
+  if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
+      (rho_phys >= props->Jeans_density_threshold)) {
+
+    const float pressure_Jeans =
+        props->Jeans_pressure_norm *
+        powf(rho_phys * props->Jeans_density_threshold_inv,
+             props->Jeans_gamma_effective);
+
+    pressure = max(pressure, pressure_Jeans);
+  }
+
+  /* Are we in the regime of the Cool equation of state? */
+  if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
+      (rho_phys >= props->Cool_density_threshold)) {
+
+    const float pressure_Cool =
+        props->Cool_pressure_norm *
+        powf(rho_phys * props->Cool_density_threshold_inv,
+             props->Cool_gamma_effective);
+
+    pressure = max(pressure, pressure_Cool);
+  }
+
+  return pressure;
+}
+
+/**
+ * @brief Compute the entropy floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * @param p The #part.
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+float entropy_floor(const struct part *p, const struct cosmology *cosmo,
+                    const struct entropy_floor_properties *props) {
+
+  /* Comoving density in internal units */
+  const float rho_com = hydro_get_comoving_density(p);
+
+  /* Physical density in internal units */
+  const float rho_phys = hydro_get_physical_density(p, cosmo);
+
+  const float pressure =
+      entropy_floor_gas_pressure(rho_phys, rho_com, cosmo, props);
+
+  /* Convert to an entropy.
+   * (Recall that the entropy is the same in co-moving and physical frames) */
+  return gas_entropy_from_pressure(rho_phys, pressure);
+}
+
+/**
+ * @brief Compute the temperature from the entropy floor at a given density
+ *
+ * This is the temperature exactly corresponding to the imposed EoS shape.
+ * It only matches the entropy returned by the entropy_floor() function
+ * for a neutral gas with primoridal abundance.
+ *
+ * @param rho_phys The physical density (internal units).
+ * @param rho_com The comoving density (internal units).
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+float entropy_floor_gas_temperature(
+    const float rho_phys, const float rho_com, const struct cosmology *cosmo,
+    const struct entropy_floor_properties *props) {
+
+  /* Mean baryon density in co-moving internal units for over-density condition
+   * (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
+   * making the over-density condition a no-op) */
+  const float rho_crit_0 = cosmo->critical_density_0;
+  const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
+
+  /* Physical */
+  float temperature = 0.f;
+
+  /* Are we in the regime of the Jeans equation of state? */
+  if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
+      (rho_phys >= props->Jeans_density_threshold)) {
+
+    const float jeans_slope = props->Jeans_gamma_effective - 1.f;
+
+    const float temperature_Jeans =
+        props->Jeans_temperature_norm *
+        pow(rho_phys * props->Jeans_density_threshold_inv, jeans_slope);
+
+    temperature = max(temperature, temperature_Jeans);
+  }
+
+  /* Are we in the regime of the Cool equation of state? */
+  if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
+      (rho_phys >= props->Cool_density_threshold)) {
+
+    const float cool_slope = props->Cool_gamma_effective - 1.f;
+
+    const float temperature_Cool =
+        props->Cool_temperature_norm *
+        pow(rho_phys * props->Cool_density_threshold_inv, cool_slope);
+
+    temperature = max(temperature, temperature_Cool);
+  }
+
+  return temperature;
+}
+
+/**
+ * @brief Compute the temperature from the entropy floor for a given #part
+ *
+ * Calculate the EoS temperature, the particle is not updated.
+ * This is the temperature exactly corresponding to the imposed EoS shape.
+ * It only matches the entropy returned by the entropy_floor() function
+ * for a neutral gas with primoridal abundance.
+ *
+ * @param p The #part.
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+float entropy_floor_temperature(const struct part *p,
+                                const struct cosmology *cosmo,
+                                const struct entropy_floor_properties *props) {
+
+  /* Comoving density in internal units */
+  const float rho_com = hydro_get_comoving_density(p);
+
+  /* Physical density in internal units */
+  const float rho_phys = hydro_get_physical_density(p, cosmo);
+
+  return entropy_floor_gas_temperature(rho_phys, rho_com, cosmo, props);
+}
+
+/**
+ * @brief Initialise the entropy floor by reading the parameters and converting
+ * to internal units.
+ *
+ * The input temperatures and number densities are converted to entropy and
+ * density assuming a neutral gas of primoridal abundance.
+ *
+ * @param params The YAML parameter file.
+ * @param us The system of units used internally.
+ * @param phys_const The physical constants.
+ * @param hydro_props The propoerties of the hydro scheme.
+ * @param props The entropy floor properties to fill.
+ */
+void entropy_floor_init(struct entropy_floor_properties *props,
+                        const struct phys_const *phys_const,
+                        const struct unit_system *us,
+                        const struct hydro_props *hydro_props,
+                        struct swift_params *params) {
+
+  /* Read the parameters in the units they are set */
+  props->Jeans_density_threshold_H_p_cm3 = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3");
+  props->Jeans_over_density_threshold = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Jeans_over_density_threshold");
+  props->Jeans_temperature_norm_K = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Jeans_temperature_norm_K");
+  props->Jeans_gamma_effective =
+      parser_get_param_float(params, "EAGLEEntropyFloor:Jeans_gamma_effective");
+
+  props->Cool_density_threshold_H_p_cm3 = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3");
+  props->Cool_over_density_threshold = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Cool_over_density_threshold");
+  props->Cool_temperature_norm_K = parser_get_param_float(
+      params, "EAGLEEntropyFloor:Cool_temperature_norm_K");
+  props->Cool_gamma_effective =
+      parser_get_param_float(params, "EAGLEEntropyFloor:Cool_gamma_effective");
+
+  /* Cross-check that the input makes sense */
+  if (props->Cool_density_threshold_H_p_cm3 >=
+      props->Jeans_density_threshold_H_p_cm3) {
+    error(
+        "Invalid values for the entrop floor density thresholds. The 'Jeans' "
+        "threshold (%e cm^-3) should be at a higher density than the 'Cool' "
+        "threshold (%e cm^-3)",
+        props->Jeans_density_threshold_H_p_cm3,
+        props->Cool_density_threshold_H_p_cm3);
+  }
+
+  /* Initial Hydrogen abundance (mass fraction) */
+  const double X_H = hydro_props->hydrogen_mass_fraction;
+
+  /* Now convert to internal units assuming primodial Hydrogen abundance */
+  props->Jeans_temperature_norm =
+      props->Jeans_temperature_norm_K /
+      units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+  props->Jeans_density_threshold =
+      props->Jeans_density_threshold_H_p_cm3 /
+      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
+      phys_const->const_proton_mass / X_H;
+
+  props->Cool_temperature_norm =
+      props->Cool_temperature_norm_K /
+      units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+  props->Cool_density_threshold =
+      props->Cool_density_threshold_H_p_cm3 /
+      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
+      phys_const->const_proton_mass / X_H;
+
+  /* We assume neutral gas */
+  const float mean_molecular_weight = hydro_props->mu_neutral;
+
+  /* Get the common terms */
+  props->Jeans_density_threshold_inv = 1.f / props->Jeans_density_threshold;
+  props->Cool_density_threshold_inv = 1.f / props->Cool_density_threshold;
+
+  /* P_norm = (k_B * T) / (m_p * mu) * rho_threshold */
+  props->Jeans_pressure_norm =
+      ((phys_const->const_boltzmann_k * props->Jeans_temperature_norm) /
+       (phys_const->const_proton_mass * mean_molecular_weight)) *
+      props->Jeans_density_threshold;
+
+  props->Cool_pressure_norm =
+      ((phys_const->const_boltzmann_k * props->Cool_temperature_norm) /
+       (phys_const->const_proton_mass * mean_molecular_weight)) *
+      props->Cool_density_threshold;
+}
+
+/**
+ * @brief Print the properties of the entropy floor to stdout.
+ *
+ * @param props The entropy floor properties.
+ */
+void entropy_floor_print(const struct entropy_floor_properties *props) {
+
+  message("Entropy floor is 'EAGLE' with:");
+  message("Jeans limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
+          props->Jeans_gamma_effective, props->Jeans_density_threshold,
+          props->Jeans_density_threshold_H_p_cm3,
+          props->Jeans_temperature_norm);
+  message(" Cool limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
+          props->Cool_gamma_effective, props->Cool_density_threshold,
+          props->Cool_density_threshold_H_p_cm3, props->Cool_temperature_norm);
+}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of entropy floor to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+void entropy_floor_write_flavour(hid_t h_grp) {
+
+  io_write_attribute_s(h_grp, "Entropy floor", "EAGLE");
+}
+#endif
+
+/**
+ * @brief Write an entropy floor struct to the given FILE as a stream of bytes.
+ *
+ * @param props the struct
+ * @param stream the file stream
+ */
+void entropy_floor_struct_dump(const struct entropy_floor_properties *props,
+                               FILE *stream) {
+
+  restart_write_blocks((void *)props, sizeof(struct entropy_floor_properties),
+                       1, stream, "entropy floor", "entropy floor properties");
+}
+
+/**
+ * @brief Restore a entropy floor struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the struct
+ * @param stream the file stream
+ */
+void entropy_floor_struct_restore(struct entropy_floor_properties *props,
+                                  FILE *stream) {
+
+  restart_read_blocks((void *)props, sizeof(struct entropy_floor_properties), 1,
+                      stream, NULL, "entropy floor properties");
+}

src/entropy_floor/EAGLE/entropy_floor.h

@@ -19,13 +19,19 @@
 #ifndef SWIFT_ENTROPY_FLOOR_EAGLE_H
 #define SWIFT_ENTROPY_FLOOR_EAGLE_H
 
-#include "adiabatic_index.h"
-#include "cosmology.h"
-#include "equation_of_state.h"
-#include "hydro.h"
-#include "hydro_properties.h"
-#include "parser.h"
-#include "units.h"
+/* Code config */
+#include <config.h>
+
+/* System include */
+#include <stdio.h>
+
+/* Pre-declarations */
+struct cosmology;
+struct part;
+struct phys_const;
+struct unit_system;
+struct hydro_props;
+struct swift_params;
 
 /**
  * @file src/entropy_floor/EAGLE/entropy_floor.h
@@ -86,303 +92,38 @@ struct entropy_floor_properties {
   float Cool_pressure_norm;
 };
 
-/**
- * @brief Compute the pressure from the entropy floor at a given density
- *
- * @param rho_phys The physical density (internal units).
- * @param rho_com The comoving density (internal units).
- * @param cosmo The cosmological model.
- * @param props The properties of the entropy floor.
- */
-static INLINE float entropy_floor_gas_pressure(
-    const float rho_phys, const float rho_com, const struct cosmology *cosmo,
-    const struct entropy_floor_properties *props) {
-
-  /* Mean baryon density in co-moving internal units for over-density condition
-   * (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
-   * making the over-density condition a no-op) */
-  const float rho_crit_0 = cosmo->critical_density_0;
-  const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
-
-  /* Physical pressure */
-  float pressure = 0.f;
-
-  /* Are we in the regime of the Jeans equation of state? */
-  if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
-      (rho_phys >= props->Jeans_density_threshold)) {
-
-    const float pressure_Jeans =
-        props->Jeans_pressure_norm *
-        powf(rho_phys * props->Jeans_density_threshold_inv,
-             props->Jeans_gamma_effective);
-
-    pressure = max(pressure, pressure_Jeans);
-  }
-
-  /* Are we in the regime of the Cool equation of state? */
-  if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
-      (rho_phys >= props->Cool_density_threshold)) {
-
-    const float pressure_Cool =
-        props->Cool_pressure_norm *
-        powf(rho_phys * props->Cool_density_threshold_inv,
-             props->Cool_gamma_effective);
-
-    pressure = max(pressure, pressure_Cool);
-  }
-
-  return pressure;
-}
-
-/**
- * @brief Compute the entropy floor of a given #part.
- *
- * Note that the particle is not updated!!
- *
- * @param p The #part.
- * @param cosmo The cosmological model.
- * @param props The properties of the entropy floor.
- */
-static INLINE float entropy_floor(
-    const struct part *p, const struct cosmology *cosmo,
-    const struct entropy_floor_properties *props) {
+float entropy_floor_gas_pressure(const float rho_phys, const float rho_com,
+                                 const struct cosmology *cosmo,
+                                 const struct entropy_floor_properties *props);
 
-  /* Comoving density in internal units */
-  const float rho_com = hydro_get_comoving_density(p);
+float entropy_floor(const struct part *p, const struct cosmology *cosmo,
+                    const struct entropy_floor_properties *props);
 
-  /* Physical density in internal units */
-  const float rho_phys = hydro_get_physical_density(p, cosmo);
-
-  const float pressure =
-      entropy_floor_gas_pressure(rho_phys, rho_com, cosmo, props);
-
-  /* Convert to an entropy.
-   * (Recall that the entropy is the same in co-moving and physical frames) */
-  return gas_entropy_from_pressure(rho_phys, pressure);
-}
-
-/**
- * @brief Compute the temperature from the entropy floor at a given density
- *
- * This is the temperature exactly corresponding to the imposed EoS shape.
- * It only matches the entropy returned by the entropy_floor() function
- * for a neutral gas with primoridal abundance.
- *
- * @param rho_phys The physical density (internal units).
- * @param rho_com The comoving density (internal units).
- * @param cosmo The cosmological model.
- * @param props The properties of the entropy floor.
- */
-static INLINE float entropy_floor_gas_temperature(
+float entropy_floor_gas_temperature(
     const float rho_phys, const float rho_com, const struct cosmology *cosmo,
-    const struct entropy_floor_properties *props) {
-
-  /* Mean baryon density in co-moving internal units for over-density condition
-   * (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
-   * making the over-density condition a no-op) */
-  const float rho_crit_0 = cosmo->critical_density_0;
-  const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
-
-  /* Physical */
-  float temperature = 0.f;
-
-  /* Are we in the regime of the Jeans equation of state? */
-  if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
-      (rho_phys >= props->Jeans_density_threshold)) {
-
-    const float jeans_slope = props->Jeans_gamma_effective - 1.f;
-
-    const float temperature_Jeans =
-        props->Jeans_temperature_norm *
-        pow(rho_phys * props->Jeans_density_threshold_inv, jeans_slope);
+    const struct entropy_floor_properties *props);
 
-    temperature = max(temperature, temperature_Jeans);
-  }
+float entropy_floor_temperature(const struct part *p,
+                                const struct cosmology *cosmo,
+                                const struct entropy_floor_properties *props);
 
-  /* Are we in the regime of the Cool equation of state? */
-  if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
-      (rho_phys >= props->Cool_density_threshold)) {
+void entropy_floor_init(struct entropy_floor_properties *props,
+                        const struct phys_const *phys_const,
+                        const struct unit_system *us,
+                        const struct hydro_props *hydro_props,
+                        struct swift_params *params);
 
-    const float cool_slope = props->Cool_gamma_effective - 1.f;
-
-    const float temperature_Cool =
-        props->Cool_temperature_norm *
-        pow(rho_phys * props->Cool_density_threshold_inv, cool_slope);
-
-    temperature = max(temperature, temperature_Cool);
-  }
-
-  return temperature;
-}
-
-/**
- * @brief Compute the temperature from the entropy floor for a given #part
- *
- * Calculate the EoS temperature, the particle is not updated.
- * This is the temperature exactly corresponding to the imposed EoS shape.
- * It only matches the entropy returned by the entropy_floor() function
- * for a neutral gas with primoridal abundance.
- *
- * @param p The #part.
- * @param cosmo The cosmological model.
- * @param props The properties of the entropy floor.
- */
-static INLINE float entropy_floor_temperature(
-    const struct part *p, const struct cosmology *cosmo,
-    const struct entropy_floor_properties *props) {
-
-  /* Comoving density in internal units */
-  const float rho_com = hydro_get_comoving_density(p);
-
-  /* Physical density in internal units */
-  const float rho_phys = hydro_get_physical_density(p, cosmo);
-
-  return entropy_floor_gas_temperature(rho_phys, rho_com, cosmo, props);
-}
-
-/**
- * @brief Initialise the entropy floor by reading the parameters and converting
- * to internal units.
- *
- * The input temperatures and number densities are converted to entropy and
- * density assuming a neutral gas of primoridal abundance.
- *
- * @param params The YAML parameter file.
- * @param us The system of units used internally.
- * @param phys_const The physical constants.
- * @param hydro_props The propoerties of the hydro scheme.
- * @param props The entropy floor properties to fill.
- */
-static INLINE void entropy_floor_init(struct entropy_floor_properties *props,
-                                      const struct phys_const *phys_const,
-                                      const struct unit_system *us,
-                                      const struct hydro_props *hydro_props,
-                                      struct swift_params *params) {
-
-  /* Read the parameters in the units they are set */
-  props->Jeans_density_threshold_H_p_cm3 = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3");
-  props->Jeans_over_density_threshold = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Jeans_over_density_threshold");
-  props->Jeans_temperature_norm_K = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Jeans_temperature_norm_K");
-  props->Jeans_gamma_effective =
-      parser_get_param_float(params, "EAGLEEntropyFloor:Jeans_gamma_effective");
-
-  props->Cool_density_threshold_H_p_cm3 = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3");
-  props->Cool_over_density_threshold = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Cool_over_density_threshold");
-  props->Cool_temperature_norm_K = parser_get_param_float(
-      params, "EAGLEEntropyFloor:Cool_temperature_norm_K");
-  props->Cool_gamma_effective =
-      parser_get_param_float(params, "EAGLEEntropyFloor:Cool_gamma_effective");
-
-  /* Cross-check that the input makes sense */
-  if (props->Cool_density_threshold_H_p_cm3 >=
-      props->Jeans_density_threshold_H_p_cm3) {
-    error(
-        "Invalid values for the entrop floor density thresholds. The 'Jeans' "
-        "threshold (%e cm^-3) should be at a higher density than the 'Cool' "
-        "threshold (%e cm^-3)",
-        props->Jeans_density_threshold_H_p_cm3,
-        props->Cool_density_threshold_H_p_cm3);
-  }
-
-  /* Initial Hydrogen abundance (mass fraction) */
-  const double X_H = hydro_props->hydrogen_mass_fraction;
-
-  /* Now convert to internal units assuming primodial Hydrogen abundance */
-  props->Jeans_temperature_norm =
-      props->Jeans_temperature_norm_K /
-      units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
-  props->Jeans_density_threshold =
-      props->Jeans_density_threshold_H_p_cm3 /
-      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
-      phys_const->const_proton_mass / X_H;
-
-  props->Cool_temperature_norm =
-      props->Cool_temperature_norm_K /
-      units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
-  props->Cool_density_threshold =
-      props->Cool_density_threshold_H_p_cm3 /
-      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
-      phys_const->const_proton_mass / X_H;
-
-  /* We assume neutral gas */
-  const float mean_molecular_weight = hydro_props->mu_neutral;
-
-  /* Get the common terms */
-  props->Jeans_density_threshold_inv = 1.f / props->Jeans_density_threshold;
-  props->Cool_density_threshold_inv = 1.f / props->Cool_density_threshold;
-
-  /* P_norm = (k_B * T) / (m_p * mu) * rho_threshold */
-  props->Jeans_pressure_norm =
-      ((phys_const->const_boltzmann_k * props->Jeans_temperature_norm) /
-       (phys_const->const_proton_mass * mean_molecular_weight)) *
-      props->Jeans_density_threshold;
-
-  props->Cool_pressure_norm =
-      ((phys_const->const_boltzmann_k * props->Cool_temperature_norm) /
-       (phys_const->const_proton_mass * mean_molecular_weight)) *
-      props->Cool_density_threshold;
-}
-
-/**
- * @brief Print the properties of the entropy floor to stdout.
- *
- * @param props The entropy floor properties.
- */
-static INLINE void entropy_floor_print(
-    const struct entropy_floor_properties *props) {
-
-  message("Entropy floor is 'EAGLE' with:");
-  message("Jeans limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
-          props->Jeans_gamma_effective, props->Jeans_density_threshold,
-          props->Jeans_density_threshold_H_p_cm3,
-          props->Jeans_temperature_norm);
-  message(" Cool limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
-          props->Cool_gamma_effective, props->Cool_density_threshold,
-          props->Cool_density_threshold_H_p_cm3, props->Cool_temperature_norm);
-}
+void entropy_floor_print(const struct entropy_floor_properties *props);
 
 #ifdef HAVE_HDF5
 
-/**
- * @brief Writes the current model of entropy floor to the file
- * @param h_grp The HDF5 group in which to write
- */
-INLINE static void entropy_floor_write_flavour(hid_t h_grp) {
-
-  io_write_attribute_s(h_grp, "Entropy floor", "EAGLE");
-}
+void entropy_floor_write_flavour(hid_t h_grp);
 #endif
 
-/**
- * @brief Write an entropy floor struct to the given FILE as a stream of bytes.
- *
- * @param props the struct
- * @param stream the file stream
- */
-static INLINE void entropy_floor_struct_dump(
-    const struct entropy_floor_properties *props, FILE *stream) {
-
-  restart_write_blocks((void *)props, sizeof(struct entropy_floor_properties),
-                       1, stream, "entropy floor", "entropy floor properties");
-}
-
-/**
- * @brief Restore a entropy floor struct from the given FILE as a stream of
- * bytes.
- *
- * @param props the struct
- * @param stream the file stream
- */
-static INLINE void entropy_floor_struct_restore(
-    struct entropy_floor_properties *props, FILE *stream) {
+void entropy_floor_struct_dump(const struct entropy_floor_properties *props,
+                               FILE *stream);
 
-  restart_read_blocks((void *)props, sizeof(struct entropy_floor_properties), 1,
-                      stream, NULL, "entropy floor properties");
-}
+void entropy_floor_struct_restore(struct entropy_floor_properties *props,
+                                  FILE *stream);
 
 #endif /* SWIFT_ENTROPY_FLOOR_EAGLE_H */

src/equation_of_state.h

@@ -38,6 +38,8 @@
 #include "./equation_of_state/planetary/equation_of_state.h"
 #elif defined(EOS_BAROTROPIC_GAS)
 #include "./equation_of_state/barotropic/equation_of_state.h"
+#elif defined(EOS_CONTINUOUS_BAROTROPIC_GAS)
+#include "./equation_of_state/continuous_barotropic/equation_of_state.h"
 #else
 #error "Invalid choice of equation of state"
 #endif

src/equation_of_state/barotropic/equation_of_state.h

@@ -56,7 +56,7 @@ struct eos_parameters {
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -76,7 +76,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_entropy(float density, float entropy) {
+gas_pressure_from_entropy(const float density, const float entropy) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -96,7 +96,7 @@ gas_pressure_from_entropy(float density, float entropy) {
  * @return The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_pressure(float density, float pressure) {
+gas_entropy_from_pressure(const float density, const float pressure) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -116,7 +116,7 @@ gas_entropy_from_pressure(float density, float pressure) {
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_entropy(float density, float entropy) {
+gas_soundspeed_from_entropy(const float density, const float entropy) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -135,7 +135,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -155,7 +155,7 @@ gas_entropy_from_internal_energy(float density, float u) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -175,7 +175,7 @@ gas_pressure_from_internal_energy(float density, float u) {
  * @return The internal energy \f$u\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -195,7 +195,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -214,7 +214,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
  * @param P The pressure \f$P\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_pressure(float density, float P) {
+gas_soundspeed_from_pressure(const float density, const float P) {
 
   const float density_frac = density * eos.inverse_core_density;
   const float density_factor = pow_gamma(density_frac);
@@ -241,7 +241,7 @@ INLINE static void eos_init(struct eos_parameters *e,
       parser_get_param_float(params, "EoS:barotropic_vacuum_sound_speed");
   e->vacuum_sound_speed2 = vacuum_sound_speed * vacuum_sound_speed;
   e->inverse_core_density =
-      1. / parser_get_param_float(params, "EoS:barotropic_core_density");
+      1.f / parser_get_param_float(params, "EoS:barotropic_core_density");
 }
 /**
  * @brief Print the equation of state

src/equation_of_state/continuous_barotropic/equation_of_state.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2025  Matthieu Schaller (schaller@strw.leidenuniv.nl).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser 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 Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CONTINUOUS_BAROTROPIC_GAS_EQUATION_OF_STATE_H
+#define SWIFT_CONTINUOUS_BAROTROPIC_GAS_EQUATION_OF_STATE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "inline.h"
+#include "physical_constants.h"
+
+extern struct eos_parameters eos;
+
+/**
+ * @brief The parameters of the equation of state for the gas.
+ *
+ * Barotropic equation of state for the Leiden jet comparison project.
+ */
+struct eos_parameters {
+
+  /*! Square of barotropic sound speed in vacuum */
+  float vacuum_sound_speed2;
+
+  /*! Inverse of the core density */
+  float inverse_core_density;
+
+  /*! External pressure */
+  float external_pressure;
+};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Since we are using a barotropic EoS, the entropy value is ignored.
+ * Computes \f$u = c_0^2 \frac{1 + \left(\frac{\rho}{\rho_c}\right)^\gamma
+ * }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$A\f$.
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_internal_energy_from_entropy(const float density, const float entropy) {
+  error("Undefined");
+  return -1.;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Since we are using a barotropic EoS, the entropy value is ignored.
+ * Computes \f$P = c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$A\f$.
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_pressure_from_entropy(const float density, const float entropy) {
+
+  const float parenthesis =
+      1 + pow_gamma_minus_one(density * eos.inverse_core_density);
+  const float P = eos.vacuum_sound_speed2 * density * parenthesis;
+  return fmaxf(P, eos.external_pressure);
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * Since we are using a barotropic EoS, the pressure value is ignored.
+ * Computes \f$A = \rho^{1-\gamma}c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_entropy_from_pressure(const float density, const float pressure) {
+  error("Undefined");
+  return -1.;
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Since we are using a barotropic EoS, the entropy is ignored.
+ * Computes \f$c = \sqrt{c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$A\f$.
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_soundspeed_from_entropy(const float density, const float entropy) {
+
+  const float P = gas_pressure_from_entropy(density, entropy);
+  return sqrtf(hydro_gamma * P / density);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Since we are using a barotropic EoS, the internal energy value is ignored.
+ * Computes \f$A = \rho^{1-\gamma}c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_entropy_from_internal_energy(const float density, const float u) {
+  error("Undefined");
+  return -1.;
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Since we are using a barotropic EoS, the internal energy value is ignored.
+ * Computes \f$P = c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_pressure_from_internal_energy(const float density, const float u) {
+
+  const float parenthesis =
+      1 + pow_gamma_minus_one(density * eos.inverse_core_density);
+  const float P = eos.vacuum_sound_speed2 * density * parenthesis;
+  return fmaxf(P, eos.external_pressure);
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * Since we are using a barotropic EoS, the pressure value is ignored.
+ * Computes \f$u = c_0^2 \frac{1 + \left(\frac{\rho}{\rho_c}\right)^\gamma
+ * }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The internal energy \f$u\f$.
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_internal_energy_from_pressure(const float density, const float pressure) {
+
+  error("Undefined");
+  return -1.;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Since we are using a barotropic EoS, the internal energy value is ignored.
+ * Computes \f$c = \sqrt{c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_soundspeed_from_internal_energy(const float density, const float u) {
+
+  const float P = gas_pressure_from_internal_energy(density, u);
+  return sqrtf(hydro_gamma * P / density);
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Since we are using a barotropic EoS, the pressure value is ignored.
+ * Computes \f$c = \sqrt{c_0^2 \left(1 +
+ * \left(\frac{\rho}{\rho_c}\right)^\gamma\right)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline, const)) INLINE static float
+gas_soundspeed_from_pressure(const float density, const float P) {
+
+  return sqrtf(hydro_gamma * P / density);
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * Read the vacuum sound speed and core density from the parameter file.
+ *
+ * @param e The #eos_parameters.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ */
+INLINE static void eos_init(struct eos_parameters *e,
+                            const struct phys_const *phys_const,
+                            const struct unit_system *us,
+                            struct swift_params *params) {
+
+  const float vacuum_sound_speed =
+      parser_get_param_float(params, "EoS:barotropic_vacuum_sound_speed");
+  e->vacuum_sound_speed2 = vacuum_sound_speed * vacuum_sound_speed;
+  e->inverse_core_density =
+      1.f / parser_get_param_float(params, "EoS:barotropic_core_density");
+  e->external_pressure =
+      parser_get_param_float(params, "EoS:external_pressure");
+}
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+INLINE static void eos_print(const struct eos_parameters *e) {
+
+  message(
+      "Equation of state: Continuous barotropic gas with vacuum sound speed "
+      "set to %f, core density set to %f. and external pressure set to %f.",
+      sqrtf(e->vacuum_sound_speed2), 1. / e->inverse_core_density,
+      e->external_pressure);
+
+  message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+INLINE static void eos_print_snapshot(hid_t h_grpsph,
+                                      const struct eos_parameters *e) {
+
+  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+  io_write_attribute_s(h_grpsph, "Equation of state",
+                       "Continuous barotropic gas");
+}
+#endif
+
+#endif /* SWIFT_CONTINUOUS_BAROTROPIC_GAS_EQUATION_OF_STATE_H */

src/equation_of_state/ideal_gas/equation_of_state.h

@@ -46,7 +46,7 @@ struct eos_parameters {};
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
 
   return entropy * pow_gamma_minus_one(density) *
          hydro_one_over_gamma_minus_one;
@@ -61,7 +61,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_entropy(float density, float entropy) {
+gas_pressure_from_entropy(const float density, const float entropy) {
 
   return entropy * pow_gamma(density);
 }
@@ -76,7 +76,7 @@ gas_pressure_from_entropy(float density, float entropy) {
  * @return The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_pressure(float density, float pressure) {
+gas_entropy_from_pressure(const float density, const float pressure) {
 
   return pressure * pow_minus_gamma(density);
 }
@@ -90,7 +90,7 @@ gas_entropy_from_pressure(float density, float pressure) {
  * @param entropy The entropy \f$A\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_entropy(float density, float entropy) {
+gas_soundspeed_from_entropy(const float density, const float entropy) {
 
   return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
 }
@@ -104,7 +104,7 @@ gas_soundspeed_from_entropy(float density, float entropy) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
 
   return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
 }
@@ -118,7 +118,7 @@ gas_entropy_from_internal_energy(float density, float u) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
 
   return hydro_gamma_minus_one * u * density;
 }
@@ -133,7 +133,7 @@ gas_pressure_from_internal_energy(float density, float u) {
  * @return The internal energy \f$u\f$.
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
 
   return hydro_one_over_gamma_minus_one * pressure / density;
 }
@@ -147,7 +147,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
 
   return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
 }
@@ -161,7 +161,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
  * @param P The pressure \f$P\f$
  */
 __attribute__((always_inline, const)) INLINE static float
-gas_soundspeed_from_pressure(float density, float P) {
+gas_soundspeed_from_pressure(const float density, const float P) {
 
   return sqrtf(hydro_gamma * P / density);
 }

src/equation_of_state/isothermal/equation_of_state.h

@@ -50,7 +50,7 @@ struct eos_parameters {
  * @param entropy The entropy \f$S\f$.
  */
 __attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
+gas_internal_energy_from_entropy(const float density, const float entropy) {
 
   return eos.isothermal_internal_energy;
 }
@@ -65,7 +65,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
  * @param entropy The entropy \f$S\f$.
  */
 __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
-    float density, float entropy) {
+    const float density, const float entropy) {
 
   return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
 }
@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
  * @return The entropy \f$A\f$.
  */
 __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
-    float density, float pressure) {
+    const float density, const float pressure) {
 
   return hydro_gamma_minus_one * eos.isothermal_internal_energy *
          pow_minus_gamma_minus_one(density);
@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
  * @param entropy The entropy \f$S\f$.
  */
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
-    float density, float entropy) {
+    const float density, const float entropy) {
 
   return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);
@@ -114,7 +114,7 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
+gas_entropy_from_internal_energy(const float density, const float u) {
 
   return hydro_gamma_minus_one * eos.isothermal_internal_energy *
          pow_minus_gamma_minus_one(density);
@@ -130,7 +130,7 @@ gas_entropy_from_internal_energy(float density, float u) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
+gas_pressure_from_internal_energy(const float density, const float u) {
 
   return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
 }
@@ -145,7 +145,7 @@ gas_pressure_from_internal_energy(float density, float u) {
  * @return The internal energy \f$u\f$ (which is constant).
  */
 __attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
+gas_internal_energy_from_pressure(const float density, const float pressure) {
   return eos.isothermal_internal_energy;
 }
 
@@ -160,7 +160,7 @@ gas_internal_energy_from_pressure(float density, float pressure) {
  * @param u The internal energy \f$u\f$
  */
 __attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
+gas_soundspeed_from_internal_energy(const float density, const float u) {
 
   return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);
@@ -177,7 +177,7 @@ gas_soundspeed_from_internal_energy(float density, float u) {
  * @param P The pressure \f$P\f$
  */
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
-    float density, float P) {
+    const float density, const float P) {
 
   return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
                hydro_gamma_minus_one);

src/equation_of_state/planetary/equation_of_state.h

@@ -102,6 +102,10 @@ enum eos_planetary_material_id {
   eos_planetary_id_Til_basalt =
       eos_planetary_type_Til * eos_planetary_type_factor + 3,
 
+  /*! Tillotson ice */
+  eos_planetary_id_Til_ice =
+      eos_planetary_type_Til * eos_planetary_type_factor + 4,
+
   /* Hubbard & MacFarlane (1980) Uranus/Neptune */
 
   /*! Hydrogen-helium atmosphere */
@@ -134,6 +138,23 @@ enum eos_planetary_material_id {
   eos_planetary_id_SS08_water =
       eos_planetary_type_SESAME * eos_planetary_type_factor + 3,
 
+  /*! AQUA (Haldemann et al. 2020) SESAME-like water */
+  eos_planetary_id_AQUA =
+      eos_planetary_type_SESAME * eos_planetary_type_factor + 4,
+
+  /*! CMS19 hydrogen (Chabrier et al. 2019) SESAME-like hydrogen */
+  eos_planetary_id_CMS19_H =
+      eos_planetary_type_SESAME * eos_planetary_type_factor + 5,
+
+  /*! CMS19 helium (Chabrier et al. 2019) SESAME-like helium */
+  eos_planetary_id_CMS19_He =
+      eos_planetary_type_SESAME * eos_planetary_type_factor + 6,
+
+  /*! CD21 hydrogen-helium (Chabrier & Debras 2021) SESAME-like H-He mixture
+     (Y=0.245) */
+  eos_planetary_id_CD21_HHe =
+      eos_planetary_type_SESAME * eos_planetary_type_factor + 7,
+
   /* ANEOS */
 
   /*! ANEOS forsterite (Stewart et al. 2019) -- in SESAME-style tables */
@@ -149,6 +170,10 @@ enum eos_planetary_material_id {
       eos_planetary_type_ANEOS * eos_planetary_type_factor + 2,
 };
 
+/* Base material ID for custom Tillotson EoS */
+#define eos_planetary_Til_custom_base_id \
+  (eos_planetary_type_Til * eos_planetary_type_factor + 90)
+
 /* Individual EOS function headers. */
 #include "hm80.h"
 #include "ideal_gas.h"
@@ -160,9 +185,11 @@ enum eos_planetary_material_id {
  */
 struct eos_parameters {
   struct idg_params idg_def;
-  struct Til_params Til_iron, Til_granite, Til_water, Til_basalt;
+  struct Til_params Til_iron, Til_granite, Til_water, Til_basalt, Til_ice,
+      Til_custom[10];
   struct HM80_params HM80_HHe, HM80_ice, HM80_rock;
-  struct SESAME_params SESAME_iron, SESAME_basalt, SESAME_water, SS08_water;
+  struct SESAME_params SESAME_iron, SESAME_basalt, SESAME_water, SS08_water,
+      AQUA, CMS19_H, CMS19_He, CD21_HHe;
   struct SESAME_params ANEOS_forsterite, ANEOS_iron, ANEOS_Fe85Si15;
   struct SESAME_params custom[10];
 };
@@ -223,8 +250,20 @@ gas_internal_energy_from_entropy(float density, float entropy,
                                                   &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_internal_energy_from_entropy(density, entropy,
+                                                  &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_internal_energy_from_entropy(density, entropy,
+                                                    &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -278,6 +317,42 @@ gas_internal_energy_from_entropy(float density, float entropy,
                                                      &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.AQUA.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_AQUA: 1");
+#endif
+          return SESAME_internal_energy_from_entropy(density, entropy,
+                                                     &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CMS19_H.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CMS19_H: 1");
+#endif
+          return SESAME_internal_energy_from_entropy(density, entropy,
+                                                     &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CMS19_He.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CMS19_He: 1");
+#endif
+          return SESAME_internal_energy_from_entropy(density, entropy,
+                                                     &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CD21_HHe.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CD21_HHe: 1");
+#endif
+          return SESAME_internal_energy_from_entropy(density, entropy,
+                                                     &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -372,8 +447,19 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
           return Til_pressure_from_entropy(density, entropy, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_pressure_from_entropy(density, entropy, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_pressure_from_entropy(density, entropy,
+                                             &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -424,6 +510,22 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
                                               &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_pressure_from_entropy(density, entropy, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_pressure_from_entropy(density, entropy, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_pressure_from_entropy(density, entropy, &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_pressure_from_entropy(density, entropy, &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -519,8 +621,19 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
           return Til_entropy_from_pressure(density, P, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_entropy_from_pressure(density, P, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_entropy_from_pressure(density, P,
+                                             &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -567,6 +680,22 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
           return SESAME_entropy_from_pressure(density, P, &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_entropy_from_pressure(density, P, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_entropy_from_pressure(density, P, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_entropy_from_pressure(density, P, &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_entropy_from_pressure(density, P, &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -659,8 +788,19 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
           return Til_soundspeed_from_entropy(density, entropy, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_soundspeed_from_entropy(density, entropy, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_soundspeed_from_entropy(density, entropy,
+                                               &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -711,6 +851,24 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
                                                 &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_soundspeed_from_entropy(density, entropy, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_soundspeed_from_entropy(density, entropy, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_soundspeed_from_entropy(density, entropy,
+                                                &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_soundspeed_from_entropy(density, entropy,
+                                                &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -805,8 +963,19 @@ gas_entropy_from_internal_energy(float density, float u,
           return Til_entropy_from_internal_energy(density, u, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_entropy_from_internal_energy(density, u, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_entropy_from_internal_energy(density, u,
+                                                    &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -857,6 +1026,22 @@ gas_entropy_from_internal_energy(float density, float u,
                                                      &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_entropy_from_internal_energy(density, u, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_entropy_from_internal_energy(density, u, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_entropy_from_internal_energy(density, u, &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_entropy_from_internal_energy(density, u, &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -978,11 +1163,26 @@ gas_pressure_from_internal_energy(float density, float u,
           return Til_pressure_from_internal_energy(density, u, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.Til_ice.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_Til_ice: 1");
+#endif
+          return Til_pressure_from_internal_energy(density, u, &eos.Til_ice);
+          break;
+
         default:
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_pressure_from_internal_energy(density, u,
+                                                     &eos.Til_custom[i_custom]);
+          } else {
 #ifdef SWIFT_DEBUG_CHECKS
-          error("Unknown material ID! mat_id = %d", mat_id);
+            error("Unknown material ID! mat_id = %d", mat_id);
 #endif
-          return -1.f;
+            return -1.f;
+          }
       };
       break;
 
@@ -1070,6 +1270,40 @@ gas_pressure_from_internal_energy(float density, float u,
                                                       &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.AQUA.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_AQUA: 1");
+#endif
+          return SESAME_pressure_from_internal_energy(density, u, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CMS19_H.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CMS19_H: 1");
+#endif
+          return SESAME_pressure_from_internal_energy(density, u, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CMS19_He.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CMS19_He: 1");
+#endif
+          return SESAME_pressure_from_internal_energy(density, u,
+                                                      &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+#ifdef SWIFT_DEBUG_CHECKS
+          if (eos.CD21_HHe.mat_id != mat_id)
+            error("EoS not enabled. Please set EoS:planetary_use_CD21_HHe: 1");
+#endif
+          return SESAME_pressure_from_internal_energy(density, u,
+                                                      &eos.CD21_HHe);
+          break;
+
         default:
 #ifdef SWIFT_DEBUG_CHECKS
           error("Unknown material ID! mat_id = %d", mat_id);
@@ -1200,8 +1434,19 @@ gas_internal_energy_from_pressure(float density, float P,
           return Til_internal_energy_from_pressure(density, P, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_internal_energy_from_pressure(density, P, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_internal_energy_from_pressure(density, P,
+                                                     &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -1252,6 +1497,24 @@ gas_internal_energy_from_pressure(float density, float P,
                                                       &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_internal_energy_from_pressure(density, P, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_internal_energy_from_pressure(density, P, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_internal_energy_from_pressure(density, P,
+                                                      &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_internal_energy_from_pressure(density, P,
+                                                      &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -1350,8 +1613,19 @@ gas_soundspeed_from_internal_energy(float density, float u,
                                                      &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_soundspeed_from_internal_energy(density, u, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_soundspeed_from_internal_energy(
+                density, u, &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -1405,6 +1679,25 @@ gas_soundspeed_from_internal_energy(float density, float u,
                                                         &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_soundspeed_from_internal_energy(density, u, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_soundspeed_from_internal_energy(density, u,
+                                                        &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_soundspeed_from_internal_energy(density, u,
+                                                        &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_soundspeed_from_internal_energy(density, u,
+                                                        &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -1499,8 +1792,19 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
           return Til_soundspeed_from_pressure(density, P, &eos.Til_basalt);
           break;
 
+        case eos_planetary_id_Til_ice:
+          return Til_soundspeed_from_pressure(density, P, &eos.Til_ice);
+          break;
+
         default:
-          return -1.f;
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_soundspeed_from_pressure(density, P,
+                                                &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
       };
       break;
 
@@ -1548,6 +1852,22 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
           return SESAME_soundspeed_from_pressure(density, P, &eos.SS08_water);
           break;
 
+        case eos_planetary_id_AQUA:
+          return SESAME_soundspeed_from_pressure(density, P, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_soundspeed_from_pressure(density, P, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_soundspeed_from_pressure(density, P, &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_soundspeed_from_pressure(density, P, &eos.CD21_HHe);
+          break;
+
         default:
           return -1.f;
       };
@@ -1591,106 +1911,740 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
 }
 
 /**
- * @brief Initialize the eos parameters
+ * @brief Returns the temperature given density and internal energy
  *
- * @param e The #eos_parameters
- * @param params The parsed parameters
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
  */
-__attribute__((always_inline)) INLINE static void eos_init(
-    struct eos_parameters *e, const struct phys_const *phys_const,
-    const struct unit_system *us, struct swift_params *params) {
+__attribute__((always_inline)) INLINE static float
+gas_temperature_from_internal_energy(float density, float u,
+                                     enum eos_planetary_material_id mat_id) {
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
-  // Prepare any/all requested EoS: Set the parameters and material IDs, load
-  // tables etc., and convert to internal units
+  /* Select the material base type */
+  switch (type) {
 
-  // Ideal gas
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_idg_def", 0)) {
-    set_idg_def(&e->idg_def, eos_planetary_id_idg_def);
-  }
+    /* Ideal gas EoS */
+    case eos_planetary_type_idg:
 
-  // Tillotson
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_iron", 0)) {
-    set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
-    convert_units_Til(&e->Til_iron, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_granite", 0)) {
-    set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
-    convert_units_Til(&e->Til_granite, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_water", 0)) {
-    set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
-    convert_units_Til(&e->Til_water, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_basalt", 0)) {
-    set_Til_basalt(&e->Til_basalt, eos_planetary_id_Til_basalt);
-    convert_units_Til(&e->Til_basalt, us);
-  }
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_idg_def:
+          return idg_temperature_from_internal_energy(density, u, &eos.idg_def);
+          break;
 
-  // Hubbard & MacFarlane (1980)
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_HHe", 0)) {
-    char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
-    set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
-    parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
-                            HM80_HHe_table_file);
-    load_table_HM80(&e->HM80_HHe, HM80_HHe_table_file);
-    prepare_table_HM80(&e->HM80_HHe);
-    convert_units_HM80(&e->HM80_HHe, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_ice", 0)) {
-    char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
-    set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
-    parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
-                            HM80_ice_table_file);
-    load_table_HM80(&e->HM80_ice, HM80_ice_table_file);
-    prepare_table_HM80(&e->HM80_ice);
-    convert_units_HM80(&e->HM80_ice, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_rock", 0)) {
-    char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
-    set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
-    parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
-                            HM80_rock_table_file);
-    load_table_HM80(&e->HM80_rock, HM80_rock_table_file);
-    prepare_table_HM80(&e->HM80_rock);
-    convert_units_HM80(&e->HM80_rock, us);
-  }
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
 
-  // SESAME
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_iron", 0)) {
-    char SESAME_iron_table_file[PARSER_MAX_LINE_SIZE];
-    set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
-    parser_get_param_string(params, "EoS:planetary_SESAME_iron_table_file",
-                            SESAME_iron_table_file);
-    load_table_SESAME(&e->SESAME_iron, SESAME_iron_table_file);
-    prepare_table_SESAME(&e->SESAME_iron);
-    convert_units_SESAME(&e->SESAME_iron, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_basalt", 0)) {
-    char SESAME_basalt_table_file[PARSER_MAX_LINE_SIZE];
-    set_SESAME_basalt(&e->SESAME_basalt, eos_planetary_id_SESAME_basalt);
-    parser_get_param_string(params, "EoS:planetary_SESAME_basalt_table_file",
-                            SESAME_basalt_table_file);
-    load_table_SESAME(&e->SESAME_basalt, SESAME_basalt_table_file);
-    prepare_table_SESAME(&e->SESAME_basalt);
-    convert_units_SESAME(&e->SESAME_basalt, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_water", 0)) {
-    char SESAME_water_table_file[PARSER_MAX_LINE_SIZE];
-    set_SESAME_water(&e->SESAME_water, eos_planetary_id_SESAME_water);
-    parser_get_param_string(params, "EoS:planetary_SESAME_water_table_file",
-                            SESAME_water_table_file);
-    load_table_SESAME(&e->SESAME_water, SESAME_water_table_file);
-    prepare_table_SESAME(&e->SESAME_water);
-    convert_units_SESAME(&e->SESAME_water, us);
-  }
-  if (parser_get_opt_param_int(params, "EoS:planetary_use_SS08_water", 0)) {
-    char SS08_water_table_file[PARSER_MAX_LINE_SIZE];
-    set_SS08_water(&e->SESAME_water, eos_planetary_id_SS08_water);
-    parser_get_param_string(params, "EoS:planetary_SS08_water_table_file",
-                            SS08_water_table_file);
-    load_table_SESAME(&e->SS08_water, SS08_water_table_file);
-    prepare_table_SESAME(&e->SS08_water);
-    convert_units_SESAME(&e->SS08_water, us);
+    /* Tillotson EoS */
+    case eos_planetary_type_Til:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_Til_iron:
+          return Til_temperature_from_internal_energy(density, u,
+                                                      &eos.Til_iron);
+          break;
+
+        case eos_planetary_id_Til_granite:
+          return Til_temperature_from_internal_energy(density, u,
+                                                      &eos.Til_granite);
+          break;
+
+        case eos_planetary_id_Til_water:
+          return Til_temperature_from_internal_energy(density, u,
+                                                      &eos.Til_water);
+          break;
+
+        case eos_planetary_id_Til_basalt:
+          return Til_temperature_from_internal_energy(density, u,
+                                                      &eos.Til_basalt);
+          break;
+
+        case eos_planetary_id_Til_ice:
+          return Til_temperature_from_internal_energy(density, u, &eos.Til_ice);
+          break;
+
+        default:
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_temperature_from_internal_energy(
+                density, u, &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
+      };
+      break;
+
+    /* Hubbard & MacFarlane (1980) EoS */
+    case eos_planetary_type_HM80:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_HM80_HHe:
+          return HM80_temperature_from_internal_energy(density, u,
+                                                       &eos.HM80_HHe);
+          break;
+
+        case eos_planetary_id_HM80_ice:
+          return HM80_temperature_from_internal_energy(density, u,
+                                                       &eos.HM80_ice);
+          break;
+
+        case eos_planetary_id_HM80_rock:
+          return HM80_temperature_from_internal_energy(density, u,
+                                                       &eos.HM80_rock);
+          break;
+
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
+
+    /* SESAME EoS */
+    case eos_planetary_type_SESAME:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_SESAME_iron:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.SESAME_iron);
+          break;
+
+        case eos_planetary_id_SESAME_basalt:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.SESAME_basalt);
+          break;
+
+        case eos_planetary_id_SESAME_water:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.SESAME_water);
+          break;
+
+        case eos_planetary_id_SS08_water:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.SS08_water);
+          break;
+
+        case eos_planetary_id_AQUA:
+          return SESAME_temperature_from_internal_energy(density, u, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.CD21_HHe);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /* ANEOS -- using SESAME-style tables */
+    case eos_planetary_type_ANEOS:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_ANEOS_forsterite:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.ANEOS_forsterite);
+          break;
+
+        case eos_planetary_id_ANEOS_iron:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.ANEOS_iron);
+          break;
+
+        case eos_planetary_id_ANEOS_Fe85Si15:
+          return SESAME_temperature_from_internal_energy(density, u,
+                                                         &eos.ANEOS_Fe85Si15);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /*! Generic user-provided custom tables */
+    case eos_planetary_type_custom: {
+      const int i_custom =
+          mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+      return SESAME_temperature_from_internal_energy(density, u,
+                                                     &eos.custom[i_custom]);
+      break;
+    }
+
+    default:
+      return -1.f;
+  }
+}
+
+/**
+ * @brief Returns the density given pressure and temperature
+ *
+ * @param P The pressure \f$P\f$
+ * @param T The temperature \f$T\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_density_from_pressure_and_temperature(
+    float P, float T, enum eos_planetary_material_id mat_id) {
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
+
+  /* Select the material base type */
+  switch (type) {
+
+    /* Ideal gas EoS */
+    case eos_planetary_type_idg:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_idg_def:
+          return idg_density_from_pressure_and_temperature(P, T, &eos.idg_def);
+          break;
+
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
+
+    /* Tillotson EoS */
+    case eos_planetary_type_Til:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_Til_iron:
+          return Til_density_from_pressure_and_temperature(P, T, &eos.Til_iron);
+          break;
+
+        case eos_planetary_id_Til_granite:
+          return Til_density_from_pressure_and_temperature(P, T,
+                                                           &eos.Til_granite);
+          break;
+
+        case eos_planetary_id_Til_water:
+          return Til_density_from_pressure_and_temperature(P, T,
+                                                           &eos.Til_water);
+          break;
+
+        case eos_planetary_id_Til_basalt:
+          return Til_density_from_pressure_and_temperature(P, T,
+                                                           &eos.Til_basalt);
+          break;
+
+        case eos_planetary_id_Til_ice:
+          return Til_density_from_pressure_and_temperature(P, T, &eos.Til_ice);
+          break;
+
+        default:
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_density_from_pressure_and_temperature(
+                P, T, &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
+      };
+      break;
+
+    /* Hubbard & MacFarlane (1980) EoS */
+    case eos_planetary_type_HM80:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_HM80_HHe:
+          return HM80_density_from_pressure_and_temperature(P, T,
+                                                            &eos.HM80_HHe);
+          break;
+
+        case eos_planetary_id_HM80_ice:
+          return HM80_density_from_pressure_and_temperature(P, T,
+                                                            &eos.HM80_ice);
+          break;
+
+        case eos_planetary_id_HM80_rock:
+          return HM80_density_from_pressure_and_temperature(P, T,
+                                                            &eos.HM80_rock);
+          break;
+
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
+
+    /* SESAME EoS */
+    case eos_planetary_type_SESAME:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_SESAME_iron:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.SESAME_iron);
+          break;
+
+        case eos_planetary_id_SESAME_basalt:
+          return SESAME_density_from_pressure_and_temperature(
+              P, T, &eos.SESAME_basalt);
+          break;
+
+        case eos_planetary_id_SESAME_water:
+          return SESAME_density_from_pressure_and_temperature(
+              P, T, &eos.SESAME_water);
+          break;
+
+        case eos_planetary_id_SS08_water:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.SS08_water);
+          break;
+
+        case eos_planetary_id_AQUA:
+          return SESAME_density_from_pressure_and_temperature(P, T, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.CD21_HHe);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /* ANEOS -- using SESAME-style tables */
+    case eos_planetary_type_ANEOS:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_ANEOS_forsterite:
+          return SESAME_density_from_pressure_and_temperature(
+              P, T, &eos.ANEOS_forsterite);
+          break;
+
+        case eos_planetary_id_ANEOS_iron:
+          return SESAME_density_from_pressure_and_temperature(P, T,
+                                                              &eos.ANEOS_iron);
+          break;
+
+        case eos_planetary_id_ANEOS_Fe85Si15:
+          return SESAME_density_from_pressure_and_temperature(
+              P, T, &eos.ANEOS_Fe85Si15);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /*! Generic user-provided custom tables */
+    case eos_planetary_type_custom: {
+      const int i_custom =
+          mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+      return SESAME_density_from_pressure_and_temperature(
+          P, T, &eos.custom[i_custom]);
+      break;
+    }
+
+    default:
+      return -1.f;
+  }
+}
+
+/**
+ * @brief Returns the density given pressure and internal energy
+ *
+ * @param P The pressure \f$P\f$
+ * @param T The temperature \f$T\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_density_from_pressure_and_internal_energy(
+    float P, float u, float rho_ref, float rho_sph,
+    enum eos_planetary_material_id mat_id) {
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
+
+  /* Select the material base type */
+  switch (type) {
+
+    /* Ideal gas EoS */
+    case eos_planetary_type_idg:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_idg_def:
+          return idg_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.idg_def);
+          break;
+
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
+
+    /* Tillotson EoS */
+    case eos_planetary_type_Til:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_Til_iron:
+          return Til_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.Til_iron);
+          break;
+
+        case eos_planetary_id_Til_granite:
+          return Til_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.Til_granite);
+          break;
+
+        case eos_planetary_id_Til_water:
+          return Til_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.Til_water);
+          break;
+
+        case eos_planetary_id_Til_basalt:
+          return Til_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.Til_basalt);
+          break;
+
+        case eos_planetary_id_Til_ice:
+          return Til_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.Til_ice);
+          break;
+
+        default:
+          // Custom user-provided Tillotson
+          if (mat_id >= eos_planetary_Til_custom_base_id) {
+            const int i_custom = mat_id - eos_planetary_Til_custom_base_id;
+            return Til_density_from_pressure_and_internal_energy(
+                P, u, rho_ref, rho_sph, &eos.Til_custom[i_custom]);
+          } else {
+            return -1.f;
+          }
+      };
+      break;
+
+    /* Hubbard & MacFarlane (1980) EoS */
+    case eos_planetary_type_HM80:
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_HM80_HHe:
+          return HM80_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.HM80_HHe);
+          break;
+
+        case eos_planetary_id_HM80_ice:
+          return HM80_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.HM80_ice);
+          break;
+
+        case eos_planetary_id_HM80_rock:
+          return HM80_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.HM80_rock);
+          break;
+
+        default:
+#ifdef SWIFT_DEBUG_CHECKS
+          error("Unknown material ID! mat_id = %d", mat_id);
+#endif
+          return -1.f;
+      };
+      break;
+
+    /* SESAME EoS */
+    case eos_planetary_type_SESAME:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_SESAME_iron:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.SESAME_iron);
+          break;
+
+        case eos_planetary_id_SESAME_basalt:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.SESAME_basalt);
+          break;
+
+        case eos_planetary_id_SESAME_water:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.SESAME_water);
+          break;
+
+        case eos_planetary_id_SS08_water:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.SS08_water);
+          break;
+
+        case eos_planetary_id_AQUA:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.AQUA);
+          break;
+
+        case eos_planetary_id_CMS19_H:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.CMS19_H);
+          break;
+
+        case eos_planetary_id_CMS19_He:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.CMS19_He);
+          break;
+
+        case eos_planetary_id_CD21_HHe:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.CD21_HHe);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /* ANEOS -- using SESAME-style tables */
+    case eos_planetary_type_ANEOS:;
+
+      /* Select the material of this type */
+      switch (mat_id) {
+        case eos_planetary_id_ANEOS_forsterite:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.ANEOS_forsterite);
+          break;
+
+        case eos_planetary_id_ANEOS_iron:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.ANEOS_iron);
+          break;
+
+        case eos_planetary_id_ANEOS_Fe85Si15:
+          return SESAME_density_from_pressure_and_internal_energy(
+              P, u, rho_ref, rho_sph, &eos.ANEOS_Fe85Si15);
+          break;
+
+        default:
+          return -1.f;
+      };
+      break;
+
+    /*! Generic user-provided custom tables */
+    case eos_planetary_type_custom: {
+      const int i_custom =
+          mat_id - eos_planetary_type_custom * eos_planetary_type_factor;
+      return SESAME_density_from_pressure_and_internal_energy(
+          P, u, rho_ref, rho_sph, &eos.custom[i_custom]);
+      break;
+    }
+
+    default:
+      return -1.f;
+  }
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_parameters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+    struct eos_parameters *e, const struct phys_const *phys_const,
+    const struct unit_system *us, struct swift_params *params) {
+
+  // Prepare any/all requested EoS: Set the parameters and material IDs, load
+  // tables etc., and convert to internal units
+
+  // Ideal gas
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_idg_def", 0)) {
+    set_idg_def(&e->idg_def, eos_planetary_id_idg_def);
+  }
+
+  // Tillotson
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_iron", 0)) {
+    set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
+    set_Til_u_cold(&e->Til_iron, eos_planetary_id_Til_iron);
+    convert_units_Til(&e->Til_iron, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_granite", 0)) {
+    set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
+    set_Til_u_cold(&e->Til_granite, eos_planetary_id_Til_granite);
+    convert_units_Til(&e->Til_granite, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_water", 0)) {
+    set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+    set_Til_u_cold(&e->Til_water, eos_planetary_id_Til_water);
+    convert_units_Til(&e->Til_water, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_basalt", 0)) {
+    set_Til_basalt(&e->Til_basalt, eos_planetary_id_Til_basalt);
+    set_Til_u_cold(&e->Til_basalt, eos_planetary_id_Til_basalt);
+    convert_units_Til(&e->Til_basalt, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_Til_ice", 0)) {
+    set_Til_ice(&e->Til_ice, eos_planetary_id_Til_ice);
+    set_Til_u_cold(&e->Til_ice, eos_planetary_id_Til_ice);
+    convert_units_Til(&e->Til_ice, us);
+  }
+
+  // Custom user-provided Tillotson
+  for (int i_custom = 0; i_custom <= 9; i_custom++) {
+    char param_name[PARSER_MAX_LINE_SIZE];
+    sprintf(param_name, "EoS:planetary_use_Til_custom_%d", i_custom);
+    if (parser_get_opt_param_int(params, param_name, 0)) {
+      char Til_custom_file[PARSER_MAX_LINE_SIZE];
+      int mat_id = eos_planetary_Til_custom_base_id + i_custom;
+
+      sprintf(param_name, "EoS:planetary_Til_custom_%d_param_file", i_custom);
+      parser_get_param_string(params, param_name, Til_custom_file);
+
+      set_Til_custom(&e->Til_custom[i_custom],
+                     (enum eos_planetary_material_id)mat_id, Til_custom_file);
+      set_Til_u_cold(&e->Til_custom[i_custom],
+                     (enum eos_planetary_material_id)mat_id);
+      convert_units_Til(&e->Til_custom[i_custom], us);
+    }
+  }
+
+  // Hubbard & MacFarlane (1980)
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_HHe", 0)) {
+    char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
+    set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
+    parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
+                            HM80_HHe_table_file);
+    load_table_HM80(&e->HM80_HHe, HM80_HHe_table_file);
+    prepare_table_HM80(&e->HM80_HHe);
+    convert_units_HM80(&e->HM80_HHe, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_ice", 0)) {
+    char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
+    set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
+    parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
+                            HM80_ice_table_file);
+    load_table_HM80(&e->HM80_ice, HM80_ice_table_file);
+    prepare_table_HM80(&e->HM80_ice);
+    convert_units_HM80(&e->HM80_ice, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80_rock", 0)) {
+    char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
+    set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
+    parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
+                            HM80_rock_table_file);
+    load_table_HM80(&e->HM80_rock, HM80_rock_table_file);
+    prepare_table_HM80(&e->HM80_rock);
+    convert_units_HM80(&e->HM80_rock, us);
+  }
+
+  // SESAME
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_iron", 0)) {
+    char SESAME_iron_table_file[PARSER_MAX_LINE_SIZE];
+    set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+    parser_get_param_string(params, "EoS:planetary_SESAME_iron_table_file",
+                            SESAME_iron_table_file);
+    load_table_SESAME(&e->SESAME_iron, SESAME_iron_table_file);
+    prepare_table_SESAME(&e->SESAME_iron);
+    convert_units_SESAME(&e->SESAME_iron, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_basalt", 0)) {
+    char SESAME_basalt_table_file[PARSER_MAX_LINE_SIZE];
+    set_SESAME_basalt(&e->SESAME_basalt, eos_planetary_id_SESAME_basalt);
+    parser_get_param_string(params, "EoS:planetary_SESAME_basalt_table_file",
+                            SESAME_basalt_table_file);
+    load_table_SESAME(&e->SESAME_basalt, SESAME_basalt_table_file);
+    prepare_table_SESAME(&e->SESAME_basalt);
+    convert_units_SESAME(&e->SESAME_basalt, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME_water", 0)) {
+    char SESAME_water_table_file[PARSER_MAX_LINE_SIZE];
+    set_SESAME_water(&e->SESAME_water, eos_planetary_id_SESAME_water);
+    parser_get_param_string(params, "EoS:planetary_SESAME_water_table_file",
+                            SESAME_water_table_file);
+    load_table_SESAME(&e->SESAME_water, SESAME_water_table_file);
+    prepare_table_SESAME(&e->SESAME_water);
+    convert_units_SESAME(&e->SESAME_water, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_SS08_water", 0)) {
+    char SS08_water_table_file[PARSER_MAX_LINE_SIZE];
+    set_SS08_water(&e->SS08_water, eos_planetary_id_SS08_water);
+    parser_get_param_string(params, "EoS:planetary_SS08_water_table_file",
+                            SS08_water_table_file);
+    load_table_SESAME(&e->SS08_water, SS08_water_table_file);
+    prepare_table_SESAME(&e->SS08_water);
+    convert_units_SESAME(&e->SS08_water, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_AQUA", 0)) {
+    char AQUA_table_file[PARSER_MAX_LINE_SIZE];
+    set_AQUA(&e->AQUA, eos_planetary_id_AQUA);
+    parser_get_param_string(params, "EoS:planetary_AQUA_table_file",
+                            AQUA_table_file);
+    load_table_SESAME(&e->AQUA, AQUA_table_file);
+    prepare_table_SESAME(&e->AQUA);
+    convert_units_SESAME(&e->AQUA, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_CMS19_H", 0)) {
+    char CMS19_H_table_file[PARSER_MAX_LINE_SIZE];
+    set_CMS19_H(&e->CMS19_H, eos_planetary_id_CMS19_H);
+    parser_get_param_string(params, "EoS:planetary_CMS19_H_table_file",
+                            CMS19_H_table_file);
+    load_table_SESAME(&e->CMS19_H, CMS19_H_table_file);
+    prepare_table_SESAME(&e->CMS19_H);
+    convert_units_SESAME(&e->CMS19_H, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_CMS19_He", 0)) {
+    char CMS19_He_table_file[PARSER_MAX_LINE_SIZE];
+    set_CMS19_He(&e->CMS19_He, eos_planetary_id_CMS19_He);
+    parser_get_param_string(params, "EoS:planetary_CMS19_He_table_file",
+                            CMS19_He_table_file);
+    load_table_SESAME(&e->CMS19_He, CMS19_He_table_file);
+    prepare_table_SESAME(&e->CMS19_He);
+    convert_units_SESAME(&e->CMS19_He, us);
+  }
+  if (parser_get_opt_param_int(params, "EoS:planetary_use_CD21_HHe", 0)) {
+    char CD21_HHe_table_file[PARSER_MAX_LINE_SIZE];
+    set_CD21_HHe(&e->CD21_HHe, eos_planetary_id_CD21_HHe);
+    parser_get_param_string(params, "EoS:planetary_CD21_HHe_table_file",
+                            CD21_HHe_table_file);
+    load_table_SESAME(&e->CD21_HHe, CD21_HHe_table_file);
+    prepare_table_SESAME(&e->CD21_HHe);
+    convert_units_SESAME(&e->CD21_HHe, us);
   }
 
   // ANEOS -- using SESAME-style tables

src/equation_of_state/planetary/get_eos_tables.sh

 #!/bin/bash
 
 # Download the tables of the publicly available planetary equations of state
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_HHe.txt
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_ice.txt
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_rock.txt
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_HHe.txt
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_ice.txt
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/planetary_HM80_rock.txt
 
 mv planetary_HM80_HHe.txt ../../../examples/
 mv planetary_HM80_ice.txt ../../../examples/

src/equation_of_state/planetary/hm80.h

@@ -179,7 +179,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat,
 
 // gas_internal_energy_from_entropy
 INLINE static float HM80_internal_energy_from_entropy(
-    float density, float entropy, const struct HM80_params *mat) {
+    const float density, const float entropy, const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -187,7 +187,8 @@ INLINE static float HM80_internal_energy_from_entropy(
 }
 
 // gas_pressure_from_entropy
-INLINE static float HM80_pressure_from_entropy(float density, float entropy,
+INLINE static float HM80_pressure_from_entropy(const float density,
+                                               const float entropy,
                                                const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -196,7 +197,8 @@ INLINE static float HM80_pressure_from_entropy(float density, float entropy,
 }
 
 // gas_entropy_from_pressure
-INLINE static float HM80_entropy_from_pressure(float density, float pressure,
+INLINE static float HM80_entropy_from_pressure(const float density,
+                                               const float pressure,
                                                const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -206,7 +208,7 @@ INLINE static float HM80_entropy_from_pressure(float density, float pressure,
 
 // gas_soundspeed_from_entropy
 INLINE static float HM80_soundspeed_from_entropy(
-    float density, float entropy, const struct HM80_params *mat) {
+    const float density, const float entropy, const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -215,29 +217,25 @@ INLINE static float HM80_soundspeed_from_entropy(
 
 // gas_entropy_from_internal_energy
 INLINE static float HM80_entropy_from_internal_energy(
-    float density, float u, const struct HM80_params *mat) {
+    const float density, const float u, const struct HM80_params *mat) {
 
   return 0.f;
 }
 
 // gas_pressure_from_internal_energy
 INLINE static float HM80_pressure_from_internal_energy(
-    float density, float u, const struct HM80_params *mat) {
-
-  float log_P, log_P_1, log_P_2, log_P_3, log_P_4;
+    const float density, const float u, const struct HM80_params *mat) {
 
   if (u <= 0.f) {
     return 0.f;
   }
 
-  int idx_rho, idx_u;
-  float intp_rho, intp_u;
   const float log_rho = logf(density);
   const float log_u = logf(u);
 
   // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u)
-  idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
-  idx_u = floor((log_u - mat->log_u_min) * mat->inv_log_u_step);
+  int idx_rho = floor((log_rho - mat->log_rho_min) * mat->inv_log_rho_step);
+  int idx_u = floor((log_u - mat->log_u_min) * mat->inv_log_u_step);
 
   // If outside the table then extrapolate from the edge and edge-but-one values
   if (idx_rho <= -1) {
@@ -251,26 +249,31 @@ INLINE static float HM80_pressure_from_internal_energy(
     idx_u = mat->num_u - 2;
   }
 
-  intp_rho = (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
-             mat->inv_log_rho_step;
-  intp_u =
+  const float intp_rho =
+      (log_rho - mat->log_rho_min - idx_rho * mat->log_rho_step) *
+      mat->inv_log_rho_step;
+  const float intp_u =
       (log_u - mat->log_u_min - idx_u * mat->log_u_step) * mat->inv_log_u_step;
 
   // Table values
-  log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u];
-  log_P_2 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
-  log_P_3 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u];
-  log_P_4 = mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u + 1];
-
-  log_P = (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) +
-          intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * log_P_4);
+  const float log_P_1 = mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u];
+  const float log_P_2 =
+      mat->table_log_P_rho_u[idx_rho * mat->num_u + idx_u + 1];
+  const float log_P_3 =
+      mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u];
+  const float log_P_4 =
+      mat->table_log_P_rho_u[(idx_rho + 1) * mat->num_u + idx_u + 1];
+
+  const float log_P =
+      (1.f - intp_rho) * ((1.f - intp_u) * log_P_1 + intp_u * log_P_2) +
+      intp_rho * ((1.f - intp_u) * log_P_3 + intp_u * log_P_4);
 
   return expf(log_P);
 }
 
 // gas_internal_energy_from_pressure
 INLINE static float HM80_internal_energy_from_pressure(
-    float density, float P, const struct HM80_params *mat) {
+    const float density, const float P, const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -279,9 +282,9 @@ INLINE static float HM80_internal_energy_from_pressure(
 
 // gas_soundspeed_from_internal_energy
 INLINE static float HM80_soundspeed_from_internal_energy(
-    float density, float u, const struct HM80_params *mat) {
+    const float density, const float u, const struct HM80_params *mat) {
 
-  float c, P;
+  float c;
 
   // Bulk modulus
   if (mat->bulk_mod != 0) {
@@ -289,10 +292,10 @@ INLINE static float HM80_soundspeed_from_internal_energy(
   }
   // Ideal gas
   else {
-    P = HM80_pressure_from_internal_energy(density, u, mat);
+    const float P = HM80_pressure_from_internal_energy(density, u, mat);
     c = sqrtf(hydro_gamma * P / density);
 
-    if (c <= 0) {
+    if (c <= 0.f) {
       c = sqrtf(hydro_gamma * mat->P_min_for_c_min / density);
     }
   }
@@ -302,7 +305,35 @@ INLINE static float HM80_soundspeed_from_internal_energy(
 
 // gas_soundspeed_from_pressure
 INLINE static float HM80_soundspeed_from_pressure(
-    float density, float P, const struct HM80_params *mat) {
+    const float density, const float P, const struct HM80_params *mat) {
+
+  error("This EOS function is not yet implemented!");
+
+  return 0.f;
+}
+
+// gas_entropy_from_internal_energy
+INLINE static float HM80_temperature_from_internal_energy(
+    const float density, const float u, const struct HM80_params *mat) {
+
+  error("This EOS function is not yet implemented!");
+
+  return 0.f;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float HM80_density_from_pressure_and_temperature(
+    const float P, const float T, const struct HM80_params *mat) {
+
+  error("This EOS function is not yet implemented!");
+
+  return 0.f;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float HM80_density_from_pressure_and_internal_energy(
+    const float P, const float u, const float rho_ref, const float rho_sph,
+    const struct HM80_params *mat) {
 
   error("This EOS function is not yet implemented!");
 

src/equation_of_state/planetary/ideal_gas.h

@@ -58,7 +58,7 @@ INLINE static void set_idg_def(struct idg_params *mat,
  * @param entropy The entropy \f$A\f$.
  */
 INLINE static float idg_internal_energy_from_entropy(
-    float density, float entropy, const struct idg_params *mat) {
+    const float density, const float entropy, const struct idg_params *mat) {
 
   return entropy * powf(density, mat->gamma - 1.f) *
          mat->one_over_gamma_minus_one;
@@ -72,7 +72,8 @@ INLINE static float idg_internal_energy_from_entropy(
  * @param density The density \f$\rho\f$.
  * @param entropy The entropy \f$A\f$.
  */
-INLINE static float idg_pressure_from_entropy(float density, float entropy,
+INLINE static float idg_pressure_from_entropy(const float density,
+                                              const float entropy,
                                               const struct idg_params *mat) {
 
   return entropy * powf(density, mat->gamma);
@@ -87,7 +88,8 @@ INLINE static float idg_pressure_from_entropy(float density, float entropy,
  * @param pressure The pressure \f$P\f$.
  * @return The entropy \f$A\f$.
  */
-INLINE static float idg_entropy_from_pressure(float density, float pressure,
+INLINE static float idg_entropy_from_pressure(const float density,
+                                              const float pressure,
                                               const struct idg_params *mat) {
 
   return pressure * powf(density, -mat->gamma);
@@ -101,7 +103,8 @@ INLINE static float idg_entropy_from_pressure(float density, float pressure,
  * @param density The density \f$\rho\f$.
  * @param entropy The entropy \f$A\f$.
  */
-INLINE static float idg_soundspeed_from_entropy(float density, float entropy,
+INLINE static float idg_soundspeed_from_entropy(const float density,
+                                                const float entropy,
                                                 const struct idg_params *mat) {
 
   return sqrtf(mat->gamma * powf(density, mat->gamma - 1.f) * entropy);
@@ -116,7 +119,7 @@ INLINE static float idg_soundspeed_from_entropy(float density, float entropy,
  * @param u The internal energy \f$u\f$
  */
 INLINE static float idg_entropy_from_internal_energy(
-    float density, float u, const struct idg_params *mat) {
+    const float density, const float u, const struct idg_params *mat) {
 
   return (mat->gamma - 1.f) * u * powf(density, 1.f - mat->gamma);
 }
@@ -130,7 +133,7 @@ INLINE static float idg_entropy_from_internal_energy(
  * @param u The internal energy \f$u\f$
  */
 INLINE static float idg_pressure_from_internal_energy(
-    float density, float u, const struct idg_params *mat) {
+    const float density, const float u, const struct idg_params *mat) {
 
   return (mat->gamma - 1.f) * u * density;
 }
@@ -145,7 +148,7 @@ INLINE static float idg_pressure_from_internal_energy(
  * @return The internal energy \f$u\f$.
  */
 INLINE static float idg_internal_energy_from_pressure(
-    float density, float pressure, const struct idg_params *mat) {
+    const float density, const float pressure, const struct idg_params *mat) {
 
   return mat->one_over_gamma_minus_one * pressure / density;
 }
@@ -159,7 +162,7 @@ INLINE static float idg_internal_energy_from_pressure(
  * @param u The internal energy \f$u\f$
  */
 INLINE static float idg_soundspeed_from_internal_energy(
-    float density, float u, const struct idg_params *mat) {
+    const float density, const float u, const struct idg_params *mat) {
 
   return sqrtf(u * mat->gamma * (mat->gamma - 1.f));
 }
@@ -172,10 +175,37 @@ INLINE static float idg_soundspeed_from_internal_energy(
  * @param density The density \f$\rho\f$
  * @param P The pressure \f$P\f$
  */
-INLINE static float idg_soundspeed_from_pressure(float density, float P,
+INLINE static float idg_soundspeed_from_pressure(const float density,
+                                                 const float P,
                                                  const struct idg_params *mat) {
 
   return sqrtf(mat->gamma * P / density);
 }
 
+// gas_temperature_from_internal_energy
+INLINE static float idg_temperature_from_internal_energy(
+    const float density, const float u, const struct idg_params *mat) {
+
+  error("This EOS function is not yet implemented!");
+
+  return 0.f;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float idg_density_from_pressure_and_temperature(
+    const float P, const float T, const struct idg_params *mat) {
+
+  error("This EOS function is not yet implemented!");
+
+  return 0.f;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float idg_density_from_pressure_and_internal_energy(
+    const float P, const float u, const float rho_ref, const float rho_sph,
+    const struct idg_params *mat) {
+
+  return mat->one_over_gamma_minus_one * P / u;
+}
+
 #endif /* SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H */

src/equation_of_state/planetary/sesame.h

@@ -43,6 +43,7 @@
 // SESAME parameters
 struct SESAME_params {
   float *table_log_rho;
+  float *table_log_T;
   float *table_log_u_rho_T;
   float *table_P_rho_T;
   float *table_c_rho_T;
@@ -77,6 +78,30 @@ INLINE static void set_SS08_water(struct SESAME_params *mat,
   mat->mat_id = mat_id;
   mat->version_date = 20220714;
 }
+INLINE static void set_AQUA(struct SESAME_params *mat,
+                            enum eos_planetary_material_id mat_id) {
+  // Haldemann et al. (2020)
+  mat->mat_id = mat_id;
+  mat->version_date = 20220714;
+}
+INLINE static void set_CMS19_H(struct SESAME_params *mat,
+                               enum eos_planetary_material_id mat_id) {
+  // Chabrier et al. (2019)
+  mat->mat_id = mat_id;
+  mat->version_date = 20220905;
+}
+INLINE static void set_CMS19_He(struct SESAME_params *mat,
+                                enum eos_planetary_material_id mat_id) {
+  // Chabrier et al. (2019)
+  mat->mat_id = mat_id;
+  mat->version_date = 20220905;
+}
+INLINE static void set_CD21_HHe(struct SESAME_params *mat,
+                                enum eos_planetary_material_id mat_id) {
+  // Chabrier & Debras (2021)
+  mat->mat_id = mat_id;
+  mat->version_date = 20220905;
+}
 INLINE static void set_ANEOS_forsterite(struct SESAME_params *mat,
                                         enum eos_planetary_material_id mat_id) {
   // Stewart et al. (2019)
@@ -176,6 +201,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
 
   // Allocate table memory
   mat->table_log_rho = (float *)malloc(mat->num_rho * sizeof(float));
+  mat->table_log_T = (float *)malloc(mat->num_T * sizeof(float));
   mat->table_log_u_rho_T =
       (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
   mat->table_P_rho_T =
@@ -197,10 +223,16 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
     }
   }
 
-  // Temperatures (ignored)
+  // Temperatures (not log yet)
   for (int i_T = -1; i_T < mat->num_T; i_T++) {
-    c = fscanf(f, "%f", &ignore);
-    if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+    // Ignore the first elements of rho = 0, T = 0
+    if (i_T == -1) {
+      c = fscanf(f, "%f", &ignore);
+      if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+    } else {
+      c = fscanf(f, "%f", &mat->table_log_T[i_T]);
+      if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
+    }
   }
 
   // Sp. int. energies (not log yet), pressures, sound speeds, and sp.
@@ -232,6 +264,10 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
   for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
     mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]);
   }
+  // Convert temperatures to log(temperature)
+  for (int i_T = 0; i_T < mat->num_T; i_T++) {
+    mat->table_log_T[i_T] = logf(mat->table_log_T[i_T]);
+  }
 
   // Initialise tiny values
   mat->u_tiny = FLT_MAX;
@@ -297,6 +333,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
 
       mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] =
           logf(mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
+
+      // Ensure P > 0
+      if (mat->table_P_rho_T[i_rho * mat->num_T + i_T] <= 0) {
+        mat->table_P_rho_T[i_rho * mat->num_T + i_T] = mat->P_tiny;
+      }
     }
   }
 }
@@ -316,6 +357,13 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
              units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
   }
 
+  // Temperatures (log)
+  for (int i_T = 0; i_T < mat->num_T; i_T++) {
+    mat->table_log_T[i_T] +=
+        logf(units_cgs_conversion_factor(&si, UNIT_CONV_TEMPERATURE) /
+             units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE));
+  }
+
   // Sp. int. energies (log), pressures, sound speeds, and sp. entropies
   for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
     for (int i_T = 0; i_T < mat->num_T; i_T++) {
@@ -352,28 +400,26 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
 
 // gas_internal_energy_from_entropy
 INLINE static float SESAME_internal_energy_from_entropy(
-    float density, float entropy, const struct SESAME_params *mat) {
-
-  float u, log_u_1, log_u_2, log_u_3, log_u_4;
+    const float density, const float entropy, const struct SESAME_params *mat) {
 
+  // Return zero if entropy is zero
   if (entropy <= 0.f) {
     return 0.f;
   }
 
-  int idx_rho, idx_s_1, idx_s_2;
-  float intp_rho, intp_s_1, intp_s_2;
   const float log_rho = logf(density);
   const float log_s = logf(entropy);
 
-  // 2D interpolation (bilinear with log(rho), log(s)) to find u(rho, s))
+  // 2D interpolation (bilinear with log(rho), log(s) to find u(rho, s))
+
   // Density index
-  idx_rho =
+  int idx_rho =
       find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
 
   // Sp. entropy at this and the next density (in relevant slice of s array)
-  idx_s_1 = find_value_in_monot_incr_array(
+  int idx_s_1 = find_value_in_monot_incr_array(
       log_s, mat->table_log_s_rho_T + idx_rho * mat->num_T, mat->num_T);
-  idx_s_2 = find_value_in_monot_incr_array(
+  int idx_s_2 = find_value_in_monot_incr_array(
       log_s, mat->table_log_s_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
 
   // If outside the table then extrapolate from the edge and edge-but-one values
@@ -394,6 +440,7 @@ INLINE static float SESAME_internal_energy_from_entropy(
   }
 
   // Check for duplicates in SESAME tables before interpolation
+  float intp_rho, intp_s_1, intp_s_2;
   if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
     intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
                (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
@@ -420,10 +467,13 @@ INLINE static float SESAME_internal_energy_from_entropy(
   }
 
   // Table values
-  log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
-  log_u_2 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
-  log_u_3 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
-  log_u_4 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
+  const float log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
+  const float log_u_2 =
+      mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
+  const float log_u_3 =
+      mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
+  const float log_u_4 =
+      mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
 
   // If below the minimum s at this rho then just use the lowest table values
   if ((idx_rho > 0.f) && ((intp_s_1 < 0.f) || (intp_s_2 < 0.f) ||
@@ -433,18 +483,17 @@ INLINE static float SESAME_internal_energy_from_entropy(
   }
 
   // Interpolate with the log values
-  u = (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
+  const float u =
+      (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
       intp_rho * ((1.f - intp_s_2) * log_u_3 + intp_s_2 * log_u_4);
 
   // Convert back from log
-  u = expf(u);
-
-  return u;
+  return expf(u);
 }
 
 // gas_pressure_from_entropy
 INLINE static float SESAME_pressure_from_entropy(
-    float density, float entropy, const struct SESAME_params *mat) {
+    const float density, const float entropy, const struct SESAME_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -453,7 +502,8 @@ INLINE static float SESAME_pressure_from_entropy(
 
 // gas_entropy_from_pressure
 INLINE static float SESAME_entropy_from_pressure(
-    float density, float pressure, const struct SESAME_params *mat) {
+    const float density, const float pressure,
+    const struct SESAME_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -462,7 +512,7 @@ INLINE static float SESAME_entropy_from_pressure(
 
 // gas_soundspeed_from_entropy
 INLINE static float SESAME_soundspeed_from_entropy(
-    float density, float entropy, const struct SESAME_params *mat) {
+    const float density, const float entropy, const struct SESAME_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -471,35 +521,32 @@ INLINE static float SESAME_soundspeed_from_entropy(
 
 // gas_entropy_from_internal_energy
 INLINE static float SESAME_entropy_from_internal_energy(
-    float density, float u, const struct SESAME_params *mat) {
+    const float density, const float u, const struct SESAME_params *mat) {
 
   return 0.f;
 }
 
 // gas_pressure_from_internal_energy
 INLINE static float SESAME_pressure_from_internal_energy(
-    float density, float u, const struct SESAME_params *mat) {
-
-  float P, P_1, P_2, P_3, P_4;
+    const float density, const float u, const struct SESAME_params *mat) {
 
   if (u <= 0.f) {
     return 0.f;
   }
 
-  int idx_rho, idx_u_1, idx_u_2;
-  float intp_rho, intp_u_1, intp_u_2;
   const float log_rho = logf(density);
   const float log_u = logf(u);
 
-  // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u))
+  // 2D interpolation (bilinear with log(rho), log(u) to find P(rho, u))
+
   // Density index
-  idx_rho =
+  int idx_rho =
       find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
 
   // Sp. int. energy at this and the next density (in relevant slice of u array)
-  idx_u_1 = find_value_in_monot_incr_array(
+  int idx_u_1 = find_value_in_monot_incr_array(
       log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
-  idx_u_2 = find_value_in_monot_incr_array(
+  int idx_u_2 = find_value_in_monot_incr_array(
       log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
 
   // If outside the table then extrapolate from the edge and edge-but-one values
@@ -520,18 +567,19 @@ INLINE static float SESAME_pressure_from_internal_energy(
   }
 
   // Check for duplicates in SESAME tables before interpolation
+  float intp_rho, intp_u_1, intp_u_2;
+  const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
   if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
     intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
                (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
   } else {
     intp_rho = 1.f;
   }
-  if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
-      mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
-    intp_u_1 =
-        (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
-        (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
-         mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+  if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+      mat->table_log_u_rho_T[idx_u_rho_T]) {
+    intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+               (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+                mat->table_log_u_rho_T[idx_u_rho_T]);
   } else {
     intp_u_1 = 1.f;
   }
@@ -546,10 +594,10 @@ INLINE static float SESAME_pressure_from_internal_energy(
   }
 
   // Table values
-  P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1];
-  P_2 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
-  P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
-  P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+  float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+  float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+  float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+  float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
 
   // If below the minimum u at this rho then just use the lowest table values
   if ((idx_rho > 0.f) &&
@@ -582,19 +630,18 @@ INLINE static float SESAME_pressure_from_internal_energy(
   P_2 = logf(P_2);
   P_3 = logf(P_3);
   P_4 = logf(P_4);
+  const float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+  const float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
 
-  P = (1.f - intp_rho) * ((1.f - intp_u_1) * P_1 + intp_u_1 * P_2) +
-      intp_rho * ((1.f - intp_u_2) * P_3 + intp_u_2 * P_4);
+  const float P = (1.f - intp_rho) * P_1_2 + intp_rho * P_3_4;
 
   // Convert back from log
-  P = expf(P);
-
-  return P;
+  return expf(P);
 }
 
 // gas_internal_energy_from_pressure
 INLINE static float SESAME_internal_energy_from_pressure(
-    float density, float P, const struct SESAME_params *mat) {
+    const float density, const float P, const struct SESAME_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -603,28 +650,26 @@ INLINE static float SESAME_internal_energy_from_pressure(
 
 // gas_soundspeed_from_internal_energy
 INLINE static float SESAME_soundspeed_from_internal_energy(
-    float density, float u, const struct SESAME_params *mat) {
-
-  float c, c_1, c_2, c_3, c_4;
+    const float density, const float u, const struct SESAME_params *mat) {
 
+  // Return zero if internal energy is non-positive
   if (u <= 0.f) {
     return 0.f;
   }
 
-  int idx_rho, idx_u_1, idx_u_2;
-  float intp_rho, intp_u_1, intp_u_2;
   const float log_rho = logf(density);
   const float log_u = logf(u);
 
-  // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u))
+  // 2D interpolation (bilinear with log(rho), log(u) to find c(rho, u))
+
   // Density index
-  idx_rho =
+  int idx_rho =
       find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
 
   // Sp. int. energy at this and the next density (in relevant slice of u array)
-  idx_u_1 = find_value_in_monot_incr_array(
+  int idx_u_1 = find_value_in_monot_incr_array(
       log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
-  idx_u_2 = find_value_in_monot_incr_array(
+  int idx_u_2 = find_value_in_monot_incr_array(
       log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
 
   // If outside the table then extrapolate from the edge and edge-but-one values
@@ -645,18 +690,19 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
   }
 
   // Check for duplicates in SESAME tables before interpolation
+  float intp_rho, intp_u_1, intp_u_2;
   if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
     intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
                (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
   } else {
     intp_rho = 1.f;
   }
-  if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
-      mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
-    intp_u_1 =
-        (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
-        (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
-         mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
+  const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+  if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+      mat->table_log_u_rho_T[idx_u_rho_T]) {
+    intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+               (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+                mat->table_log_u_rho_T[idx_u_rho_T]);
   } else {
     intp_u_1 = 1.f;
   }
@@ -671,10 +717,10 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
   }
 
   // Table values
-  c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];
-  c_2 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1 + 1];
-  c_3 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
-  c_4 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+  float c_1 = mat->table_c_rho_T[idx_u_rho_T];
+  float c_2 = mat->table_c_rho_T[idx_u_rho_T + 1];
+  float c_3 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+  float c_4 = mat->table_c_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
 
   // If more than two table values are non-positive then return zero
   int num_non_pos = 0;
@@ -710,22 +756,536 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
     intp_u_2 = 0;
   }
 
-  c = (1.f - intp_rho) * ((1.f - intp_u_1) * c_1 + intp_u_1 * c_2) +
-      intp_rho * ((1.f - intp_u_2) * c_3 + intp_u_2 * c_4);
+  const float c = (1.f - intp_rho) * ((1.f - intp_u_1) * c_1 + intp_u_1 * c_2) +
+                  intp_rho * ((1.f - intp_u_2) * c_3 + intp_u_2 * c_4);
 
   // Convert back from log
-  c = expf(c);
-
-  return c;
+  return expf(c);
 }
 
 // gas_soundspeed_from_pressure
 INLINE static float SESAME_soundspeed_from_pressure(
-    float density, float P, const struct SESAME_params *mat) {
+    const float density, const float P, const struct SESAME_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
   return 0.f;
 }
 
+// gas_temperature_from_internal_energy
+INLINE static float SESAME_temperature_from_internal_energy(
+    const float density, const float u, const struct SESAME_params *mat) {
+
+  // Return zero if zero internal energy
+  if (u <= 0.f) {
+    return 0.f;
+  }
+
+  const float log_rho = logf(density);
+  const float log_u = logf(u);
+
+  // 2D interpolation (bilinear with log(rho), log(u) to find T(rho, u)
+
+  // Density index
+  int idx_rho =
+      find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
+
+  // Sp. int. energy at this and the next density (in relevant slice of u array)
+  int idx_u_1 = find_value_in_monot_incr_array(
+      log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+  int idx_u_2 = find_value_in_monot_incr_array(
+      log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+  // If outside the table then extrapolate from the edge and edge-but-one values
+  if (idx_rho <= -1) {
+    idx_rho = 0;
+  } else if (idx_rho >= mat->num_rho) {
+    idx_rho = mat->num_rho - 2;
+  }
+  if (idx_u_1 <= -1) {
+    idx_u_1 = 0;
+  } else if (idx_u_1 >= mat->num_T) {
+    idx_u_1 = mat->num_T - 2;
+  }
+  if (idx_u_2 <= -1) {
+    idx_u_2 = 0;
+  } else if (idx_u_2 >= mat->num_T) {
+    idx_u_2 = mat->num_T - 2;
+  }
+
+  // Check for duplicates in SESAME tables before interpolation
+  float intp_u_1, intp_u_2;
+  const int idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+  if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+      mat->table_log_u_rho_T[idx_u_rho_T]) {
+    intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+               (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+                mat->table_log_u_rho_T[idx_u_rho_T]);
+  } else {
+    intp_u_1 = 1.f;
+  }
+  if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+      mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+    intp_u_2 =
+        (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+        (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+         mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+  } else {
+    intp_u_2 = 1.f;
+  }
+
+  // Compute line points
+  float log_rho_1 = mat->table_log_rho[idx_rho];
+  float log_rho_2 = mat->table_log_rho[idx_rho + 1];
+  float log_T_1 = mat->table_log_T[idx_u_1];
+  log_T_1 +=
+      intp_u_1 * (mat->table_log_T[idx_u_1 + 1] - mat->table_log_T[idx_u_1]);
+  float log_T_2 = mat->table_log_T[idx_u_2];
+  log_T_2 +=
+      intp_u_2 * (mat->table_log_T[idx_u_2 + 1] - mat->table_log_T[idx_u_2]);
+
+  // Intersect line passing through (log_rho_1, log_T_1), (log_rho_2, log_T_2)
+  // with line density = log_rho
+
+  // Check for log_T_1 == log_T_2
+  float log_T;
+  if (log_T_1 == log_T_2) {
+    log_T = log_T_1;
+  } else {
+    // log_rho = slope*log_T + intercept
+    const float slope = (log_rho_1 - log_rho_2) / (log_T_1 - log_T_2);
+    const float intercept = log_rho_1 - slope * log_T_1;
+    log_T = (log_rho - intercept) / slope;
+  }
+
+  // Convert back from log
+  return expf(log_T);
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float SESAME_density_from_pressure_and_temperature(
+    float P, float T, const struct SESAME_params *mat) {
+
+  // Return zero if pressure is non-positive
+  if (P <= 0.f) {
+    return 0.f;
+  }
+
+  const float log_T = logf(T);
+
+  // 2D interpolation (bilinear with log(T), P to find rho(T, P))
+
+  // Temperature index
+  int idx_T =
+      find_value_in_monot_incr_array(log_T, mat->table_log_T, mat->num_T);
+
+  // If outside the table then extrapolate from the edge and edge-but-one values
+  if (idx_T <= -1) {
+    idx_T = 0;
+  } else if (idx_T >= mat->num_T) {
+    idx_T = mat->num_T - 2;
+  }
+
+  // Pressure at this and the next temperature (in relevant vertical slice of P
+  // array)
+  int idx_P_1 = vertical_find_value_in_monot_incr_array(
+      P, mat->table_P_rho_T, mat->num_rho, mat->num_T, idx_T);
+  int idx_P_2 = vertical_find_value_in_monot_incr_array(
+      P, mat->table_P_rho_T, mat->num_rho, mat->num_T, idx_T + 1);
+
+  // If outside the table then extrapolate from the edge and edge-but-one values
+  if (idx_P_1 <= -1) {
+    idx_P_1 = 0;
+  } else if (idx_P_1 >= mat->num_rho) {
+    idx_P_1 = mat->num_rho - 2;
+  }
+  if (idx_P_2 <= -1) {
+    idx_P_2 = 0;
+  } else if (idx_P_2 >= mat->num_rho) {
+    idx_P_2 = mat->num_rho - 2;
+  }
+
+  // Check for duplicates in SESAME tables before interpolation
+  float intp_P_1, intp_P_2;
+  if (mat->table_P_rho_T[(idx_P_1 + 1) * mat->num_T + idx_T] !=
+      mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]) {
+    intp_P_1 = (P - mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]) /
+               (mat->table_P_rho_T[(idx_P_1 + 1) * mat->num_T + idx_T] -
+                mat->table_P_rho_T[idx_P_1 * mat->num_T + idx_T]);
+  } else {
+    intp_P_1 = 1.f;
+  }
+  if (mat->table_P_rho_T[(idx_P_2 + 1) * mat->num_T + (idx_T + 1)] !=
+      mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]) {
+    intp_P_2 = (P - mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]) /
+               (mat->table_P_rho_T[(idx_P_2 + 1) * mat->num_T + (idx_T + 1)] -
+                mat->table_P_rho_T[idx_P_2 * mat->num_T + (idx_T + 1)]);
+  } else {
+    intp_P_2 = 1.f;
+  }
+
+  // Compute line points
+  const float log_T_1 = mat->table_log_T[idx_T];
+  const float log_T_2 = mat->table_log_T[idx_T + 1];
+  const float log_rho_1 = (1.f - intp_P_1) * mat->table_log_rho[idx_P_1] +
+                          intp_P_1 * mat->table_log_rho[idx_P_1 + 1];
+  const float log_rho_2 = (1.f - intp_P_2) * mat->table_log_rho[idx_P_2] +
+                          intp_P_2 * mat->table_log_rho[idx_P_2 + 1];
+
+  // Intersect line passing through (log_rho_1, log_T_1), (log_rho_2, log_T_2)
+  // with line temperature = log_T
+
+  // Check for log_rho_1 == log_rho_2
+  float log_rho;
+  if (log_rho_1 == log_rho_2) {
+    log_rho = log_rho_1;
+  } else {
+    // log_T = slope*log_rho + intercept
+    const float slope = (log_T_1 - log_T_2) / (log_rho_1 - log_rho_2);
+    const float intercept = log_T_1 - slope * log_rho_1;
+    log_rho = (log_T - intercept) / slope;
+  }
+
+  // Convert back from log
+  return expf(log_rho);
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float SESAME_density_from_pressure_and_internal_energy(
+    float P, float u, float rho_ref, float rho_sph,
+    const struct SESAME_params *mat) {
+
+  // Return the unchanged density if u or P is non-positive
+  if (u <= 0.f || P <= 0.f) {
+    return rho_sph;
+  }
+
+  // Convert inputs to log
+  const float log_u = logf(u);
+  const float log_P = logf(P);
+  const float log_rho_ref = logf(rho_ref);
+
+  // Find rounded down index of reference density. This is where we start our
+  // search
+  int idx_rho_ref = find_value_in_monot_incr_array(
+      log_rho_ref, mat->table_log_rho, mat->num_rho);
+
+  // If no roots are found in the current search range, we increase search range
+  // by search_factor_log_rho above and below the reference density each
+  // iteration.
+  const float search_factor_log_rho = logf(10.f);
+
+  // Initialise the minimum and maximum densities we're searching to at the
+  // reference density. These will change before the first iteration.
+  float log_rho_min = log_rho_ref;
+  float log_rho_max = log_rho_ref;
+
+  // Initialise search indices around rho_ref
+  int idx_rho_below_min, idx_rho_above_max;
+
+  // If we find a root, it will get stored as closest_root
+  float closest_root = -1.f;
+  float root_below;
+
+  // Initialise pressures
+  float P_above_lower, P_above_upper;
+  float P_below_lower, P_below_upper;
+  P_above_upper = 0.f;
+  P_below_lower = 0.f;
+
+  // Increase search range by search_factor_log_rho
+  log_rho_max += search_factor_log_rho;
+  idx_rho_above_max = find_value_in_monot_incr_array(
+      log_rho_max, mat->table_log_rho, mat->num_rho);
+  log_rho_min -= search_factor_log_rho;
+  idx_rho_below_min = find_value_in_monot_incr_array(
+      log_rho_min, mat->table_log_rho, mat->num_rho);
+
+  // When searching above/below, we are looking for where the pressure P(rho, u)
+  // of the table densities changes from being less than to more than, or vice
+  // versa, the desired pressure. If this is the case, there is a root between
+  // these table values of rho.
+
+  // First look for roots above rho_ref
+  int idx_u_rho_T;
+  float intp_rho;
+  for (int idx_rho = idx_rho_ref; idx_rho <= idx_rho_above_max; idx_rho++) {
+
+    // This is similar to P_u_rho, but we're not interested in intp_rho,
+    // and instead calculate the pressure for both intp_rho=0 and intp_rho=1
+
+    // Sp. int. energy at this and the next density (in relevant slice of u
+    // array)
+    int idx_u_1 = find_value_in_monot_incr_array(
+        log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+    int idx_u_2 = find_value_in_monot_incr_array(
+        log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+    // If outside the table then extrapolate from the edge and edge-but-one
+    // values
+    if (idx_rho <= -1) {
+      idx_rho = 0;
+    } else if (idx_rho >= mat->num_rho) {
+      idx_rho = mat->num_rho - 2;
+    }
+    if (idx_u_1 <= -1) {
+      idx_u_1 = 0;
+    } else if (idx_u_1 >= mat->num_T) {
+      idx_u_1 = mat->num_T - 2;
+    }
+    if (idx_u_2 <= -1) {
+      idx_u_2 = 0;
+    } else if (idx_u_2 >= mat->num_T) {
+      idx_u_2 = mat->num_T - 2;
+    }
+
+    float intp_u_1, intp_u_2;
+    idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+    if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+        mat->table_log_u_rho_T[idx_u_rho_T]) {
+      intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+                 (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+                  mat->table_log_u_rho_T[idx_u_rho_T]);
+    } else {
+      intp_u_1 = 1.f;
+    }
+    if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+        mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+      intp_u_2 =
+          (log_u -
+           mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+          (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+           mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+    } else {
+      intp_u_2 = 1.f;
+    }
+
+    // Table values
+    float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+    float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+    float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+    float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+    // If below the minimum u at this rho then just use the lowest table values
+    if ((idx_rho > 0.f) &&
+        ((intp_u_1 < 0.f) || (intp_u_2 < 0.f) || (P_1 > P_2) || (P_3 > P_4))) {
+      intp_u_1 = 0;
+      intp_u_2 = 0;
+    }
+
+    // If more than two table values are non-positive then return zero
+    int num_non_pos = 0;
+    if (P_1 <= 0.f) num_non_pos++;
+    if (P_2 <= 0.f) num_non_pos++;
+    if (P_3 <= 0.f) num_non_pos++;
+    if (P_4 <= 0.f) num_non_pos++;
+    if (num_non_pos > 0) {
+      // If just one or two are non-positive then replace them with a tiny value
+      // Unless already trying to extrapolate in which case return zero
+      if ((num_non_pos > 2) || (mat->P_tiny == 0.f) || (intp_u_1 < 0.f) ||
+          (intp_u_2 < 0.f)) {
+        break;  // return rho_sph;
+      }
+      if (P_1 <= 0.f) P_1 = mat->P_tiny;
+      if (P_2 <= 0.f) P_2 = mat->P_tiny;
+      if (P_3 <= 0.f) P_3 = mat->P_tiny;
+      if (P_4 <= 0.f) P_4 = mat->P_tiny;
+    }
+
+    // Interpolate with the log values
+    P_1 = logf(P_1);
+    P_2 = logf(P_2);
+    P_3 = logf(P_3);
+    P_4 = logf(P_4);
+    float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+    float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
+
+    // Pressure for intp_rho = 0
+    P_above_lower = expf(P_1_2);
+
+    // Because of linear interpolation, pressures are not exactly continuous
+    // as we go from one side of a grid point to another. See if there is
+    // a root between the last P_above_upper and the new P_above_lower,
+    // which are approx the same.
+    if (idx_rho != idx_rho_ref) {
+      if ((P_above_lower - P) * (P_above_upper - P) <= 0) {
+        closest_root = expf(mat->table_log_rho[idx_rho]);
+        break;
+      }
+    }
+
+    // Pressure for intp_rho = 1
+    P_above_upper = expf(P_3_4);
+
+    // Does the pressure of the adjacent table densities switch from being
+    // above to below the desired pressure, or vice versa? If so, there is a
+    // root.
+    if ((P_above_lower - P) * (P_above_upper - P) <= 0.f) {
+
+      // If there is a root, interpolate between the table values:
+      intp_rho = (log_P - P_1_2) / (P_3_4 - P_1_2);
+
+      closest_root = expf(mat->table_log_rho[idx_rho] +
+                          intp_rho * (mat->table_log_rho[idx_rho + 1] -
+                                      mat->table_log_rho[idx_rho]));
+
+      // If the root is between the same table values as the reference value,
+      // then this is the closest root, so we can return it without further
+      // searching
+      if (idx_rho == idx_rho_ref) {
+        return closest_root;
+      }
+
+      // Found a root, so no need to search higher densities
+      break;
+    }
+  }
+
+  // If we found a root above, change search range below so that we're only
+  // looking for closer (in log) roots than the one we found
+  if (closest_root > 0.f) {
+    log_rho_min = log_rho_ref - (logf(closest_root) - log_rho_ref);
+    idx_rho_below_min = find_value_in_monot_incr_array(
+        log_rho_min, mat->table_log_rho, mat->num_rho);
+  }
+
+  // Now look for roots below rho_ref
+  for (int idx_rho = idx_rho_ref; idx_rho >= idx_rho_below_min; idx_rho--) {
+
+    // Sp. int. energy at this and the next density (in relevant slice of u
+    // array)
+    int idx_u_1 = find_value_in_monot_incr_array(
+        log_u, mat->table_log_u_rho_T + idx_rho * mat->num_T, mat->num_T);
+    int idx_u_2 = find_value_in_monot_incr_array(
+        log_u, mat->table_log_u_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
+
+    // If outside the table then extrapolate from the edge and edge-but-one
+    // values
+    if (idx_rho <= -1) {
+      idx_rho = 0;
+    } else if (idx_rho >= mat->num_rho) {
+      idx_rho = mat->num_rho - 2;
+    }
+    if (idx_u_1 <= -1) {
+      idx_u_1 = 0;
+    } else if (idx_u_1 >= mat->num_T) {
+      idx_u_1 = mat->num_T - 2;
+    }
+    if (idx_u_2 <= -1) {
+      idx_u_2 = 0;
+    } else if (idx_u_2 >= mat->num_T) {
+      idx_u_2 = mat->num_T - 2;
+    }
+
+    idx_u_rho_T = idx_rho * mat->num_T + idx_u_1;
+    float intp_u_1, intp_u_2;
+    if (mat->table_log_u_rho_T[idx_u_rho_T + 1] !=
+        mat->table_log_u_rho_T[idx_u_rho_T]) {
+      intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_u_rho_T]) /
+                 (mat->table_log_u_rho_T[idx_u_rho_T + 1] -
+                  mat->table_log_u_rho_T[idx_u_rho_T]);
+    } else {
+      intp_u_1 = 1.f;
+    }
+    if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
+        mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
+      intp_u_2 =
+          (log_u -
+           mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
+          (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
+           mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
+    } else {
+      intp_u_2 = 1.f;
+    }
+
+    // Table values
+    float P_1 = mat->table_P_rho_T[idx_u_rho_T];
+    float P_2 = mat->table_P_rho_T[idx_u_rho_T + 1];
+    float P_3 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2];
+    float P_4 = mat->table_P_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2 + 1];
+
+    // If below the minimum u at this rho then just use the lowest table values
+    if ((idx_rho > 0.f) &&
+        ((intp_u_1 < 0.f) || (intp_u_2 < 0.f) || (P_1 > P_2) || (P_3 > P_4))) {
+      intp_u_1 = 0;
+      intp_u_2 = 0;
+    }
+
+    // If more than two table values are non-positive then return zero
+    int num_non_pos = 0;
+    if (P_1 <= 0.f) num_non_pos++;
+    if (P_2 <= 0.f) num_non_pos++;
+    if (P_3 <= 0.f) num_non_pos++;
+    if (P_4 <= 0.f) num_non_pos++;
+    if (num_non_pos > 0) {
+      // If just one or two are non-positive then replace them with a tiny value
+      // Unless already trying to extrapolate in which case return zero
+      if ((num_non_pos > 2) || (mat->P_tiny == 0.f) || (intp_u_1 < 0.f) ||
+          (intp_u_2 < 0.f)) {
+        break;
+      }
+      if (P_1 <= 0.f) P_1 = mat->P_tiny;
+      if (P_2 <= 0.f) P_2 = mat->P_tiny;
+      if (P_3 <= 0.f) P_3 = mat->P_tiny;
+      if (P_4 <= 0.f) P_4 = mat->P_tiny;
+    }
+
+    // Interpolate with the log values
+    P_1 = logf(P_1);
+    P_2 = logf(P_2);
+    P_3 = logf(P_3);
+    P_4 = logf(P_4);
+    const float P_1_2 = (1.f - intp_u_1) * P_1 + intp_u_1 * P_2;
+    const float P_3_4 = (1.f - intp_u_2) * P_3 + intp_u_2 * P_4;
+
+    // Pressure for intp_rho = 1
+    P_below_upper = expf(P_3_4);
+    // Because of linear interpolation, pressures are not exactly continuous
+    // as we go from one side of a grid point to another. See if there is
+    // a root between the last P_below_lower and the new P_below_upper,
+    // which are approx the same.
+    if (idx_rho != idx_rho_ref) {
+      if ((P_below_lower - P) * (P_below_upper - P) <= 0) {
+        closest_root = expf(mat->table_log_rho[idx_rho + 1]);
+        break;
+      }
+    }
+    // Pressure for intp_rho = 0
+    P_below_lower = expf(P_1_2);
+
+    // Does the pressure of the adjacent table densities switch from being
+    // above to below the desired pressure, or vice versa? If so, there is a
+    // root.
+    if ((P_below_lower - P) * (P_below_upper - P) <= 0.f) {
+
+      // If there is a root, interpolate between the table values:
+      intp_rho = (log_P - P_1_2) / (P_3_4 - P_1_2);
+
+      root_below = expf(mat->table_log_rho[idx_rho] +
+                        intp_rho * (mat->table_log_rho[idx_rho + 1] -
+                                    mat->table_log_rho[idx_rho]));
+
+      // If we found a root above, which one is closer to the reference rho?
+      if (closest_root > 0.f) {
+        if (fabsf(logf(root_below) - logf(rho_ref)) <
+            fabsf(logf(closest_root) - logf(rho_ref))) {
+          closest_root = root_below;
+        }
+      } else {
+        closest_root = root_below;
+      }
+      // Found a root, so no need to search higher densities
+      break;
+    }
+  }
+
+  // Return the root if we found one
+  if (closest_root > 0.f) {
+    return closest_root;
+  }
+
+  // If we don't find a root before we reach max_counter, return rho_ref. Maybe
+  // we should give an error here?
+  return rho_sph;
+}
 #endif /* SWIFT_SESAME_EQUATION_OF_STATE_H */

src/equation_of_state/planetary/tillotson.h

@@ -37,12 +37,16 @@
 #include "equation_of_state.h"
 #include "inline.h"
 #include "physical_constants.h"
+#include "sesame.h"
 #include "units.h"
 
 // Tillotson parameters
 struct Til_params {
-  float rho_0, a, b, A, B, u_0, u_iv, u_cv, alpha, beta, eta_min, eta_zero,
-      P_min;
+  float rho_0, a, b, A, B, u_0, u_iv, u_cv, alpha, beta;
+  float eta_min, eta_zero, P_min, C_V;
+  float *A1_u_cold;
+  float rho_cold_min, rho_cold_max;
+  float rho_min, rho_max;
   enum eos_planetary_material_id mat_id;
 };
 
@@ -62,7 +66,12 @@ INLINE static void set_Til_iron(struct Til_params *mat,
   mat->beta = 5.0f;
   mat->eta_min = 0.0f;
   mat->eta_zero = 0.0f;
-  mat->P_min = 0.0f;
+  mat->P_min = 0.01f;
+  mat->C_V = 449.0f;
+  mat->rho_cold_min = 100.0f;
+  mat->rho_cold_max = 1.0e5f;
+  mat->rho_min = 1.0f;
+  mat->rho_max = 1.0e5f;
 }
 INLINE static void set_Til_granite(struct Til_params *mat,
                                    enum eos_planetary_material_id mat_id) {
@@ -79,7 +88,12 @@ INLINE static void set_Til_granite(struct Til_params *mat,
   mat->beta = 5.0f;
   mat->eta_min = 0.0f;
   mat->eta_zero = 0.0f;
-  mat->P_min = 0.0f;
+  mat->P_min = 0.01f;
+  mat->C_V = 790.0f;
+  mat->rho_cold_min = 100.0f;
+  mat->rho_cold_max = 1.0e5f;
+  mat->rho_min = 1.0f;
+  mat->rho_max = 1.0e5f;
 }
 INLINE static void set_Til_basalt(struct Til_params *mat,
                                   enum eos_planetary_material_id mat_id) {
@@ -96,7 +110,12 @@ INLINE static void set_Til_basalt(struct Til_params *mat,
   mat->beta = 5.0f;
   mat->eta_min = 0.0f;
   mat->eta_zero = 0.0f;
-  mat->P_min = 0.0f;
+  mat->P_min = 0.01f;
+  mat->C_V = 790.0f;
+  mat->rho_cold_min = 100.0f;
+  mat->rho_cold_max = 1.0e5f;
+  mat->rho_min = 1.0f;
+  mat->rho_max = 1.0e5f;
 }
 INLINE static void set_Til_water(struct Til_params *mat,
                                  enum eos_planetary_material_id mat_id) {
@@ -113,7 +132,77 @@ INLINE static void set_Til_water(struct Til_params *mat,
   mat->beta = 5.0f;
   mat->eta_min = 0.925f;
   mat->eta_zero = 0.875f;
+  mat->P_min = 0.01f;
+  mat->C_V = 4186.0f;
+  mat->rho_cold_min = 100.0f;
+  mat->rho_cold_max = 1.0e5f;
+  mat->rho_min = 1.0f;
+  mat->rho_max = 1.0e5f;
+}
+INLINE static void set_Til_ice(struct Til_params *mat,
+                               enum eos_planetary_material_id mat_id) {
+  mat->mat_id = mat_id;
+  mat->rho_0 = 1293.0f;
+  mat->a = 0.3f;
+  mat->b = 0.1f;
+  mat->A = 1.07e10f;
+  mat->B = 6.5e10f;
+  mat->u_0 = 1.0e7f;
+  mat->u_iv = 7.73e5f;
+  mat->u_cv = 3.04e6f;
+  mat->alpha = 10.0f;
+  mat->beta = 5.0f;
+  mat->eta_min = 0.925f;
+  mat->eta_zero = 0.875f;
   mat->P_min = 0.0f;
+  mat->C_V = 2093.0f;
+  mat->rho_cold_min = 100.0f;
+  mat->rho_cold_max = 1.0e5f;
+  mat->rho_min = 1.0f;
+  mat->rho_max = 1.0e5f;
+}
+
+/*
+    Read the parameters from a file.
+
+    File contents
+    -------------
+    # header (5 lines)
+    rho_0 (kg/m3)  a (-)  b (-)  A (Pa)  B (Pa)
+    u_0 (J)  u_iv (J)  u_cv (J)  alpha (-)  beta (-)
+    eta_min (-)  eta_zero (-)  P_min (Pa)  C_V (J kg^-1 K^-1)
+    rho_cold_min (kg/m3)  rho_cold_max (kg/m3)  rho_min (kg/m3)  rho_max (kg/m3)
+*/
+INLINE static void set_Til_custom(struct Til_params *mat,
+                                  enum eos_planetary_material_id mat_id,
+                                  char *param_file) {
+  mat->mat_id = mat_id;
+
+  // Load table contents from file
+  FILE *f = fopen(param_file, "r");
+  if (f == NULL)
+    error("Failed to open the Tillotson EoS file '%s'", param_file);
+
+  // Skip header lines
+  skip_lines(f, 5);
+
+  // Read parameters (SI)
+  int c;
+  c = fscanf(f, "%f %f %f %f %f", &mat->rho_0, &mat->a, &mat->b, &mat->A,
+             &mat->B);
+  if (c != 5) error("Failed to read the Tillotson EoS file %s", param_file);
+
+  c = fscanf(f, "%f %f %f %f %f", &mat->u_0, &mat->u_iv, &mat->u_cv,
+             &mat->alpha, &mat->beta);
+  if (c != 5) error("Failed to read the Tillotson EoS file %s", param_file);
+
+  c = fscanf(f, "%f %f %f %f", &mat->eta_min, &mat->eta_zero, &mat->P_min,
+             &mat->C_V);
+  if (c != 4) error("Failed to read the Tillotson EoS file %s", param_file);
+
+  c = fscanf(f, "%f %f %f %f", &mat->rho_cold_min, &mat->rho_cold_max,
+             &mat->rho_min, &mat->rho_max);
+  if (c != 4) error("Failed to read the Tillotson EoS file %s", param_file);
 }
 
 // Convert to internal units
@@ -123,6 +212,8 @@ INLINE static void convert_units_Til(struct Til_params *mat,
   struct unit_system si;
   units_init_si(&si);
 
+  const int N = 10000;
+
   // SI to cgs
   mat->rho_0 *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
   mat->A *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
@@ -132,6 +223,19 @@ INLINE static void convert_units_Til(struct Til_params *mat,
   mat->u_cv *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   mat->P_min *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE);
 
+  for (int i = 0; i < N; i++) {
+    mat->A1_u_cold[i] *=
+        units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  }
+
+  mat->C_V *= units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  // Entropy units don't work? using internal kelvin
+
+  mat->rho_cold_min *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+  mat->rho_cold_max *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+  mat->rho_min *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+  mat->rho_max *= units_cgs_conversion_factor(&si, UNIT_CONV_DENSITY);
+
   // cgs to internal
   mat->rho_0 /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
   mat->A /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
@@ -140,11 +244,24 @@ INLINE static void convert_units_Til(struct Til_params *mat,
   mat->u_iv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   mat->u_cv /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   mat->P_min /= units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
+
+  for (int i = 0; i < N; i++) {
+    mat->A1_u_cold[i] /=
+        units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  }
+
+  mat->C_V /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+  // Entropy units don't work? using internal kelvin
+
+  mat->rho_cold_min /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  mat->rho_cold_max /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  mat->rho_min /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+  mat->rho_max /= units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
 }
 
 // gas_internal_energy_from_entropy
 INLINE static float Til_internal_energy_from_entropy(
-    float density, float entropy, const struct Til_params *mat) {
+    const float density, const float entropy, const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -152,7 +269,8 @@ INLINE static float Til_internal_energy_from_entropy(
 }
 
 // gas_pressure_from_entropy
-INLINE static float Til_pressure_from_entropy(float density, float entropy,
+INLINE static float Til_pressure_from_entropy(const float density,
+                                              const float entropy,
                                               const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -161,7 +279,8 @@ INLINE static float Til_pressure_from_entropy(float density, float entropy,
 }
 
 // gas_entropy_from_pressure
-INLINE static float Til_entropy_from_pressure(float density, float pressure,
+INLINE static float Til_entropy_from_pressure(const float density,
+                                              const float pressure,
                                               const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -170,7 +289,8 @@ INLINE static float Til_entropy_from_pressure(float density, float pressure,
 }
 
 // gas_soundspeed_from_entropy
-INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
+INLINE static float Til_soundspeed_from_entropy(const float density,
+                                                const float entropy,
                                                 const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -180,14 +300,14 @@ INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
 
 // gas_entropy_from_internal_energy
 INLINE static float Til_entropy_from_internal_energy(
-    float density, float u, const struct Til_params *mat) {
+    const float density, const float u, const struct Til_params *mat) {
 
   return 0.f;
 }
 
 // gas_pressure_from_internal_energy
 INLINE static float Til_pressure_from_internal_energy(
-    float density, float u, const struct Til_params *mat) {
+    const float density, const float u, const struct Til_params *mat) {
 
   const float eta = density / mat->rho_0;
   const float eta_sq = eta * eta;
@@ -195,10 +315,9 @@ INLINE static float Til_pressure_from_internal_energy(
   const float nu = 1.f / eta - 1.f;
   const float w = u / (mat->u_0 * eta_sq) + 1.f;
   const float w_inv = 1.f / w;
-  float P_c, P_e, P;
 
   // Condensed or cold
-  P_c =
+  float P_c =
       (mat->a + mat->b * w_inv) * density * u + mat->A * mu + mat->B * mu * mu;
   if (eta < mat->eta_zero) {
     P_c = 0.f;
@@ -206,11 +325,12 @@ INLINE static float Til_pressure_from_internal_energy(
     P_c *= (eta - mat->eta_zero) / (mat->eta_min - mat->eta_zero);
   }
   // Expanded and hot
-  P_e = mat->a * density * u +
-        (mat->b * density * u * w_inv + mat->A * mu * expf(-mat->beta * nu)) *
-            expf(-mat->alpha * nu * nu);
+  float P_e = mat->a * density * u + (mat->b * density * u * w_inv +
+                                      mat->A * mu * expf(-mat->beta * nu)) *
+                                         expf(-mat->alpha * nu * nu);
 
   // Condensed or cold state
+  float P;
   if ((1.f < eta) || (u < mat->u_iv)) {
     P = P_c;
   }
@@ -234,7 +354,7 @@ INLINE static float Til_pressure_from_internal_energy(
 
 // gas_internal_energy_from_pressure
 INLINE static float Til_internal_energy_from_pressure(
-    float density, float P, const struct Til_params *mat) {
+    const float density, const float P, const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
 
@@ -243,7 +363,7 @@ INLINE static float Til_internal_energy_from_pressure(
 
 // gas_soundspeed_from_internal_energy
 INLINE static float Til_soundspeed_from_internal_energy(
-    float density, float u, const struct Til_params *mat) {
+    const float density, const float u, const struct Til_params *mat) {
 
   const float rho_0_inv = 1.f / mat->rho_0;
   const float eta = density * rho_0_inv;
@@ -256,25 +376,25 @@ INLINE static float Til_soundspeed_from_internal_energy(
   const float w_inv_sq = w_inv * w_inv;
   const float exp_beta = expf(-mat->beta * nu);
   const float exp_alpha = expf(-mat->alpha * nu * nu);
-  float P_c, P_e, c_sq_c, c_sq_e, c_sq;
 
   // Condensed or cold
-  P_c =
+  float P_c =
       (mat->a + mat->b * w_inv) * density * u + mat->A * mu + mat->B * mu * mu;
   if (eta < mat->eta_zero) {
     P_c = 0.f;
   } else if (eta < mat->eta_min) {
     P_c *= (eta - mat->eta_zero) / (mat->eta_min - mat->eta_zero);
   }
-  c_sq_c = P_c * rho_inv * (1.f + mat->a + mat->b * w_inv) +
-           mat->b * (w - 1.f) * w_inv_sq * (2.f * u - P_c * rho_inv) +
-           rho_inv * (mat->A + mat->B * (eta_sq - 1.f));
+  float c_sq_c = P_c * rho_inv * (1.f + mat->a + mat->b * w_inv) +
+                 mat->b * (w - 1.f) * w_inv_sq * (2.f * u - P_c * rho_inv) +
+                 rho_inv * (mat->A + mat->B * (eta_sq - 1.f));
 
   // Expanded and hot
-  P_e = mat->a * density * u +
-        (mat->b * density * u * w_inv + mat->A * mu * exp_beta) * exp_alpha;
+  float P_e =
+      mat->a * density * u +
+      (mat->b * density * u * w_inv + mat->A * mu * exp_beta) * exp_alpha;
 
-  c_sq_e =
+  float c_sq_e =
       P_e * rho_inv * (1.f + mat->a + mat->b * w_inv * exp_alpha) +
       (mat->b * density * u * w_inv_sq / eta_sq *
            (rho_inv / mat->u_0 * (2.f * u - P_e * rho_inv) +
@@ -285,6 +405,7 @@ INLINE static float Til_soundspeed_from_internal_energy(
           exp_alpha;
 
   // Condensed or cold state
+  float c_sq;
   if ((1.f < eta) || (u < mat->u_iv)) {
     c_sq = c_sq_c;
   }
@@ -304,7 +425,8 @@ INLINE static float Til_soundspeed_from_internal_energy(
 }
 
 // gas_soundspeed_from_pressure
-INLINE static float Til_soundspeed_from_pressure(float density, float P,
+INLINE static float Til_soundspeed_from_pressure(const float density,
+                                                 const float P,
                                                  const struct Til_params *mat) {
 
   error("This EOS function is not yet implemented!");
@@ -312,4 +434,494 @@ INLINE static float Til_soundspeed_from_pressure(float density, float P,
   return 0.f;
 }
 
+// Compute u cold
+INLINE static float compute_u_cold(
+    const float density, const struct Til_params *mat,
+    const enum eos_planetary_material_id mat_id) {
+  const int N = 10000;
+  const float rho_0 = mat->rho_0;
+  const float drho = (density - rho_0) / N;
+
+  float x = rho_0;
+  float u_cold = 1e-9;
+  for (int i = 0; i < N; i++) {
+    x += drho;
+    u_cold +=
+        Til_pressure_from_internal_energy(x, u_cold, mat) * drho / (x * x);
+  }
+
+  return u_cold;
+}
+
+// Compute A1_u_cold
+INLINE static void set_Til_u_cold(struct Til_params *mat,
+                                  enum eos_planetary_material_id mat_id) {
+
+  const int N = 10000;
+  const float rho_min = 100.f;
+  const float rho_max = 100000.f;
+
+  // Allocate table memory
+  mat->A1_u_cold = (float *)malloc(N * sizeof(float));
+
+  float rho = rho_min;
+  const float drho = (rho_max - rho_min) / (N - 1);
+
+  for (int i = 0; i < N; i++) {
+    mat->A1_u_cold[i] = compute_u_cold(rho, mat, mat_id);
+    rho += drho;
+  }
+}
+
+// Compute u cold fast from precomputed values
+INLINE static float compute_fast_u_cold(const float density,
+                                        const struct Til_params *mat) {
+
+  const int N = 10000;
+  const float rho_min = mat->rho_cold_min;
+  const float rho_max = mat->rho_cold_max;
+
+  const float drho = (rho_max - rho_min) / (N - 1);
+
+  const int a = (int)((density - rho_min) / drho);
+  const int b = a + 1;
+
+  float u_cold;
+  if (a >= 0 && a < (N - 1)) {
+    u_cold = mat->A1_u_cold[a];
+    u_cold += ((mat->A1_u_cold[b] - mat->A1_u_cold[a]) / drho) *
+              (density - rho_min - a * drho);
+  } else if (density < rho_min) {
+    u_cold = mat->A1_u_cold[0];
+  } else {
+    u_cold = mat->A1_u_cold[N - 1];
+    u_cold += ((mat->A1_u_cold[N - 1] - mat->A1_u_cold[N - 2]) / drho) *
+              (density - rho_max);
+  }
+  return u_cold;
+}
+
+// gas_temperature_from_internal_energy
+INLINE static float Til_temperature_from_internal_energy(
+    const float density, const float u, const struct Til_params *mat) {
+
+  const float u_cold = compute_fast_u_cold(density, mat);
+
+  float T = (u - u_cold) / (mat->C_V);
+  if (T < 0.f) {
+    T = 0.f;
+  }
+
+  return T;
+}
+
+// gas_pressure_from_density_and_temperature
+INLINE static float Til_pressure_from_temperature(
+    const float density, const float T, const struct Til_params *mat) {
+
+  const float u = compute_fast_u_cold(density, mat) + mat->C_V * T;
+  const float P = Til_pressure_from_internal_energy(density, u, mat);
+
+  return P;
+}
+
+// gas_density_from_pressure_and_temperature
+INLINE static float Til_density_from_pressure_and_temperature(
+    const float P, const float T, const struct Til_params *mat) {
+
+  float rho_min = mat->rho_min;
+  float rho_max = mat->rho_max;
+  float rho_mid = (rho_min + rho_max) / 2.f;
+
+  // Check for P == 0 or T == 0
+  float P_des;
+  if (P <= mat->P_min) {
+    P_des = mat->P_min;
+  } else {
+    P_des = P;
+  }
+
+  float P_min = Til_pressure_from_temperature(rho_min, T, mat);
+  float P_mid = Til_pressure_from_temperature(rho_mid, T, mat);
+  float P_max = Til_pressure_from_temperature(rho_max, T, mat);
+
+  // quick fix?
+  if (P_des < P_min) {
+    P_des = P_min;
+  }
+
+  const float tolerance = 0.001 * rho_min;
+  int counter = 0;
+  const int max_counter = 200;
+  if (P_des >= P_min && P_des <= P_max) {
+    while ((rho_max - rho_min) > tolerance && counter < max_counter) {
+
+      P_min = Til_pressure_from_temperature(rho_min, T, mat);
+      P_mid = Til_pressure_from_temperature(rho_mid, T, mat);
+      P_max = Til_pressure_from_temperature(rho_max, T, mat);
+
+      const float f0 = P_des - P_min;
+      const float f2 = P_des - P_mid;
+
+      if ((f0 * f2) > 0) {
+        rho_min = rho_mid;
+      } else {
+        rho_max = rho_mid;
+      }
+
+      rho_mid = (rho_min + rho_max) / 2.f;
+      counter += 1;
+    }
+  } else {
+    error("Error in Til_density_from_pressure_and_temperature");
+    return 0.f;
+  }
+  return rho_mid;
+}
+
+// gas_density_from_pressure_and_internal_energy
+INLINE static float Til_density_from_pressure_and_internal_energy(
+    const float P, const float u, const float rho_ref, const float rho_sph,
+    const struct Til_params *mat) {
+
+  if (P <= mat->P_min || u == 0) {
+    return rho_sph;
+  }
+
+  // We start search on the same curve as rho_ref, since this is most likely
+  // curve to find rho on
+  float eta_ref = rho_ref / mat->rho_0;
+
+  /*
+  There are 5 possible curves:
+    1: cold_min
+    2: cold
+    3: hybrid_min
+    4: hybrid
+    5: hot
+
+  These curves cover different eta ranges within three different regions of u:
+  u REGION 1 (u < u_iv):
+      eta < eta_min:           cold_min
+      eta_min < eta:           cold
+
+  u REGION 2 (u_iv < u < u_cv):
+      eta < eta_min:           hybrid_min
+      eta_min < eta < 1:       hybrid
+      1 < eta:                 cold
+
+  u REGION 3 (u_cv < u):
+      eta < 1:                 hot
+      1 < eta:                 cold
+
+  NOTE: for a lot of EoS, eta_min = 0, so search this region last if given the
+  option to save time for most EoS
+  */
+  enum Til_region {
+    Til_region_none,
+    Til_region_cold_min,
+    Til_region_cold,
+    Til_region_hybrid_min,
+    Til_region_hybrid,
+    Til_region_hot
+  };
+
+  // Based on our u region, what possible curves can rho be on? Ordered based on
+  // order we search for roots. Numbers correspond to curves in order given
+  // above. 0 is a dummy which breaks the loop.
+  enum Til_region possible_curves[3];
+  // u REGION 1
+  if (u <= mat->u_iv) {
+    if (eta_ref <= mat->eta_min) {
+      possible_curves[0] = Til_region_cold_min;
+      possible_curves[1] = Til_region_cold;
+      possible_curves[2] = Til_region_none;
+    } else {
+      possible_curves[0] = Til_region_cold;
+      possible_curves[1] = Til_region_cold_min;
+      possible_curves[2] = Til_region_none;
+    }
+    // u REGION 2
+  } else if (u <= mat->u_cv) {
+    if (eta_ref <= mat->eta_min) {
+      possible_curves[0] = Til_region_hybrid_min;
+      possible_curves[1] = Til_region_hybrid;
+      possible_curves[2] = Til_region_cold;
+    } else if (eta_ref <= 1) {
+      possible_curves[0] = Til_region_hybrid;
+      possible_curves[1] = Til_region_cold;
+      possible_curves[2] = Til_region_hybrid_min;
+    } else {
+      possible_curves[0] = Til_region_cold;
+      possible_curves[1] = Til_region_hybrid;
+      possible_curves[2] = Til_region_hybrid_min;
+    }
+    // u REGION 3
+  } else {
+    if (eta_ref <= 1) {
+      possible_curves[0] = Til_region_hot;
+      possible_curves[1] = Til_region_cold;
+      possible_curves[2] = Til_region_none;
+    } else {
+      possible_curves[0] = Til_region_cold;
+      possible_curves[1] = Til_region_hot;
+      possible_curves[2] = Til_region_none;
+    }
+  }
+  // Newton-Raphson
+  const int max_iter = 10;
+  const float tol = 1e-5;
+
+  // loops over possible curves
+  for (int i = 0; i < 3; i++) {
+
+    enum Til_region curve = possible_curves[i];
+
+    // if there are only two possible curves, break when we get to three and
+    // haven't found a root
+    if (curve == Til_region_none) {
+      break;
+    }
+
+    // Start iteration at reference value
+    float rho_iter = rho_ref;
+
+    // Constrain our initial guess to be on the curve we're currently looking at
+    // in the first loop, this is already satisfied.
+    if (i > 0) {
+
+      // u REGION 1
+      if (u <= mat->u_iv) {
+        if (curve == Til_region_cold_min) {
+          if (rho_iter > mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+        // u REGION 2
+      } else if (u <= mat->u_cv) {
+        if (curve == Til_region_hybrid_min) {
+          if (rho_iter > mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else if (curve == Til_region_hybrid) {
+          if (rho_iter < mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+          if (rho_iter > mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+
+        // u REGION 3
+      } else {
+        if (curve == Til_region_hot) {
+          if (rho_iter > mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+      }
+    }
+
+    // Set this to an arbitrary number so we definitely dont't think we converge
+    // straigt away
+    float last_rho_iter = -1e5;
+
+    // Now iterate
+    for (int j = 0; j < max_iter; j++) {
+
+      const float eta_iter = rho_iter / mat->rho_0;
+      const float eta_iter_sq = eta_iter * eta_iter;
+      const float mu_iter = eta_iter - 1.0f;
+      const float nu_iter = 1.0f / eta_iter - 1.0f;
+      const float w_iter = u / (mat->u_0 * eta_iter_sq) + 1.0f;
+      const float w_iter_inv = 1.0f / w_iter;
+      const float exp1 = expf(-mat->beta * nu_iter);
+      const float exp2 = expf(-mat->alpha * nu_iter * nu_iter);
+
+      // Derivatives
+      const float dw_inv_drho_iter =
+          (2.f * mat->u_0 * u * eta_iter / mat->rho_0) /
+          ((u + mat->u_0 * eta_iter_sq) * (u + mat->u_0 * eta_iter_sq));
+      const float dmu_drho_iter = 1.f / mat->rho_0;
+      const float dmu_sq_drho_iter =
+          2.f * rho_iter / (mat->rho_0 * mat->rho_0) - 2.f / mat->rho_0;
+      const float dexp1_drho_iter =
+          mat->beta * mat->rho_0 * exp1 / (rho_iter * rho_iter);
+      const float dexp2_drho_iter = 2.f * mat->alpha * mat->rho_0 *
+                                    (mat->rho_0 - rho_iter) * exp2 /
+                                    (rho_iter * rho_iter * rho_iter);
+
+      // Use P_fraction to determine whether we've converged on a root
+      float P_fraction;
+
+      // Newton-Raphson
+      float P_c_iter, dP_c_drho_iter, P_h_iter, dP_h_drho_iter;
+
+      // if "cold" or "hybrid"
+      if (curve == Til_region_cold || curve == Til_region_hybrid) {
+        P_c_iter = (mat->a + mat->b * w_iter_inv) * rho_iter * u +
+                   mat->A * mu_iter + mat->B * mu_iter * mu_iter - P;
+        dP_c_drho_iter = (mat->a + mat->b * w_iter_inv) * u +
+                         mat->b * u * rho_iter * dw_inv_drho_iter +
+                         mat->A * dmu_drho_iter + mat->B * dmu_sq_drho_iter;
+        P_fraction = P_c_iter / P;
+
+        // if curve is cold then we've got everything we need
+        if (curve == Til_region_cold) {
+          rho_iter -= P_c_iter / dP_c_drho_iter;
+          // Don't use these:
+          P_h_iter = 0.f;
+          dP_h_drho_iter = 0.f;
+        }
+        // if "cold_min" or "hybrid_min"
+        // Can only have one version of the cold curve, therefore either use the
+        // min version or the normal version hence "else if"
+      } else if (curve == Til_region_cold_min ||
+                 curve == Til_region_hybrid_min) {
+        P_c_iter = ((mat->a + mat->b * w_iter_inv) * rho_iter * u +
+                    mat->A * mu_iter + mat->B * mu_iter * mu_iter) *
+                       (eta_iter - mat->eta_zero) /
+                       (mat->eta_min - mat->eta_zero) -
+                   P;
+        dP_c_drho_iter =
+            ((mat->a + mat->b * w_iter_inv) * u +
+             mat->b * u * rho_iter * dw_inv_drho_iter + mat->A * dmu_drho_iter +
+             mat->B * dmu_sq_drho_iter) *
+                (eta_iter - mat->eta_zero) / (mat->eta_min - mat->eta_zero) +
+            ((mat->a + mat->b * w_iter_inv) * rho_iter * u + mat->A * mu_iter +
+             mat->B * mu_iter * mu_iter) *
+                (1 / (mat->rho_0 * (mat->eta_min - mat->eta_zero)));
+        P_fraction = P_c_iter / P;
+
+        // if curve is cold_min then we've got everything we need
+        if (curve == Til_region_cold_min) {
+          rho_iter -= P_c_iter / dP_c_drho_iter;
+          // Don't use these:
+          P_c_iter = 0.f;
+          dP_c_drho_iter = 0.f;
+        }
+      }
+
+      // if "hybrid_min" or "hybrid" or "hot"
+      if (curve == Til_region_hybrid_min || curve == Til_region_hybrid ||
+          curve == Til_region_hot) {
+        P_h_iter =
+            mat->a * rho_iter * u +
+            (mat->b * rho_iter * u * w_iter_inv + mat->A * mu_iter * exp1) *
+                exp2 -
+            P;
+        dP_h_drho_iter =
+            mat->a * u +
+            (mat->b * u * w_iter_inv +
+             mat->b * u * rho_iter * dw_inv_drho_iter +
+             mat->A * mu_iter * dexp1_drho_iter +
+             mat->A * exp1 * dmu_drho_iter) *
+                exp2 +
+            (mat->b * rho_iter * u * w_iter_inv + mat->A * mu_iter * exp1) *
+                dexp2_drho_iter;
+        P_fraction = P_h_iter / P;
+
+        // if curve is hot then we've got everything we need
+        if (curve == Til_region_hot) {
+          rho_iter -= P_h_iter / dP_h_drho_iter;
+        }
+      }
+
+      // If we are on a hybrid or hybrid_min curve, we combie hot and cold
+      // curves
+      if (curve == Til_region_hybrid_min || curve == Til_region_hybrid) {
+        const float P_hybrid_iter =
+            ((u - mat->u_iv) * P_h_iter + (mat->u_cv - u) * P_c_iter) /
+            (mat->u_cv - mat->u_iv);
+        const float dP_hybrid_drho_iter = ((u - mat->u_iv) * dP_h_drho_iter +
+                                           (mat->u_cv - u) * dP_c_drho_iter) /
+                                          (mat->u_cv - mat->u_iv);
+        rho_iter -= P_hybrid_iter / dP_hybrid_drho_iter;
+        P_fraction = P_hybrid_iter / P;
+      }
+
+      // Now we have to constrain the new rho_iter to the curve we're on
+      // u REGION 1
+      if (u <= mat->u_iv) {
+        if (curve == Til_region_cold_min) {
+          if (rho_iter > mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+        // u REGION 2
+      } else if (u <= mat->u_cv) {
+        if (curve == Til_region_hybrid_min) {
+          if (rho_iter > mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+        } else if (curve == Til_region_hybrid) {
+          if (rho_iter < mat->eta_min * mat->rho_0) {
+            rho_iter = mat->eta_min * mat->rho_0;
+          }
+          if (rho_iter > mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+
+        // u REGION 3
+      } else {
+        if (curve == Til_region_hot) {
+          if (rho_iter > mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else if (curve == Til_region_cold) {
+          if (rho_iter < mat->rho_0) {
+            rho_iter = mat->rho_0;
+          }
+        } else {
+          error("Error in Til_density_from_pressure_and_internal_energy");
+        }
+      }
+
+      // Either we've converged ...
+      if (fabsf(P_fraction) < tol) {
+        return rho_iter;
+
+        // ... or we're stuck at the boundary ...
+      } else if (rho_iter == last_rho_iter) {
+        break;
+      }
+
+      // ... or we loop again
+      last_rho_iter = rho_iter;
+    }
+  }
+  return rho_sph;
+}
+
 #endif /* SWIFT_TILLOTSON_EQUATION_OF_STATE_H */

src/extra_io/EAGLE/extra.h

@@ -21,10 +21,11 @@
 
 #include "chemistry.h"
 #include "cooling.h"
+#include "cooling/PS2020/cooling_tables.h"
 #include "engine.h"
 #include "star_formation.h"
 
-#define xray_table_date_string 20230110
+#define xray_table_date_string 20240406
 
 #define xray_emission_N_temperature 46
 #define xray_emission_N_density 71
@@ -139,9 +140,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   hid_t dataset = H5Dopen(tempfile_id, "Date_String", H5P_DEFAULT);
   herr_t status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
                           H5P_DEFAULT, &datestring);
-  if (status < 0) printf("error reading the date string");
+  if (status < 0) error("error reading the date string");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
   if (datestring != xray_table_date_string)
     error(
         "The table and code version do not match, please use table version %i",
@@ -155,9 +156,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/Temperature_bins", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->Temperatures);
-  if (status < 0) printf("error reading temperatures");
+  if (status < 0) error("error reading temperatures");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Read density bins */
   if (posix_memalign((void **)&xrays->Densities, SWIFT_STRUCT_ALIGNMENT,
@@ -167,9 +168,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/Density_bins", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->Densities);
-  if (status < 0) printf("error reading densities");
+  if (status < 0) error("error reading densities");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Read Helium bins */
   if (posix_memalign((void **)&xrays->He_bins, SWIFT_STRUCT_ALIGNMENT,
@@ -179,9 +180,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/He_bins", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->He_bins);
-  if (status < 0) printf("error reading Helium massfractions");
+  if (status < 0) error("error reading Helium massfractions");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Read solar metallicity */
   if (posix_memalign((void **)&xrays->Log10_solar_metallicity,
@@ -192,9 +193,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/Solar_metallicities", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->Log10_solar_metallicity);
-  if (status < 0) printf("error reading solar metalicities");
+  if (status < 0) error("error reading solar metalicities");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Get Solar metallicities from log solar metallicities */
   if (posix_memalign((void **)&xrays->Solar_metallicity, SWIFT_STRUCT_ALIGNMENT,
@@ -212,9 +213,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/Redshift_bins", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->Redshifts);
-  if (status < 0) printf("error reading redshift bins");
+  if (status < 0) error("error reading redshift bins");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Read element mass */
   if (posix_memalign((void **)&xrays->element_mass, SWIFT_STRUCT_ALIGNMENT,
@@ -224,9 +225,9 @@ INLINE static void read_xray_header(struct xray_properties *xrays,
   dataset = H5Dopen(tempfile_id, "Bins/Element_masses", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->element_mass);
-  if (status < 0) printf("error reading element masses");
+  if (status < 0) error("error reading element masses");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 }
 
 /**
@@ -261,9 +262,9 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   herr_t status =
       H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
               xrays->emissivity_erosita_low_intrinsic_photons);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* erosita-high intrinsic photons */
   if (swift_memalign("xrays_table_erosita_high_photons",
@@ -280,9 +281,9 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   dataset = H5Dopen(file_id, "erosita-high/photons_intrinsic", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->emissivity_erosita_high_intrinsic_photons);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* ROSAT intrinsic photons */
   if (swift_memalign("xray_table_ROSAT_photons",
@@ -297,9 +298,9 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   dataset = H5Dopen(file_id, "ROSAT/photons_intrinsic", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->emissivity_ROSAT_intrinsic_photons);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   // erosita-low intrinsic energies
   if (swift_memalign("xrays_table_erosita_low_energies",
@@ -316,9 +317,9 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   dataset = H5Dopen(file_id, "erosita-low/energies_intrinsic", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->emissivity_erosita_low_intrinsic_energies);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* erosita-high intrinsic energies */
   if (swift_memalign(
@@ -336,9 +337,9 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   dataset = H5Dopen(file_id, "erosita-high/energies_intrinsic", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->emissivity_erosita_high_intrinsic_energies);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* ROSAT intrinsic energies */
   if (swift_memalign("xray_table_ROSAT_energies",
@@ -354,13 +355,13 @@ INLINE static void read_xray_table(struct xray_properties *xrays,
   dataset = H5Dopen(file_id, "ROSAT/energies_intrinsic", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    xrays->emissivity_ROSAT_intrinsic_energies);
-  if (status < 0) printf("error reading X-Ray table\n");
+  if (status < 0) error("error reading X-Ray table\n");
   status = H5Dclose(dataset);
-  if (status < 0) printf("error closing dataset");
+  if (status < 0) error("error closing dataset");
 
   /* Close file */
   status = H5Fclose(file_id);
-  if (status < 0) printf("error closing file");
+  if (status < 0) error("error closing file");
 }
 
 /**

src/extra_io/EAGLE/extra_io.h

@@ -63,15 +63,15 @@ INLINE static int extra_io_write_particles(const struct part *parts,
                                            struct io_props *list,
                                            const int with_cosmology) {
 
-  list[0] = io_make_output_field_convert_part(
+  list[0] = io_make_physical_output_field_convert_part(
       "XrayPhotonLuminosities", DOUBLE, 3, UNIT_CONV_PHOTONS_PER_TIME, 0.f,
-      parts, xparts, convert_part_Xray_photons,
+      parts, xparts, /*can convert to comoving=*/0, convert_part_Xray_photons,
       "Intrinsic X-ray photon luminosities in various bands. This is 0 for "
       "star-forming particles.");
 
-  list[1] = io_make_output_field_convert_part(
+  list[1] = io_make_physical_output_field_convert_part(
       "XrayLuminosities", DOUBLE, 3, UNIT_CONV_POWER, 0.f, parts, xparts,
-      convert_part_Xray_energies,
+      /*can convert to comoving=*/0, convert_part_Xray_energies,
       "Intrinsic X-ray luminosities in various bands. This is 0 for "
       "star-forming particles.");
 

src/feedback/AGORA/feedback.c

@@ -159,6 +159,9 @@ void feedback_will_do_feedback(
     const struct unit_system* us, const struct phys_const* phys_const,
     const integertime_t ti_current, const double time_base) {
 
+  /* quit if the birth_scale_factor or birth_time is negative */
+  if (sp->birth_scale_factor < 0.0 || sp->birth_time < 0.0) return;
+
   /* skip if the particle is idle for feedback (it already exploded) */
   if (sp->feedback_data.idle == 1) {
     sp->feedback_data.will_do_feedback = 0;
@@ -223,6 +226,10 @@ void feedback_will_do_feedback(
  */
 int feedback_is_active(const struct spart* sp, const struct engine* e) {
 
+  /* the particle is inactive if its birth_scale_factor or birth_time is
+   * negative */
+  if (sp->birth_scale_factor < 0.0 || sp->birth_time < 0.0) return 0;
+
   return sp->feedback_data.will_do_feedback;
 }
 
@@ -237,14 +244,18 @@ void feedback_init_spart(struct spart* sp) {
 }
 
 /**
- * @brief Reset the feedback field when the spart is not
- * in a correct state for feeedback_will_do_feedback.
+ * @brief Prepare the feedback fields after a star is born.
+ *
+ * This function is called in the functions sink_copy_properties_to_star() and
+ * star_formation_copy_properties().
  *
- * This function is called in the timestep task.
+ * @param sp The #spart to act upon.
+ * @param feedback_props The feedback perties to use.
+ * @param star_type The stellar particle type.
  */
 void feedback_init_after_star_formation(
-    struct spart* sp, const struct feedback_props* feedback_props) {
-  feedback_init_spart(sp);
+    struct spart* sp, const struct feedback_props* feedback_props,
+    const enum stellar_type star_type) {
 
   /* Zero the energy of supernovae */
   sp->feedback_data.energy_ejected = 0;
@@ -254,6 +265,10 @@ void feedback_init_after_star_formation(
 
   /* The particle is not idle */
   sp->feedback_data.idle = 0;
+
+  /* Give to the star its appropriate type: single star, continuous IMF star or
+     single population star */
+  sp->star_type = star_type;
 }
 
 /**

src/feedback/AGORA/feedback.h

@@ -24,6 +24,7 @@
 #include "feedback_properties.h"
 #include "hydro_properties.h"
 #include "part.h"
+#include "stars.h"
 #include "units.h"
 
 /**
@@ -59,7 +60,8 @@ __attribute__((always_inline)) INLINE static void feedback_reset_part(
 void feedback_init_spart(struct spart* sp);
 
 void feedback_init_after_star_formation(
-    struct spart* sp, const struct feedback_props* feedback_props);
+    struct spart* sp, const struct feedback_props* feedback_props,
+    enum stellar_type star_type);
 
 /**
  * @brief Should we do feedback for this star?

src/feedback/EAGLE_thermal/feedback.h

@@ -208,7 +208,8 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_feedback(
     const int with_cosmology) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
+  if (sp->birth_time == -1. && dt > 0.)
+    error("Evolving a star particle that should not!");
 #endif
 
   /* Start by finishing the loops over neighbours */
@@ -217,7 +218,8 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_feedback(
   const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
 
   sp->feedback_data.to_collect.ngb_rho *= h_inv_dim;
-  const float rho_inv = 1.f / sp->feedback_data.to_collect.ngb_rho;
+  const float rho = sp->feedback_data.to_collect.ngb_rho;
+  const float rho_inv = rho != 0.f ? 1.f / rho : 0.f;
   sp->feedback_data.to_collect.ngb_Z *= h_inv_dim * rho_inv;
 
   /* Compute amount of enrichment and feedback that needs to be done in this