diff --git a/configure.ac b/configure.ac
index 9199f5cca29a5932c17ca34eb96bc497876bd532..c94f5c23b9cc5868a199ad4d177053c6f94f639d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -765,7 +765,7 @@ esac
 #  External potential
 AC_ARG_WITH([ext-potential],
    [AS_HELP_STRING([--with-ext-potential=<pot>],
-      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@]
+      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@]
    )],
    [with_potential="$withval"],
    [with_potential="none"]
@@ -783,6 +783,9 @@ case "$with_potential" in
    disc-patch)
       AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
    ;; 
+   sine-wave)
+      AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
+   ;; 
    *)
       AC_MSG_ERROR([Unknown external potential: $with_potential])
    ;;
diff --git a/examples/SineWavePotential_1D/makeIC.py b/examples/SineWavePotential_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..321af0714219dfa0c7cbb3d80845881dcbb8416d
--- /dev/null
+++ b/examples/SineWavePotential_1D/makeIC.py
@@ -0,0 +1,99 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart = 100
+boxSize = 1.
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart):
+  coords[i,0] = (i+np.random.random())/numPart
+
+for i in range(numPart):
+  # reasonable mass estimate (actually not really good, but better than assuming
+  # a constant volume)
+  if i == 0:
+    V = 0.5*(coords[1,0]-coords[-1,0]+1.)
+  elif i == numPart-1:
+    V = 0.5*(coords[0,0]+1.-coords[-2,0])
+  else:
+    V = 0.5*(coords[i+1,0] - coords[i-1,0])
+  rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+  m[i] = rho[i]*V
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 1
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_1D/plotSolution.py b/examples/SineWavePotential_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..65e981e4648fe0fe5d1da6cf3e753fb8a34f0fb4
--- /dev/null
+++ b/examples/SineWavePotential_1D/plotSolution.py
@@ -0,0 +1,77 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+  print "Need to provide a filename argument!"
+  exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+ids_reverse = np.argsort(ids)
+
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+  iself = int(ids[i])
+  corr = 0.
+  im1 = iself-1
+  if im1 < 0:
+    im1 = len(P)-1
+    corr = 1.
+  ip1 = iself+1
+  if ip1 == len(P):
+    ip1 = 0
+    corr = 1.
+  idxp1 = ids_reverse[ip1]
+  idxm1 = ids_reverse[im1]
+  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_1D/run.sh b/examples/SineWavePotential_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+  echo "Generating initial conditions for the 1D SineWavePotential example..."
+  python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+  python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_1D/sineWavePotential.yml b/examples/SineWavePotential_1D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..1f98b47a0e3c72e77a2024aebf4dcd863ed02e9b
--- /dev/null
+++ b/examples/SineWavePotential_1D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.
+  UnitLength_in_cgs:   1.
+  UnitVelocity_in_cgs: 1.
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1.
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   10.   # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sineWavePotential  # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          1.          # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  sineWavePotential.hdf5       # The file to read
+ 
+# External potential parameters
+SineWavePotential:
+  amplitude: 10.
+  timestep_limit: 1.
diff --git a/examples/SineWavePotential_2D/makeIC.py b/examples/SineWavePotential_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..62ae89f8f52bff9c0db37cd537f286ab817da3fe
--- /dev/null
+++ b/examples/SineWavePotential_2D/makeIC.py
@@ -0,0 +1,95 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart_1D = 50
+boxSize = [1., 1.]
+numPart = numPart_1D**2
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart_1D):
+  for j in range(numPart_1D):
+    coords[numPart_1D*i+j,0] = (i+0.5)/numPart_1D
+    coords[numPart_1D*i+j,1] = (j+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+  rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+  m[i] = rho[i]*V
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 2
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_2D/plotSolution.py b/examples/SineWavePotential_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..ee02f59c404db87a790465d2786e6296525e36b0
--- /dev/null
+++ b/examples/SineWavePotential_2D/plotSolution.py
@@ -0,0 +1,83 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+  print "Need to provide a filename argument!"
+  exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = int(np.sqrt(len(P)))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+  iself = int(ids[i]/n1D)
+  jself = int(ids[i]-n1D*iself)
+  corr = 0.
+  im1 = iself-1
+  if im1 < 0:
+    im1 = n1D-1
+    corr = 1.
+  ip1 = iself+1
+  if ip1 == n1D:
+    ip1 = 0
+    corr = 1.
+  idxp1 = ids_reverse[ip1*n1D+jself]
+  idxm1 = ids_reverse[im1*n1D+jself]
+  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_2D/run.sh b/examples/SineWavePotential_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_2D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+  echo "Generating initial conditions for the 1D SineWavePotential example..."
+  python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+  python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_2D/sineWavePotential.yml b/examples/SineWavePotential_2D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..f5a565ce4db46676c1f1f340f00025decea19644
--- /dev/null
+++ b/examples/SineWavePotential_2D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.
+  UnitLength_in_cgs:   1.
+  UnitVelocity_in_cgs: 1.
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1.
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.   # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sineWavePotential  # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.1          # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  sineWavePotential.hdf5       # The file to read
+ 
+# External potential parameters
+SineWavePotential:
+  amplitude: 10.
+  timestep_limit: 1.
diff --git a/examples/SineWavePotential_3D/makeIC.py b/examples/SineWavePotential_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..4833ec1b055e27b63751136f0491e972fb9e492a
--- /dev/null
+++ b/examples/SineWavePotential_3D/makeIC.py
@@ -0,0 +1,106 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart_1D = 20
+boxSize = [1., 1.]
+numPart = numPart_1D**3
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart_1D):
+#  coords[i,0] = (i+np.random.random())/numPart
+  for j in range(numPart_1D):
+    for k in range(numPart_1D):
+      coords[numPart_1D**2*i+numPart_1D*j+k,0] = (i+0.5)/numPart_1D
+      coords[numPart_1D**2*i+numPart_1D*j+k,1] = (j+0.5)/numPart_1D
+      coords[numPart_1D**2*i+numPart_1D*j+k,2] = (k+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+  # reasonable mass estimate (actually not really good, but better than assuming
+  # a constant volume)
+#  if i == 0:
+#    V = 0.5*(coords[1,0]-coords[-1,0]+1.)
+#  elif i == numPart-1:
+#    V = 0.5*(coords[0,0]+1.-coords[-2,0])
+#  else:
+#    V = 0.5*(coords[i+1,0] - coords[i-1,0])
+  rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+  m[i] = rho[i]*V
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_3D/plotSolution.py b/examples/SineWavePotential_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ae5dcd2a52235812c3c11cfcf4d01f4fde647af
--- /dev/null
+++ b/examples/SineWavePotential_3D/plotSolution.py
@@ -0,0 +1,84 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+  print "Need to provide a filename argument!"
+  exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = int(np.cbrt(len(P)))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+  iself = int(ids[i]/n1D/n1D)
+  jself = int(int(ids[i]-n1D*iself)/n1D)
+  kself = int(ids[i]-n1D**2*iself-n1D*jself)
+  corr = 0.
+  im1 = iself-1
+  if im1 < 0:
+    im1 = n1D-1
+    corr = 1.
+  ip1 = iself+1
+  if ip1 == n1D:
+    ip1 = 0
+    corr = 1.
+  idxp1 = ids_reverse[ip1*n1D**2+jself*n1D+kself]
+  idxm1 = ids_reverse[im1*n1D**2+jself*n1D+kself]
+  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_3D/run.sh b/examples/SineWavePotential_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_3D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+  echo "Generating initial conditions for the 1D SineWavePotential example..."
+  python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+  python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_3D/sineWavePotential.yml b/examples/SineWavePotential_3D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..391383568a5a94bd492cb228da4a7d4d24db413f
--- /dev/null
+++ b/examples/SineWavePotential_3D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.
+  UnitLength_in_cgs:   1.
+  UnitVelocity_in_cgs: 1.
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1.
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.1   # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sineWavePotential  # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.01          # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  sineWavePotential.hdf5       # The file to read
+ 
+# External potential parameters
+SineWavePotential:
+  amplitude: 10.
+  timestep_limit: 1.
diff --git a/src/Makefile.am b/src/Makefile.am
index f2d641637554fc71af61246e434b90fa334b8b32..29c110218416cc4bece0d766516eeb0a97fe4810 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -84,6 +84,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
 		 stars/Default/star_debug.h stars/Default/star_part.h  \
 	         potential/none/potential.h potential/point_mass/potential.h \
                  potential/isothermal/potential.h potential/disc_patch/potential.h \
+                 potential/sine_wave/potential.h \
 		 cooling/none/cooling.h cooling/none/cooling_struct.h \
 	         cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
                  cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
diff --git a/src/const.h b/src/const.h
index 09ec7636d1bcedc50af21d13ebd34f6eff940ee8..88c1a1af1cfc36cc401fdfea0b077f79fcd13bc0 100644
--- a/src/const.h
+++ b/src/const.h
@@ -55,6 +55,7 @@
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
 //#define GIZMO_FIX_PARTICLES
+//#define GIZMO_TOTAL_ENERGY
 
 /* Source terms */
 #define SOURCETERMS_NONE
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 71a082c6228591a64a4b71b950106670e7513c78..643489912a6c6b1db921e73b508910cc670d49ae 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -24,6 +24,7 @@
 #include "equation_of_state.h"
 #include "hydro_gradients.h"
 #include "minmax.h"
+#include "riemann.h"
 
 /**
  * @brief Computes the hydro time-step of a given particle
@@ -38,7 +39,15 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
 
   const float CFL_condition = hydro_properties->CFL_condition;
 
-  return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+  if (p->timestepvars.vmax == 0.) {
+    /* vmax can be zero in vacuum cells that only have vacuum neighbours */
+    /* in this case, the time step should be limited by the maximally
+       allowed time step. Since we do not know what that value is here, we set
+       the time step to a very large value */
+    return FLT_MAX;
+  } else {
+    return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+  }
 }
 
 /**
@@ -58,6 +67,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
 __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
     struct part* p, float dt) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (dt == 0.) {
+    error("Zero time step assigned to particle!");
+  }
+
+  if (dt != dt) {
+    error("NaN time step assigned to particle!");
+  }
+#endif
+
   p->force.dt = dt;
   p->force.active = 0;
 }
@@ -91,18 +110,33 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   p->conserved.momentum[1] = mass * p->primitives.v[1];
   p->conserved.momentum[2] = mass * p->primitives.v[2];
 
-  /* and the thermal energy */
+/* and the thermal energy */
+/* remember that we store the total thermal energy, not the specific thermal
+   energy (as in Gadget) */
+#if defined(EOS_ISOTHERMAL_GAS)
+  /* this overwrites the internal energy from the initial condition file */
+  p->conserved.energy = mass * const_isothermal_internal_energy;
+#else
   p->conserved.energy *= mass;
+#endif
+
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the total kinetic energy */
+  p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
+                                 p->conserved.momentum[1] * p->primitives.v[1] +
+                                 p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
 
 #if defined(GIZMO_FIX_PARTICLES)
+  /* make sure the particles are initially at rest */
   p->v[0] = 0.;
   p->v[1] = 0.;
   p->v[2] = 0.;
-#else
+#endif
+
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
-#endif
 }
 
 /**
@@ -190,6 +224,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   /* compute primitive variables */
   /* eqns (3)-(5) */
   const float m = p->conserved.mass;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (m == 0.) {
+    error("Mass is 0!");
+  }
+
+  if (volume == 0.) {
+    error("Volume is 0!");
+  }
+#endif
+
   float momentum[3];
   momentum[0] = p->conserved.momentum[0];
   momentum[1] = p->conserved.momentum[1];
@@ -198,10 +243,31 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->primitives.v[0] = momentum[0] / m;
   p->primitives.v[1] = momentum[1] / m;
   p->primitives.v[2] = momentum[2] / m;
-  const float energy = p->conserved.energy;
+
+#ifdef EOS_ISOTHERMAL_GAS
+  /* although the pressure is not formally used anywhere if an isothermal eos
+     has been selected, we still make sure it is set to the correct value */
+  p->primitives.P = const_isothermal_soundspeed * const_isothermal_soundspeed *
+                    p->primitives.rho;
+#else
+
+  float energy = p->conserved.energy;
+
+#ifdef GIZMO_TOTAL_ENERGY
+  /* subtract the kinetic energy; we want the thermal energy */
+  energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
+                    momentum[1] * p->primitives.v[1] +
+                    momentum[2] * p->primitives.v[2]);
+#endif
+
+  /* energy contains the total thermal energy, we want the specific energy.
+     this is why we divide by the volume, and not by the density */
   p->primitives.P = hydro_gamma_minus_one * energy / volume;
+#endif
 
   /* sanity checks */
+  /* it would probably be safer to throw a warning if netive densities or
+     pressures occur */
   if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
     p->primitives.rho = 0.0f;
     p->primitives.P = 0.0f;
@@ -231,6 +297,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->timestepvars.vmax = 0.0f;
 
   /* Set the actual velocity of the particle */
+  /* if GIZMO_FIX_PARTICLES has been selected, v_full will always be zero */
   p->force.v_full[0] = xp->v_full[0];
   p->force.v_full[1] = xp->v_full[1];
   p->force.v_full[2] = xp->v_full[2];
@@ -321,19 +388,33 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   else
     p->h *= expf(w1);
 
-  const float w2 = -hydro_dimension * w1;
-  if (fabsf(w2) < 0.2f) {
-    p->primitives.rho *= approx_expf(w2);
-  } else {
-    p->primitives.rho *= expf(w2);
+/* we temporarily disabled the primitive variable drift.
+   This should be reenabled once gravity works, and we have time to check that
+   drifting works properly. */
+//  const float w2 = -hydro_dimension * w1;
+//  if (fabsf(w2) < 0.2f) {
+//    p->primitives.rho *= approx_expf(w2);
+//  } else {
+//    p->primitives.rho *= expf(w2);
+//  }
+
+//  p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
+//  p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
+//  p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
+
+//#if !defined(EOS_ISOTHERMAL_GAS)
+//  if (p->conserved.mass > 0.) {
+//    const float u = p->conserved.energy + p->du_dt * dt;
+//    p->primitives.P =
+//        hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+//  }
+//#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (p->h <= 0.) {
+    error("Zero or negative smoothing length (%g)!", p->h);
   }
-
-  p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
-  p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
-  p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
-  const float u = p->conserved.energy + p->du_dt * dt;
-  p->primitives.P =
-      hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+#endif
 }
 
 /**
@@ -351,41 +432,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p) {
 
+  /* set the variables that are used to drift the primitive variables */
+
   /* Add normalization to h_dt. */
   p->force.h_dt *= p->h * hydro_dimension_inv;
 
-  /* Set the hydro acceleration, based on the new momentum and mass */
-  /* NOTE: the momentum and mass are only correct for active particles, since
-           only active particles have received flux contributions from all their
-           neighbours. Since this method is only called for active particles,
-           this is indeed the case. */
   if (p->force.dt) {
-    float mnew;
-    float vnew[3];
-
-    mnew = p->conserved.mass + p->conserved.flux.mass;
-    vnew[0] = (p->conserved.momentum[0] + p->conserved.flux.momentum[0]) / mnew;
-    vnew[1] = (p->conserved.momentum[1] + p->conserved.flux.momentum[1]) / mnew;
-    vnew[2] = (p->conserved.momentum[2] + p->conserved.flux.momentum[2]) / mnew;
-
-    p->a_hydro[0] = (vnew[0] - p->force.v_full[0]) / p->force.dt;
-    p->a_hydro[1] = (vnew[1] - p->force.v_full[1]) / p->force.dt;
-    p->a_hydro[2] = (vnew[2] - p->force.v_full[2]) / p->force.dt;
-
     p->du_dt = p->conserved.flux.energy / p->force.dt;
   } else {
-    p->a_hydro[0] = 0.0f;
-    p->a_hydro[1] = 0.0f;
-    p->a_hydro[2] = 0.0f;
-
     p->du_dt = 0.0f;
   }
 
 #if defined(GIZMO_FIX_PARTICLES)
-  p->a_hydro[0] = 0.0f;
-  p->a_hydro[1] = 0.0f;
-  p->a_hydro[2] = 0.0f;
-
   p->du_dt = 0.0f;
 
   /* disable the smoothing length update, since the smoothing lengths should
@@ -407,63 +465,54 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part* p, struct xpart* xp, float dt) {
 
-  float oldm, oldp[3], anew[3];
-  const float half_dt = 0.5f * dt;  // MATTHIEU
-
-  /* Retrieve the current value of the gravitational acceleration from the
-     gpart. We are only allowed to do this because this is the kick. We still
-     need to check whether gpart exists though.*/
-  if (p->gpart) {
-    anew[0] = p->gpart->a_grav[0];
-    anew[1] = p->gpart->a_grav[1];
-    anew[2] = p->gpart->a_grav[2];
-
-    /* Copy the old mass and momentum before updating the conserved variables */
-    oldm = p->conserved.mass;
-    oldp[0] = p->conserved.momentum[0];
-    oldp[1] = p->conserved.momentum[1];
-    oldp[2] = p->conserved.momentum[2];
-  }
+  float a_grav[3];
 
   /* Update conserved variables. */
   p->conserved.mass += p->conserved.flux.mass;
   p->conserved.momentum[0] += p->conserved.flux.momentum[0];
   p->conserved.momentum[1] += p->conserved.flux.momentum[1];
   p->conserved.momentum[2] += p->conserved.flux.momentum[2];
+#if defined(EOS_ISOTHERMAL_GAS)
+  p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
+#else
   p->conserved.energy += p->conserved.flux.energy;
+#endif
 
   /* Add gravity. We only do this if we have gravity activated. */
   if (p->gpart) {
-    p->conserved.momentum[0] +=
-        half_dt * (oldm * p->gravity.old_a[0] + p->conserved.mass * anew[0]);
-    p->conserved.momentum[1] +=
-        half_dt * (oldm * p->gravity.old_a[1] + p->conserved.mass * anew[1]);
-    p->conserved.momentum[2] +=
-        half_dt * (oldm * p->gravity.old_a[2] + p->conserved.mass * anew[2]);
-
-    float paold, panew;
-    paold = oldp[0] * p->gravity.old_a[0] + oldp[1] * p->gravity.old_a[1] +
-            oldp[2] * p->gravity.old_a[2];
-    panew = p->conserved.momentum[0] * anew[0] +
-            p->conserved.momentum[1] * anew[1] +
-            p->conserved.momentum[2] * anew[2];
-    p->conserved.energy += half_dt * (paold + panew);
-
-    float fluxaold, fluxanew;
-    fluxaold = p->gravity.old_a[0] * p->gravity.old_mflux[0] +
-               p->gravity.old_a[1] * p->gravity.old_mflux[1] +
-               p->gravity.old_a[2] * p->gravity.old_mflux[2];
-    fluxanew = anew[0] * p->gravity.mflux[0] + anew[1] * p->gravity.mflux[1] +
-               anew[2] * p->gravity.mflux[2];
-    p->conserved.energy += half_dt * (fluxaold + fluxanew);
-
-    /* Store gravitational acceleration and mass flux for next step */
-    p->gravity.old_a[0] = anew[0];
-    p->gravity.old_a[1] = anew[1];
-    p->gravity.old_a[2] = anew[2];
-    p->gravity.old_mflux[0] = p->gravity.mflux[0];
-    p->gravity.old_mflux[1] = p->gravity.mflux[1];
-    p->gravity.old_mflux[2] = p->gravity.mflux[2];
+    /* Retrieve the current value of the gravitational acceleration from the
+       gpart. We are only allowed to do this because this is the kick. We still
+       need to check whether gpart exists though.*/
+    a_grav[0] = p->gpart->a_grav[0];
+    a_grav[1] = p->gpart->a_grav[1];
+    a_grav[2] = p->gpart->a_grav[2];
+
+    /* Store the gravitational acceleration for later use. */
+    /* This is currently only used for output purposes. */
+    p->gravity.old_a[0] = a_grav[0];
+    p->gravity.old_a[1] = a_grav[1];
+    p->gravity.old_a[2] = a_grav[2];
+
+    /* Make sure the gpart knows the mass has changed. */
+    p->gpart->mass = p->conserved.mass;
+
+    /* Kick the momentum for half a time step */
+    /* Note that this also affects the particle movement, as the velocity for
+       the particles is set after this. */
+    p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
+    p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
+    p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
+
+#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
+    /* This part still needs to be tested! */
+    p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
+                                 p->conserved.momentum[1] * a_grav[1] +
+                                 p->conserved.momentum[2] * a_grav[2]);
+
+    p->conserved.energy += dt * (a_grav[0] * p->gravity.mflux[0] +
+                                 a_grav[1] * p->gravity.mflux[1] +
+                                 a_grav[2] * p->gravity.mflux[2]);
+#endif
   }
 
   /* reset fluxes */
@@ -474,6 +523,47 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.flux.momentum[1] = 0.0f;
   p->conserved.flux.momentum[2] = 0.0f;
   p->conserved.flux.energy = 0.0f;
+
+#if defined(GIZMO_FIX_PARTICLES)
+  xp->v_full[0] = 0.;
+  xp->v_full[1] = 0.;
+  xp->v_full[2] = 0.;
+
+  p->v[0] = 0.;
+  p->v[1] = 0.;
+  p->v[2] = 0.;
+
+  if (p->gpart) {
+    p->gpart->v_full[0] = 0.;
+    p->gpart->v_full[1] = 0.;
+    p->gpart->v_full[2] = 0.;
+  }
+#else
+  /* Set particle movement */
+  if (p->conserved.mass > 0.) {
+    xp->v_full[0] = p->conserved.momentum[0] / p->conserved.mass;
+    xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
+    xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
+  } else {
+    /* vacuum particles don't move */
+    xp->v_full[0] = 0.;
+    xp->v_full[1] = 0.;
+    xp->v_full[2] = 0.;
+  }
+  p->v[0] = xp->v_full[0];
+  p->v[1] = xp->v_full[1];
+  p->v[2] = xp->v_full[2];
+
+  /* Update gpart! */
+  /* This is essential, as the gpart drift is done independently from the part
+     drift, and we don't want the gpart and the part to have different
+     positions! */
+  if (p->gpart) {
+    p->gpart->v_full[0] = xp->v_full[0];
+    p->gpart->v_full[1] = xp->v_full[1];
+    p->gpart->v_full[2] = xp->v_full[2];
+  }
+#endif
 }
 
 /**
@@ -484,7 +574,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part* restrict p) {
 
-  return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+  if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+    return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#else
+    return const_isothermal_internal_energy;
+#endif
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -495,7 +593,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part* restrict p) {
 
-  return p->primitives.P / pow_gamma(p->primitives.rho);
+  if (p->primitives.rho > 0.) {
+    return p->primitives.P / pow_gamma(p->primitives.rho);
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -506,7 +608,16 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part* restrict p) {
 
-  return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+  if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+    return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+                 hydro_gamma_minus_one);
+#else
+    return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+#endif
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -558,6 +669,13 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
   /* conserved.energy is NOT the specific energy (u), but the total thermal
      energy (u*m) */
   p->conserved.energy = u * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the kinetic energy */
+  p->conserved.energy += 0.5f * p->conserved.mass *
+                         (p->conserved.momentum[0] * p->primitives.v[0] +
+                          p->conserved.momentum[1] * p->primitives.v[1] +
+                          p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
   p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
 }
 
@@ -573,7 +691,14 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
 __attribute__((always_inline)) INLINE static void hydro_set_entropy(
     struct part* restrict p, float S) {
 
-  p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
-                        p->conserved.mass;
-  p->primitives.P = gas_pressure_from_entropy(p->primitives.rho, S);
+  p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
+                        hydro_one_over_gamma_minus_one * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+  /* add the kinetic energy */
+  p->conserved.energy += 0.5f * p->conserved.mass *
+                         (p->conserved.momentum[0] * p->primitives.v[0] +
+                          p->conserved.momentum[1] * p->primitives.v[1] +
+                          p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+  p->primitives.P = S * pow_gamma(p->primitives.rho);
 }
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 65f4cb8cf892d1a9948bc043254ddb1581875ed3..a5c1e9038d0d3de6896afe773e3193a2304a6b6b 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -1,3 +1,4 @@
+
 /*******************************************************************************
  * This file is part of SWIFT.
  * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
@@ -21,6 +22,7 @@
 #define SWIFT_HYDRO_GRADIENTS_H
 
 #include "hydro_slope_limiters.h"
+#include "riemann.h"
 
 #if defined(GRADIENTS_SPH)
 
@@ -141,6 +143,33 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
 
   /* time */
   if (Wi[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+                             Wi[2] * pi->primitives.gradients.rho[1] +
+                             Wi[3] * pi->primitives.gradients.rho[2] +
+                             Wi[0] * (pi->primitives.gradients.v[0][0] +
+                                      pi->primitives.gradients.v[1][1] +
+                                      pi->primitives.gradients.v[2][2]));
+    dWi[1] -= 0.5 * mindt *
+              (Wi[1] * pi->primitives.gradients.v[0][0] +
+               Wi[2] * pi->primitives.gradients.v[0][1] +
+               Wi[3] * pi->primitives.gradients.v[0][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pi->primitives.gradients.rho[0] / Wi[0]);
+    dWi[2] -= 0.5 * mindt *
+              (Wi[1] * pi->primitives.gradients.v[1][0] +
+               Wi[2] * pi->primitives.gradients.v[1][1] +
+               Wi[3] * pi->primitives.gradients.v[1][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pi->primitives.gradients.rho[1] / Wi[0]);
+    dWi[3] -= 0.5 * mindt *
+              (Wi[1] * pi->primitives.gradients.v[2][0] +
+               Wi[2] * pi->primitives.gradients.v[2][1] +
+               Wi[3] * pi->primitives.gradients.v[2][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pi->primitives.gradients.rho[2] / Wi[0]);
+/* we don't care about P in this case */
+#else
     dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
                              Wi[2] * pi->primitives.gradients.rho[1] +
                              Wi[3] * pi->primitives.gradients.rho[2] +
@@ -166,9 +195,36 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
                                       pi->primitives.gradients.v[1][1] +
                                       pi->primitives.gradients.v[2][2]));
+#endif
   }
 
   if (Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+                             Wj[2] * pj->primitives.gradients.rho[1] +
+                             Wj[3] * pj->primitives.gradients.rho[2] +
+                             Wj[0] * (pj->primitives.gradients.v[0][0] +
+                                      pj->primitives.gradients.v[1][1] +
+                                      pj->primitives.gradients.v[2][2]));
+    dWj[1] -= 0.5 * mindt *
+              (Wj[1] * pj->primitives.gradients.v[0][0] +
+               Wj[2] * pj->primitives.gradients.v[0][1] +
+               Wj[3] * pj->primitives.gradients.v[0][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pj->primitives.gradients.rho[0] / Wj[0]);
+    dWj[2] -= 0.5 * mindt *
+              (Wj[1] * pj->primitives.gradients.v[1][0] +
+               Wj[2] * pj->primitives.gradients.v[1][1] +
+               Wj[3] * pj->primitives.gradients.v[1][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pj->primitives.gradients.rho[1] / Wj[0]);
+    dWj[3] -= 0.5 * mindt *
+              (Wj[1] * pj->primitives.gradients.v[2][0] +
+               Wj[2] * pj->primitives.gradients.v[2][1] +
+               Wj[3] * pj->primitives.gradients.v[2][2] +
+               const_isothermal_soundspeed * const_isothermal_soundspeed *
+                   pj->primitives.gradients.rho[2] / Wj[0]);
+#else
     dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
                              Wj[2] * pj->primitives.gradients.rho[1] +
                              Wj[3] * pj->primitives.gradients.rho[2] +
@@ -194,19 +250,36 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
                                       pj->primitives.gradients.v[1][1] +
                                       pj->primitives.gradients.v[2][2]));
+#endif
   }
 
-  Wi[0] += dWi[0];
+  if (-dWi[0] > Wi[0]) {
+    Wi[0] = 0.0f;
+  } else {
+    Wi[0] += dWi[0];
+  }
   Wi[1] += dWi[1];
   Wi[2] += dWi[2];
   Wi[3] += dWi[3];
-  Wi[4] += dWi[4];
+  if (-dWi[4] > Wi[4]) {
+    Wi[4] = 0.0f;
+  } else {
+    Wi[4] += dWi[4];
+  }
 
-  Wj[0] += dWj[0];
+  if (-dWj[0] > Wj[0]) {
+    Wj[0] = 0.0f;
+  } else {
+    Wj[0] += dWj[0];
+  }
   Wj[1] += dWj[1];
   Wj[2] += dWj[2];
   Wj[3] += dWj[3];
-  Wj[4] += dWj[4];
+  if (-dWj[4] > Wj[4]) {
+    Wj[4] = 0.0f;
+  } else {
+    Wj[4] += dWj[4];
+  }
 }
 
 #endif  // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 657119939d95c97ed049df95e9a4133fab135908..bdcea49e4fdc2a9ce17c43b99dab55d61f4095c8 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -232,8 +232,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* calculate the maximal signal velocity */
   if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    /* we use a value that is slightly higher than necessary, since the correct
+       value does not always work */
+    vmax = 2.5 * const_isothermal_soundspeed;
+#else
     vmax =
         sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+#endif
   } else {
     vmax = 0.0f;
   }
@@ -274,8 +280,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   float wi_dr = hidp1 * wi_dx;
   float wj_dr = hjdp1 * wj_dx;
   dvdr *= ri;
-  pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
-  if (mode == 1) {
+  if (pj->primitives.rho > 0.) {
+    pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+  }
+  if (mode == 1 && pi->primitives.rho > 0.) {
     pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
   }
 
@@ -391,6 +399,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   pi->conserved.flux.momentum[2] -= dti * Anorm * totflux[3];
   pi->conserved.flux.energy -= dti * Anorm * totflux[4];
 
+#ifndef GIZMO_TOTAL_ENERGY
   float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
                        pi->primitives.v[1] * pi->primitives.v[1] +
                        pi->primitives.v[2] * pi->primitives.v[2]);
@@ -398,6 +407,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
   pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
   pi->conserved.flux.energy -= dti * Anorm * totflux[0] * ekin;
+#endif
 
   /* here is how it works:
      Mode will only be 1 if both particles are ACTIVE and they are in the same
@@ -425,6 +435,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
     pj->conserved.flux.momentum[2] += dtj * Anorm * totflux[3];
     pj->conserved.flux.energy += dtj * Anorm * totflux[4];
 
+#ifndef GIZMO_TOTAL_ENERGY
     ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
                    pj->primitives.v[1] * pj->primitives.v[1] +
                    pj->primitives.v[2] * pj->primitives.v[2]);
@@ -432,6 +443,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
     pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
     pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
     pj->conserved.flux.energy += dtj * Anorm * totflux[0] * ekin;
+#endif
   }
 }
 
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 818217baee18c2128c8f004e3f2b14ecd454acfc..236106a1fb04cc2e5b84f997a2389d583ce17cff 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -70,7 +70,15 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
  * @return Internal energy of the particle
  */
 float convert_u(struct engine* e, struct part* p) {
-  return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+  if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+    return const_isothermal_internal_energy;
+#else
+    return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#endif
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -81,7 +89,11 @@ float convert_u(struct engine* e, struct part* p) {
  * @return Entropic function of the particle
  */
 float convert_A(struct engine* e, struct part* p) {
-  return p->primitives.P / pow_gamma(p->primitives.rho);
+  if (p->primitives.rho > 0.) {
+    return p->primitives.P / pow_gamma(p->primitives.rho);
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -92,6 +104,9 @@ float convert_A(struct engine* e, struct part* p) {
  * @return Total energy of the particle
  */
 float convert_Etot(struct engine* e, struct part* p) {
+#ifdef GIZMO_TOTAL_ENERGY
+  return p->conserved.energy;
+#else
   float momentum2;
 
   momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
@@ -99,6 +114,7 @@ float convert_Etot(struct engine* e, struct part* p) {
               p->conserved.momentum[2] * p->conserved.momentum[2];
 
   return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
+#endif
 }
 
 /**
@@ -111,7 +127,7 @@ float convert_Etot(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 13;
+  *num_fields = 14;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -142,6 +158,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
   list[12] =
       io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
                                         parts, conserved.energy, convert_Etot);
+  list[13] = io_make_output_field("GravAcceleration", FLOAT, 3,
+                                  UNIT_CONV_ACCELERATION, parts, gravity.old_a);
 }
 
 /**
diff --git a/src/multipole.h b/src/multipole.h
index 95b7450cd9c7aa60dba5c43ebd3c77e6316845a5..b5c9335ee8fabf46740cefe310fcfecbea3fd77e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -360,7 +360,7 @@ INLINE static void gravity_L2P(const struct gravity_tensors *l,
 #ifdef SWIFT_DEBUG_CHECKS
     struct gpart *gp = &gparts[i];
 
-// if(gpart_is_active(gp, e)){
+    // if(gpart_is_active(gp, e)){
 
     gp->mass_interacted += l->mass_interacted;
 #endif
diff --git a/src/potential.h b/src/potential.h
index 4c20f4c7b4eaa7ca134b3d90ae790a4710a1f217..e6ab9a5bd6bd91801eae0b3f1e3d8f65778f5065 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -36,6 +36,8 @@
 #include "./potential/isothermal/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
 #include "./potential/disc_patch/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)
+#include "./potential/sine_wave/potential.h"
 #else
 #error "Invalid choice of external potential"
 #endif
diff --git a/src/potential/sine_wave/potential.h b/src/potential/sine_wave/potential.h
new file mode 100644
index 0000000000000000000000000000000000000000..e2e2b8ffcc170c28a5facc8373a81746811a9991
--- /dev/null
+++ b/src/potential/sine_wave/potential.h
@@ -0,0 +1,133 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_SINE_WAVE_H
+#define SWIFT_SINE_WAVE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Sine wave case
+ */
+struct external_potential {
+
+  /*! Amplitude of the sine wave. */
+  double amplitude;
+
+  /*! Time-step limiting factor. */
+  double timestep_limit;
+};
+
+/**
+ * @brief Computes the time-step from the acceleration due to a sine wave.
+ *
+ * @param time The current time.
+ * @param potential The properties of the potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float external_gravity_timestep(
+    double time, const struct external_potential* restrict potential,
+    const struct phys_const* restrict phys_const,
+    const struct gpart* restrict g) {
+
+  return potential->timestep_limit;
+}
+
+/**
+ * @brief Computes the gravitational acceleration along x given by the sine
+ * wave.
+ *
+ * @param time The current time in internal units.
+ * @param potential The properties of the potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
+    double time, const struct external_potential* restrict potential,
+    const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+  g->a_grav[0] = potential->amplitude * sin(2. * M_PI * g->x[0]) /
+                 phys_const->const_newton_G;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in the
+ * sine wave.
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param gp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+    double time, const struct external_potential* potential,
+    const struct phys_const* const phys_const, const struct gpart* gp) {
+
+  /* this potential does not really have a potential energy */
+  return 0.;
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+    const struct swift_params* parameter_file,
+    const struct phys_const* phys_const, const struct unit_system* us,
+    const struct space* s, struct external_potential* potential) {
+
+  potential->amplitude =
+      parser_get_param_double(parameter_file, "SineWavePotential:amplitude");
+  potential->timestep_limit = parser_get_param_double(
+      parameter_file, "SineWavePotential:timestep_limit");
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param  potential The external potential properties.
+ */
+static INLINE void potential_print_backend(
+    const struct external_potential* potential) {
+
+  message("External potential is a sine wave with amplitude %g",
+          potential->amplitude);
+}
+
+#endif /* SWIFT_SINE_WAVE_H */
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
index 8e8a6a4dcce456b55f85148bed7342ad4e651c1b..26fb649f02ddd62c102e31c9191c20f9003c0133 100644
--- a/src/riemann/riemann_exact_isothermal.h
+++ b/src/riemann/riemann_exact_isothermal.h
@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
   /* Currently three possibilities and not really an algorithm to decide which
      one to choose: */
   /* just the average */
-  return 0.5f * (WL[0] + WR[0]);
+  //  return 0.5f * (WL[0] + WR[0]);
 
   /* two rarefaction approximation */
   return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
@@ -276,12 +276,24 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
   vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
   vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
 
-  /* VACUUM... */
+  if (WL[0] == 0. || WR[0] == 0.) {
+    error(
+        "One of the states is vacuum, the isothermal solver cannot solve "
+        "this!");
+  }
 
   rho = 0.;
   /* obtain a first guess for p */
   rhoguess = riemann_guess_rho(WL, WR, vL, vR);
-  frho = riemann_f(rho, WL, WR, vL, vR);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (rhoguess <= 0.) {
+    error("Zero or negative initial density guess.");
+  }
+#endif
+
+  /* we know that the value of riemann_f for rho=0 is negative (it is -inf). */
+  frho = -1.;
   frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
   /* ok, rhostar is close to 0, better use Brent's method... */
   /* we use Newton-Raphson until we find a suitable interval */
@@ -289,14 +301,18 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
     /* Newton-Raphson until convergence or until suitable interval is found
        to use Brent's method */
     unsigned int counter = 0;
-    while (fabs(rho - rhoguess) > 1.e-6f * 0.5f * (rho + rhoguess) &&
+    while (fabs(rho - rhoguess) > 5.e-7f * (rho + rhoguess) &&
            frhoguess < 0.0f) {
       rho = rhoguess;
       rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
       frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
       counter++;
       if (counter > 1000) {
-        error("Stuck in Newton-Raphson!\n");
+        error(
+            "Stuck in Newton-Raphson (rho: %g, rhoguess: %g, frhoguess: %g, "
+            "fprime: %g, rho-rhoguess: %g, WL: %g %g %g, WR: %g %g %g)!\n",
+            rho, rhoguess, frhoguess, riemann_fprime(rhoguess, WL, WR),
+            (rho - rhoguess), WL[0], vL, WL[4], WR[0], vR, WR[4]);
       }
     }
   }
@@ -304,7 +320,7 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
   if (1.e6 * fabs(rho - rhoguess) > 0.5f * (rho + rhoguess) &&
       frhoguess > 0.0f) {
     rho = 0.0f;
-    frho = riemann_f(rho, WL, WR, vL, vR);
+    frho = -1.;
     /* use Brent's method to find the zeropoint */
     rho = riemann_solve_brent(rho, rhoguess, frho, frhoguess, 1.e-6, WL, WR, vL,
                               vR);
diff --git a/src/runner.c b/src/runner.c
index 04948926d584593b791e0ccf5f21cac8847338a1..6c53bf0881682f2b8551ca023c18f5f0e109b768 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1350,7 +1350,13 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the part. */
       struct part *restrict p = &parts[k];
-      if (part_is_active(p, e)) hydro_end_force(p);
+
+      if (part_is_active(p, e)) {
+
+        /* First, finish the force loop */
+        hydro_end_force(p);
+        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
+      }
     }
 
     /* Loop over the g-particles in this cell. */
@@ -1358,7 +1364,11 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the gpart. */
       struct gpart *restrict gp = &gparts[k];
-      if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
+
+      if (gp->type == swift_type_dark_matter) {
+
+        if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
+      }
 
 #ifdef SWIFT_DEBUG_CHECKS
       if (e->policy & engine_policy_self_gravity) {
@@ -1377,7 +1387,12 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the spart. */
       struct spart *restrict sp = &sparts[k];
-      if (spart_is_active(sp, e)) star_end_force(sp);
+      if (spart_is_active(sp, e)) {
+
+        /* First, finish the force loop */
+        star_end_force(sp);
+        gravity_end_force(sp->gpart, const_G);
+      }
     }
   }