diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index 177f950ea248a8a25701e68c45e6dfe94ecb469a..601210deac10350bee475b6935745cef7b7e4f23 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.0    # The end time of the simulation (in internal units).
+  time_end:   10.0    # 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).
 
diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py
index b7f1b213c4e0ea8e57af917988341400afae91b3..854bdf223cfae75203a1924b4af6136b4b7aa6cd 100644
--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
@@ -21,7 +21,7 @@ def do_binning(x,y,x_bin_edges):
     return(count,y_totals)
 
 
-n_snaps = 11
+n_snaps = 100
 
 #for the plotting
 n_radial_bins = int(sys.argv[1])
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index 09e7734b0eb29e20220c8444486028b573e2b9cd..376423fa65131761594cfa25c3beb3252c3750fb 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -90,17 +90,6 @@ grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velo
 grp.attrs["Unit current in cgs (U_I)"] = 1.
 grp.attrs["Unit temperature in cgs (U_T)"] = 1.
 
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [N, 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, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 3
 
 # Runtime parameters
 grp = file.create_group("/RuntimePars")
@@ -109,12 +98,10 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
 # set seed for random number
 np.random.seed(1234)
 
-# Particle group
-grp = file.create_group("/PartType0")
 
 # Positions
 # r^(-2) distribution corresponds to uniform distribution in radius
-radius = np.random.rand(N)
+radius = boxSize  * np.sqrt(3.) / 2.* np.random.rand(N) #the diagonal extent of the cube
 ctheta = -1. + 2 * np.random.rand(N)
 stheta = np.sqrt(1.-ctheta**2)
 phi    =  2 * math.pi * np.random.rand(N)
@@ -133,9 +120,84 @@ print np.mean(coords[:,0])
 print np.mean(coords[:,1])
 print np.mean(coords[:,2])
 
+#now find the particles which are within the box
+
+x_coords = coords[:,0]
+y_coords = coords[:,1]
+z_coords = coords[:,2]
+
+ind = np.where(x_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(x_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(y_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(y_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(z_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(z_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+#count number of particles
+
+N = x_coords.size
+
+print "Number of particles in the box = " , N
+
+#make the coords and radius arrays again
+coords_2 = np.zeros((N,3))
+coords_2[:,0] = x_coords
+coords_2[:,1] = y_coords
+coords_2[:,2] = z_coords
+
+radius = np.sqrt(coords_2[:,0]**2 + coords_2[:,1]**2 + coords_2[:,2]**2)
+
+#test we've done it right
+
+print "x range = (%f,%f)" %(np.min(coords_2[:,0]),np.max(coords_2[:,0]))
+print "y range = (%f,%f)" %(np.min(coords_2[:,1]),np.max(coords_2[:,1]))
+print "z range = (%f,%f)" %(np.min(coords_2[:,2]),np.max(coords_2[:,2]))
+
+print np.mean(coords_2[:,0])
+print np.mean(coords_2[:,1])
+print np.mean(coords_2[:,2])
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [N, 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, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+# Particle group
+grp = file.create_group("/PartType0")
+
 ds = grp.create_dataset('Coordinates', (N, 3), 'd')
-ds[()] = coords
-coords = np.zeros(1)
+ds[()] = coords_2
+coords_2 = np.zeros(1)
 
 # All velocities set to zero
 v = np.zeros((N,3))
diff --git a/examples/HydrostaticHalo/radial_profile.py b/examples/HydrostaticHalo/radial_profile.py
index 781158ba069f49e1b08fa8f53e6d03c5bdf4167b..335f7089b6835b65cf37e1bcd312a17966c295a7 100644
--- a/examples/HydrostaticHalo/radial_profile.py
+++ b/examples/HydrostaticHalo/radial_profile.py
@@ -6,7 +6,7 @@ import sys
 n_snaps = 11
 
 #for the plotting
-n_radial_bins = int(sys.argv[1])
+#n_radial_bins = int(sys.argv[1])
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
@@ -56,36 +56,45 @@ for i in range(n_snaps):
 
     r = radius_over_virial_radius
 
-    bin_width = 1./n_radial_bins
-    hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
+    # bin_width = 1./n_radial_bins
+#     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
 
-#find the mass in each radial bin
+# #find the mass in each radial bin
 
-    mass_dset = f["PartType0/Masses"]
-#mass of each particles should be equal
-    part_mass = np.array(mass_dset)[0]
-    part_mass_cgs = part_mass * unit_mass_cgs
-    part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+#     mass_dset = f["PartType0/Masses"]
+# #mass of each particles should be equal
+#     part_mass = np.array(mass_dset)[0]
+#     part_mass_cgs = part_mass * unit_mass_cgs
+#     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
 
-    mass_hist = hist * part_mass_over_virial_mass
-    radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
-#volume in each radial bin
-    volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+#     mass_hist = hist * part_mass_over_virial_mass
+#     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
+# #volume in each radial bin
+#     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
 
-#now divide hist by the volume so we have a density in each bin
+# #now divide hist by the volume so we have a density in each bin
 
-    density = mass_hist / volume
+#     density = mass_hist / volume
 
-    t = np.linspace(0.01,1.0,1000)
+    # read the densities
+
+    density_dset = f["PartType0/Density"]
+    density = np.array(density_dset)
+    density_cgs = density * unit_mass_cgs / unit_length_cgs**3
+    rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+
+    t = np.linspace(0.01,2.0,1000)
     rho_analytic = t**(-2)/(4.*np.pi)
 
-    plt.plot(radial_bin_mids,density,'ko',label = "Numerical solution")
+    plt.plot(r,rho,'x',label = "Numerical solution")
     plt.plot(t,rho_analytic,label = "Analytic Solution")
     plt.legend(loc = "upper right")
     plt.xlabel(r"$r / r_{vir}$")
     plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0.1,40))
+    #plt.ylim((0.1,40))
+    plt.xscale('log')
+    plt.yscale('log')
     plot_filename = "density_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")
     plt.close()
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index aa3f2ae3074c88c0b462b1b0185e355aae319cc5..c13a7a8ea7e152e0ea6b589aa68af24f8b61ecaa 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -2,7 +2,7 @@
 
 # Generate the initial conditions if they are not present.
 echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 10000 
+python makeIC_full_box.py 10000 
 
 ../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
 
diff --git a/examples/HydrostaticHalo/test_energy_conservation.py b/examples/HydrostaticHalo/test_energy_conservation.py
index 3404f0b5e55fbd5dd249c20fa209d866f0c17a4e..7d06b8d46bb6c58ad9072d17848410e178301375 100644
--- a/examples/HydrostaticHalo/test_energy_conservation.py
+++ b/examples/HydrostaticHalo/test_energy_conservation.py
@@ -3,7 +3,7 @@ import h5py as h5
 import matplotlib.pyplot as plt
 import sys
 
-n_snaps = 101
+n_snaps = 1000
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0
@@ -63,13 +63,13 @@ for i in range(n_snaps):
 
     vels_dset = f["PartType0/Velocities"]
     vels = np.array(vels_dset)
-    speed_squared = coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2
+    speed_squared = vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2
     total_kinetic_energy = 0.5 * np.sum(speed_squared)
     kinetic_energy_array = np.append(kinetic_energy_array,total_kinetic_energy)
 
     u_dset = f["PartType0/InternalEnergy"]
     u = np.array(u_dset)
-    total_internal_energy = 0.5 * np.sum(u)
+    total_internal_energy = np.sum(u)
     internal_energy_array = np.append(internal_energy_array,total_internal_energy)
 
 #put energies in units of v_c^2 and rescale by number of particles
@@ -79,11 +79,6 @@ ke = kinetic_energy_array / (N*v_c**2)
 ie = internal_energy_array / (N*v_c**2)
 te = pe + ke + ie
 
-print pe
-print ke
-print ie
-print te
-
 dyn_time_cgs = r_vir_cgs / v_c_cgs
 time_array = time_array_cgs / dyn_time_cgs
 
@@ -98,3 +93,4 @@ plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %
 plt.ylim((-2,2))
 #plot_filename = "density_profile_%03d.png" %i
 plt.show()
+
diff --git a/src/const.h b/src/const.h
index 19f04f0098854836ccb5873e5065b5c89bd10e2d..379bbe4b630ea401f296bdb267e24b33ba2d8f24 100644
--- a/src/const.h
+++ b/src/const.h
@@ -94,7 +94,8 @@
 
 //#define EXTERNAL_POTENTIAL_NONE
 //#define EXTERNAL_POTENTIAL_POINTMASS
-#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+#define EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL
 //#define EXTERNAL_POTENTIAL_DISC_PATCH
 
 /* Source terms */
diff --git a/src/potential.h b/src/potential.h
index de23882770c4473aabbc41da1c4f18e5d13a0c03..9b16a0e7aa09a3336d8fa8fe03b77a59ba6cee84 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -37,6 +37,8 @@
 #include "./potential/point_mass/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL)
 #include "./potential/isothermal/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL)
+#include "./potential/softened_isothermal/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
 #include "./potential/disc_patch/potential.h"
 #else
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index c6c504c5addd52de3065488f842b7039ec365f18..0919d03fbb2d667456c44763140579314c109615 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
 
   g->a_grav[0] = term * dx;
   g->a_grav[1] = term * dy;
-  g->a_grav[2] = term * dz;
+  g->a_grav[2] = term * dz; 
 }
 
 /**
diff --git a/src/task.c b/src/task.c
index 3c7e638e0e6c4072e4acec1ac04105585c99f052..e58f80d76b5bdf99b37ddcfc7b480330a3d66e93 100644
--- a/src/task.c
+++ b/src/task.c
@@ -39,7 +39,6 @@
 
 /* This object's header. */
 #include "task.h"
-
 /* Local headers. */
 #include "atomic.h"
 #include "error.h"