diff --git a/examples/test.c b/examples/test.c
index 91c0242df6bcb7fd21677b470ce90de369f57967..f46d1fdff498beb1065725817c4b2bf97924a401 100644
--- a/examples/test.c
+++ b/examples/test.c
@@ -714,6 +714,7 @@ int main ( int argc , char *argv[] ) {
struct part *parts = NULL;
struct space s;
struct engine e;
+ struct UnitSystem us;
char ICfileName[200];
float dt_max = 0.0f;
ticks tic;
@@ -826,6 +827,18 @@ int main ( int argc , char *argv[] ) {
if ( myrank == 0 )
message( "sizeof(struct part) is %li bytes." , (long int)sizeof( struct part ));
+ /* Initilaize unit system */
+ initUnitSystem(&us);
+ if ( myrank == 0 )
+ {
+ message( "Unit system: U_L = %e cm.", us.UnitLength_in_cgs );
+ message( "Unit system: U_M = %e g.", us.UnitMass_in_cgs );
+ message( "Unit system: U_t = %e s.", us.UnitTime_in_cgs );
+ message( "Unit system: U_I = %e A.", us.UnitCurrent_in_cgs );
+ message( "Unit system: U_T = %e K.", us.UnitTemperature_in_cgs );
+ /* message( "Density units: %e a^%f h^%f.", conversionFactor(&us, UNIT_CONV_ENTROPY), aFactor(&us, UNIT_CONV_ENTROPY), hFactor(&us, UNIT_CONV_ENTROPY) ); */
+ }
+
/* Read particles and space information from (GADGET) IC */
tic = getticks();
#ifdef WITH_MPI
@@ -840,7 +853,7 @@ int main ( int argc , char *argv[] ) {
if( scaling != 1.0 )
for ( k = 0 ; k < N ; k++ )
parts[k].h *= scaling;
-
+
/* Apply shift */
if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )
for ( k = 0 ; k < N ; k++ ) {
diff --git a/src/Makefile.am b/src/Makefile.am
index 89738d8f7516ff95147c2fcb431a74ec0b28d6e4..4fad1b22472eea715110f363cb314504f02fa2e1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,7 +42,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
- io.c timers.c debug.c scheduler.c proxy.c parallel_io.c
+ io.c timers.c debug.c scheduler.c proxy.c parallel_io.c units.c
# Sources and flags for regular library
libswiftsim_la_SOURCES = $(AM_SOURCES)
diff --git a/src/const.h b/src/const.h
index aadd66e7129c7907f452a4ae4998120f4c332ecd..70271a4775fb2c6fcb86fceeb6c27a2ae0797a51 100644
--- a/src/const.h
+++ b/src/const.h
@@ -43,3 +43,9 @@
/* SPH variant to use */
#define LEGACY_GADGET2_SPH
+
+
+/* System of units */
+#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
+#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
+#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
diff --git a/src/swift.h b/src/swift.h
index f202d57e37daaa85225f8128c24fc3331c4ca0f4..02cf8333be1a212ba5dd5cb140dfc2ef6359a0be 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -39,6 +39,7 @@
#include "io.h"
#include "parallel_io.h"
#include "debug.h"
+#include "units.h"
#ifdef LEGACY_GADGET2_SPH
#include "runner_iact_legacy.h"
diff --git a/src/units.c b/src/units.c
new file mode 100644
index 0000000000000000000000000000000000000000..0b354473bb3719f6429e74554982937da75553af
--- /dev/null
+++ b/src/units.c
@@ -0,0 +1,194 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Matthieu Schaller (matthieu.schaller@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"
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stddef.h>
+#include <math.h>
+
+#include "const.h"
+#include "cycle.h"
+#include "part.h"
+#include "error.h"
+#include "units.h"
+
+
+/**
+ * @brief Initialises the UnitSystem structure with the constants given in const.h
+ * @param us The UnitSystem to initialize
+ */
+
+void initUnitSystem(struct UnitSystem* us)
+{
+ us->UnitMass_in_cgs = const_unit_mass_in_cgs;
+ us->UnitLength_in_cgs = const_unit_length_in_cgs;
+ us->UnitTime_in_cgs = 1.d / (const_unit_velocity_in_cgs / const_unit_length_in_cgs );
+ us->UnitCurrent_in_cgs = 1.;
+ us->UnitTemperature_in_cgs = 1.;
+}
+
+
+void getBaseUnitExponantsArray(float baseUnitsExp[5], enum UnitConversionFactor unit)
+{
+ switch( unit )
+ {
+ case UNIT_CONV_MASS:
+ baseUnitsExp[UNIT_MASS] = 1; break;
+
+ case UNIT_CONV_LENGTH:
+ baseUnitsExp[UNIT_LENGTH] = 1; break;
+
+ case UNIT_CONV_TIME:
+ baseUnitsExp[UNIT_TIME] = 1; break;
+
+ case UNIT_CONV_FREQUENCY:
+ baseUnitsExp[UNIT_TIME] = -1; break;
+
+ case UNIT_CONV_DENSITY:
+ baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = -3; break;
+
+ case UNIT_CONV_SPEED:
+ baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -1; break;
+
+ case UNIT_CONV_ACCELERATION:
+ baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ case UNIT_CONV_FORCE:
+ baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ case UNIT_CONV_ENERGY:
+ baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ case UNIT_CONV_ENERGY_PER_UNIT_MASS:
+ baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ case UNIT_CONV_ENTROPY:
+ baseUnitsExp[UNIT_MASS] = 1. - const_hydro_gamma; baseUnitsExp[UNIT_LENGTH] = 3.*const_hydro_gamma - 1.; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ case UNIT_CONV_POWER:
+ baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -3; break;
+
+ case UNIT_CONV_PRESSURE:
+ baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = -1; baseUnitsExp[UNIT_TIME] = -2; break;
+
+ }
+}
+
+
+/**
+ * @brief Returns the conversion factor for a given unit in the chosen unit system
+ * @param us The system of units in use
+ * @param unit The unit to convert
+ */
+double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
+{
+ float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
+
+ getBaseUnitExponantsArray(baseUnitsExp, unit);
+
+ return generalConversionFactor(us, baseUnitsExp);
+}
+
+/**
+ * @brief Returns the h factor exponentiation for a given unit
+ * @param us The system of units in use
+ * @param unit The unit to convert
+ */
+double hFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
+{
+ float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
+
+ getBaseUnitExponantsArray(baseUnitsExp, unit);
+
+ return generalhFactor(us, baseUnitsExp);
+
+}
+
+
+/**
+ * @brief Returns the scaling factor exponentiation for a given unit
+ * @param us The system of units in use
+ * @param unit The unit to convert
+ */
+double aFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
+{
+ float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
+
+ getBaseUnitExponantsArray(baseUnitsExp, unit);
+
+ return generalaFactor(us, baseUnitsExp);
+
+}
+
+
+
+/**
+ * @brief Returns the conversion factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
+ * @param us The unit system used
+ * @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
+ */
+double generalConversionFactor(struct UnitSystem* us, float baseUnitsExponants[5])
+{
+ double factor = 1.;
+ int i;
+
+ for(i = 0 ; i < 5 ; ++i )
+ if(baseUnitsExponants[i] != 0)
+ factor *= pow( *( &( us->UnitMass_in_cgs ) + i ) , baseUnitsExponants[i] );
+
+ return factor;
+}
+
+
+/**
+ * @brief Returns the h factor exponentiation for a given unit (expressed in terms of the 5 fundamental units)
+ * @param us The unit system used
+ * @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
+ */
+double generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5])
+{
+ float factor_exp = 0;
+
+ factor_exp += -baseUnitsExponants[UNIT_MASS];
+ factor_exp += -baseUnitsExponants[UNIT_LENGTH];
+ factor_exp += -baseUnitsExponants[UNIT_TIME];
+
+
+ return factor_exp;
+}
+
+/**
+ * @brief Returns the scaling factor exponentiation for a given unit (expressed in terms of the 5 fundamental units)
+ * @param us The unit system used
+ * @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
+ */
+double generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5])
+{
+ float factor_exp = 0;
+
+ factor_exp += baseUnitsExponants[UNIT_LENGTH];
+
+ return factor_exp;
+}
diff --git a/src/units.h b/src/units.h
new file mode 100644
index 0000000000000000000000000000000000000000..1216c396cce156764995bffbecb6a29df80a5e13
--- /dev/null
+++ b/src/units.h
@@ -0,0 +1,114 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@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/>.
+ *
+ ******************************************************************************/
+
+
+/**
+ * @brief The unit system used internally.
+ *
+ * This structure contains the conversion factors to the 7 cgs base units to the internal units.
+ * It is used everytime a conversion is performed or an i/o function is called.
+ *
+ **/
+struct UnitSystem
+{
+ double UnitMass_in_cgs; /*< Conversion factor from grams to internal mass units */
+
+ double UnitLength_in_cgs; /*< Conversion factor from centimeters to internal length units. */
+
+ double UnitTime_in_cgs; /*< Conversion factor from seconds to internal time units. */
+
+ double UnitCurrent_in_cgs; /*< Conversion factor from Ampere to internal current units. */
+
+ double UnitTemperature_in_cgs; /*< Conversion factor from Kelvins to internal temperature units. */
+};
+
+/**
+ * @brief The base units used in the cgs (and internal) system. All units are derived from those.
+ */
+enum BaseUnits
+ {
+ UNIT_MASS,
+ UNIT_LENGTH,
+ UNIT_TIME,
+ UNIT_CURRENT,
+ UNIT_TEMPERATURE
+ };
+
+
+/**
+ * @brief The different conversion factors supported by default
+ */
+enum UnitConversionFactor
+ {
+ UNIT_CONV_MASS,
+ UNIT_CONV_LENGTH,
+ UNIT_CONV_TIME,
+ UNIT_CONV_DENSITY,
+ UNIT_CONV_SPEED,
+ UNIT_CONV_ACCELERATION,
+ UNIT_CONV_FORCE,
+ UNIT_CONV_ENERGY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ UNIT_CONV_ENTROPY,
+ UNIT_CONV_POWER,
+ UNIT_CONV_PRESSURE,
+ UNIT_CONV_FREQUENCY
+ };
+
+
+/**
+ * @brief Initialises the UnitSystem structure with the constants given in const.h
+ */
+void initUnitSystem(struct UnitSystem*);
+
+
+/**
+ * @brief Returns the conversion factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
+ */
+double generalConversionFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
+
+
+/**
+ * @brief Returns the conversion factor for a given unit in the chosen unit system
+ */
+double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
+
+
+/**
+ * @brief Returns the h factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
+ */
+double generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
+
+
+/**
+ * @brief Returns the h factor for a given unit in the chosen unit system
+ */
+double hFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
+
+
+/**
+ * @brief Returns the scaling factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
+ */
+double generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
+
+
+/**
+ * @brief Returns the scaling factor for a given unit in the chosen unit system
+ */
+double aFactor(struct UnitSystem* us, enum UnitConversionFactor unit);