diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..a3b4714fe597b3f4098e6f335a0a43026f514ae5
--- /dev/null
+++ b/src/cooling/const_du/cooling.h
@@ -0,0 +1,139 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@durham.ac.uk)
+ *
+ * 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_COOLING_CONST_DU
+#define SWIFT_COOLING_CONST_DU
+
+/**
+ * @file src/cooling/const/cooling.h
+ * @brief Routines related to the "constant cooling" cooling function.
+ *
+ * This is the simplest possible cooling function. A constant cooling rate with
+ * a minimal energy floor is applied.
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_data {
+
+ /*! Cooling rate in internal units. du_dt = -cooling_rate */
+ float cooling_rate;
+
+ /*! Minimally allowed internal energy of the particles */
+ float min_energy;
+
+ /*! Constant multiplication factor for time-step criterion */
+ float cooling_tstep_mult;
+};
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+ const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct cooling_data* cooling, struct part* p, float dt) {
+
+ /* Get current internal energy (dt=0) */
+ const float u_old = hydro_get_internal_energy(p, 0.f);
+
+ /* Get cooling function properties */
+ const float du_dt = -cooling->cooling_rate;
+ const float u_floor = cooling->min_energy;
+
+ /* Constant cooling with a minimal floor */
+ float u_new;
+ if (u_old - du_dt * dt > u_floor) {
+ u_new = u_old + du_dt * dt;
+ } else {
+ u_new = u_floor;
+ }
+
+ /* Update the internal energy */
+ hydro_set_internal_energy(p, u_new);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double cooling_timestep(
+ const struct cooling_data* cooling,
+ const struct phys_const* const phys_const, const struct part* const p) {
+
+ const float cooling_rate = cooling->cooling_rate;
+ const float internal_energy =
+ hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
+ return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+INLINE void cooling_init(const struct swift_params* parameter_file,
+ const struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ struct cooling_data* cooling) {
+
+ cooling->cooling_rate = parser_get_param_double(parameter_file, "Cooling:cooling_rate");
+ cooling->min_energy =
+ parser_get_param_double(parameter_file, "Cooling:min_energy");
+ cooling->cooling_tstep_mult =
+ parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+INLINE void cooling_print(const struct cooling_data* cooling) {
+
+ message("Cooling function is 'Constant cooling' with rate %f and floor %f",
+ cooling->cooling_rate, cooling->min_energy);
+}
+
+#endif /* SWIFT_CONST_COOLING_H */
diff --git a/src/cooling/const_lambda/cooling.c b/src/cooling/const_lambda/cooling.c
new file mode 100644
index 0000000000000000000000000000000000000000..1daa5d52d5794583d1818ffa0d1f7acec83998b5
--- /dev/null
+++ b/src/cooling/const_lambda/cooling.c
@@ -0,0 +1,150 @@
+
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@durham.ac.uk)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "adiabatic_index.h"
+#include "cooling.h"
+#include "hydro.h"
+
+/**
+ * @brief Initialises the cooling properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param us The current internal system of units
+ * @param cooling The cooling properties to initialize
+ */
+void cooling_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us,
+ const struct phys_const* const phys_const,
+ struct cooling_data* cooling) {
+
+ cooling->lambda = parser_get_param_double(parameter_file, "Cooling:Lambda");
+ cooling->min_temperature = parser_get_param_double(parameter_file, "Cooling:minimum_temperature");
+ cooling->mean_molecular_weight = parser_get_param_double(parameter_file, "Cooling:mean_molecular_weight");
+ cooling->hydrogen_mass_abundance = parser_get_param_double(parameter_file, "Cooling:hydrogen_mass_abundance");
+ cooling->cooling_tstep_mult = parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
+
+ /*convert minimum temperature into minimum internal energy*/
+ float u_floor =
+ phys_const->const_boltzmann_k * cooling->min_temperature /
+ (hydro_gamma_minus_one * cooling->mean_molecular_weight * phys_const->const_proton_mass);
+ float u_floor_cgs =
+ u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+ cooling->min_internal_energy = u_floor;
+ cooling->min_internal_energy_cgs = u_floor_cgs;
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The cooling properties.
+ */
+void cooling_print(const struct cooling_data* cooling) {
+
+ message(
+ "Cooling properties are (lambda, min_temperature, "
+ "hydrogen_mass_abundance, mean_molecular_weight, tstep multiplier) %g %g "
+ "%g %g %g",
+ cooling->lambda,
+ cooling->min_temperature,
+ cooling->hydrogen_mass_abundance,
+ cooling->mean_molecular_weight,
+ cooling->cooling_tstep_mult);
+}
+
+void update_entropy(const struct phys_const* const phys_const,
+ const struct UnitSystem* us,
+ const struct cooling_data* cooling, struct part* p,
+ float dt) {
+
+ /*updates the entropy of a particle after integrating the cooling equation*/
+ float u_old;
+ float u_new;
+ float new_entropy;
+ // float old_entropy = p->entropy;
+ float rho = p->rho;
+
+ // u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1);
+ u_old =
+ hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
+ u_new = calculate_new_thermal_energy(u_old, rho, dt, cooling, phys_const, us);
+ new_entropy = u_new * pow_minus_gamma_minus_one(rho) * hydro_gamma_minus_one;
+ p->entropy = new_entropy;
+}
+
+/*This function integrates the cooling equation, given the initial
+ thermal energy, density and the timestep dt. Returns the final internal
+ energy*/
+
+float calculate_new_thermal_energy(float u_old, float rho, float dt,
+ const struct cooling_data* cooling,
+ const struct phys_const* const phys_const,
+ const struct UnitSystem* us) {
+#ifdef CONST_COOLING
+ /*du/dt = -lambda, independent of density*/
+ float du_dt = -cooling->const_cooling.lambda;
+ float u_floor = cooling->const_cooling.min_energy;
+ float u_new;
+ if (u_old - du_dt * dt > u_floor) {
+ u_new = u_old + du_dt * dt;
+ } else {
+ u_new = u_floor;
+ }
+#endif /*CONST_COOLING*/
+
+#ifdef CREASEY_COOLING
+ /* rho*du/dt = -lambda*n_H^2 */
+ float u_new;
+ float X_H = cooling->hydrogen_mass_abundance;
+ float lambda_cgs = cooling->lambda; // this is always in cgs
+ float u_floor_cgs = cooling->min_internal_energy_cgs;
+
+ /*convert from internal code units to cgs*/
+ float dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+ float rho_cgs = rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ float m_p_cgs = phys_const->const_proton_mass *
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ float n_H_cgs = X_H * rho_cgs / m_p_cgs;
+ float u_old_cgs =
+ u_old * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
+ float u_new_cgs;
+
+ if (u_old_cgs + du_dt_cgs * dt_cgs > u_floor_cgs) {
+ u_new_cgs = u_old_cgs + du_dt_cgs * dt_cgs;
+ } else {
+ u_new_cgs = u_floor_cgs;
+ }
+ /*convert back to internal code units when returning new internal energy*/
+
+ u_new = u_new_cgs /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+#endif /*CREASEY_COOLING*/
+ return u_new;
+}
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..d5d7fe136605bb2acd44ea162388afab81f568e1
--- /dev/null
+++ b/src/cooling/const_lambda/cooling.h
@@ -0,0 +1,200 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@durham.ac.uk)
+ *
+ * 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_COOLING_CONST_LAMBDA
+#define SWIFT_COOLING_CONST_LAMBDA
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* Cooling Properties */
+struct cooling_data {
+
+ /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
+ float lambda;
+
+ /*! Minimum temperature (in Kelvin) for all gas particles*/
+ float min_temperature;
+
+ /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H*/
+ float hydrogen_mass_abundance;
+
+ /* 'mu', used to convert min_temperature to min_internal energy*/
+ float mean_molecular_weight;
+
+ /*! Minimally allowed internal energy of the particles */
+ float min_energy;
+ float min_energy_cgs;
+
+ /*! Constant multiplication factor for time-step criterion */
+ float cooling_tstep_mult;
+};
+
+/**
+ * @brief Calculates du/dt in code units for a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data..
+ */
+__attribute__((always_inline)) INLINE static float cooling_rate(
+ const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct cooling_data* cooling, const struct part* p) {
+
+ /* Get particle properties */
+ /* Density */
+ const float rho = p->rho;
+ /* Get cooling function properties */
+ const float X_H = cooling->hydrogen_mass_abundance;
+ /* lambda should always be set in cgs units */
+ const float lambda_cgs = cooling->lambda;
+
+ /*convert from internal code units to cgs*/
+ const float rho_cgs = rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ const float m_p_cgs = phys_const->const_proton_mass *
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ const float n_H_cgs = X_H * rho_cgs / m_p_cgs;
+
+ /* Calculate du_dt */
+ const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
+
+ /* Convert du/dt back to internal code units */
+ const float du_dt = du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME)/
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+ return du_dt;
+}
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+ const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct cooling_data* cooling, struct part* p, float dt) {
+
+
+ /* Get current internal energy (dt=0) */
+ const float u_old = hydro_get_internal_energy(p, 0.f);
+
+ /* Internal energy floor */
+ const float u_floor = cooling->min_energy;
+
+ /* Calculate du_dt */
+ const float du_dt = cooling_rate(phys_const,us,cooling,p);
+ /* Intergrate cooling equation, but enforce energy floor */
+ float u_new;
+ if (u_old + du_dt * dt > u_floor) {
+ u_new = u_old + du_dt * dt;
+ } else {
+ u_new = u_floor;
+ }
+
+ /* Update the internal energy */
+ hydro_set_internal_energy(p, u_new);
+ //const float u_new_test = hydro_get_internal_energy(p, 0.f);
+ /* if (-(u_new_test - u_old)/u_old > 1.0e-6){ */
+ /* printf("Particle has successfully cooled: u_old = %g , du_dt = %g , dt = %g , du_dt*dt = %g, u_old + du_dt*dt = %g, u_new = %g\n",u_old,du_dt,dt,du_dt*dt,u_new,u_new_test); */
+ /* exit(-1); */
+ /* } */
+}
+
+
+/**
+ * @brief Computes the time-step due to cooling
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float cooling_timestep(
+ const struct cooling_data* cooling, const struct phys_const* const phys_const,
+ const struct UnitSystem* us, const struct part* const p) {
+
+ /* Get du_dt */
+ const float du_dt = cooling_rate(phys_const,us,cooling,p);
+
+ /* Get current internal energy (dt=0) */
+ const float u = hydro_get_internal_energy(p, 0.f);
+
+ return u / du_dt;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+INLINE void cooling_init(const struct swift_params* parameter_file,
+ const struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ struct cooling_data* cooling) {
+
+ cooling->lambda = parser_get_param_double(parameter_file, "Cooling:lambda");
+ cooling->min_temperature = parser_get_param_double(parameter_file, "Cooling:minimum_temperature");
+ cooling->hydrogen_mass_abundance = parser_get_param_double(parameter_file, "Cooling:hydrogen_mass_abundance");
+ cooling->mean_molecular_weight = parser_get_param_double(parameter_file, "Cooling:mean_molecular_weight");
+ cooling->cooling_tstep_mult = parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
+
+ /*convert minimum temperature into minimum internal energy*/
+ const float u_floor =
+ phys_const->const_boltzmann_k * cooling->min_temperature /
+ (hydro_gamma_minus_one * cooling->mean_molecular_weight * phys_const->const_proton_mass);
+ const float u_floor_cgs =
+ u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
+ cooling->min_energy = u_floor;
+ cooling->min_energy_cgs = u_floor_cgs;
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+INLINE void cooling_print(const struct cooling_data* cooling) {
+
+ message("Cooling function is 'Constant lambda' with (lambda,min_temperature,hydrogen_mass_abundance,mean_molecular_weight) = (%g,%g,%g,%g)",
+ cooling->lambda, cooling->min_temperature,cooling->hydrogen_mass_abundance,cooling->mean_molecular_weight);
+}
+
+#endif /* SWIFT_COOLING_CONST_LAMBDA */
+