diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index 0ac0564116f8a6ceb57b4f41d23eb9907df0440d..c2c68a5b93e7ea2067eecbf7bc66f097f0642f22 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -43,14 +43,14 @@ vol = boxSize[0] * boxSize[1] * boxSize[2]
 glass1 = h5py.File("glass_001.hdf5")
 pos1 = glass1["/PartType0/Coordinates"][:,:]
 pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
 
 #Read in high density glass
 # glass2 = h5py.File("../Glass/glass_50000.hdf5")
 glass2 = h5py.File("glass_002.hdf5")
 pos2 = glass2["/PartType0/Coordinates"][:,:]
 pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
 
 #Generate high density region
 rho1 = 1.
@@ -61,9 +61,10 @@ coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
 # coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
 N1 = size(coord1)/3
 v1 = zeros((N1, 3))
-h1 = ones(N1) * 2.251 * 0.5 * vol / (size(pos1)/3)**(1./3.)
 u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
 m1 = ones(N1) * vol * 0.5 * rho1 / N1
+h1 = append(glass_h1, glass_h1, 0)
+h1 = append(h1, h1, 0)
 
 #Generate low density region
 rho2 = 0.25
@@ -74,9 +75,10 @@ coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
 # coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
 N2 = size(coord2)/3
 v2 = zeros((N2, 3))
-h2 = ones(N2) * 2.251 * 0.5 * vol / (size(pos2)/3)**(1./3.)
 u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
 m2 = ones(N2) * vol * 0.5 * rho2 / N2
+h2 = append(glass_h2, glass_h2, 0)
+h2 = append(h2, h2, 0)
 
 #Merge arrays
 numPart = N1 + N2
@@ -89,8 +91,8 @@ ids = zeros(numPart, dtype='L')
 for i in range(1, numPart+1):
     ids[i-1] = i
 
-#Final operations
-h /= 2
+#Final operation since we come from Gadget-2 cubic spline ICs
+h /= 1.825752
 
 #File
 file = h5py.File(fileName, 'w')
diff --git a/src/debug.c b/src/debug.c
index 4c1434118c98aab7def28d3a53493767d249d774..53a03d66aee2c169a555ed00a2efa2d5b984066a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -60,7 +60,7 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
   /* Look for the particle. */
   for (size_t i = 0; i < N; i++)
     if (parts[i].id == id) {
-      printf("## Particle[%zd]:\n id=%lld", i, parts[i].id);
+      printf("## Particle[%zd]:\n id=%lld ", i, parts[i].id);
       hydro_debug_particle(&parts[i], &xparts[i]);
       found = 1;
       break;
@@ -76,12 +76,12 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
   /* Look for the particle. */
   for (size_t i = 0; i < N; i++)
     if (gparts[i].id == -id) {
-      printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
+      printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
       gravity_debug_particle(&gparts[i]);
       found = 1;
       break;
     } else if (gparts[i].id > 0 && gparts[i].part->id == id) {
-      printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
+      printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
       gravity_debug_particle(&gparts[i]);
       found = 1;
       break;
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 46e156bb99015069f9958aeea05954e2be6db5e0..a4d1f7dd4397ebfc850b582e1ca81fc0d4edb76a 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -23,13 +23,13 @@ __attribute__((always_inline))
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "h=%.3e, "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
+      "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
       "dS/dt=%.3e, c=%.3e\n"
       "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e]  \n "
       "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
+      p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
       p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed,
       p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1],
       p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);