diff --git a/python_sample/build_nbody_skymap.py b/python_sample/build_nbody_skymap.py
new file mode 100644
index 0000000..1aa948a
--- /dev/null
+++ b/python_sample/build_nbody_skymap.py
@@ -0,0 +1,19 @@
+import cosmotool as ct
+import h5py as h5
+
+L=600.
+Nside=128
+
+with h5.File("fields.h5", mode="r") as f:
+    density = f["density"][:]
+
+N = density.shape[0]
+ix = (np.arange(N)-0.5)*L/N - 0.5 * L
+
+
+dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
+
+flux = density / dist2
+projsky1 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=1)
+projsky0 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=0)
+
diff --git a/python_sample/test_spheric_proj.py b/python_sample/test_spheric_proj.py
index f3f80cc..06a1398 100644
--- a/python_sample/test_spheric_proj.py
+++ b/python_sample/test_spheric_proj.py
@@ -14,7 +14,7 @@ zz = ii[None,None,:].repeat(256,axis=0).repeat(256,axis=1).reshape(256**3)
 d_high = ct.interp3d(xx, yy, zz, d, 64, periodic=True)
 d_high = d_high.reshape((256,256,256))
 
-#proj0 = ct.spherical_projection(64, d, 0, 20, integrator_id=0, shifter=np.array([0.5,0.5,0.5]))
-#proj1 = ct.spherical_projection(64, d, 0, 20, integrator_id=1)
+proj0 = ct.spherical_projection(64, d, 0, 20, integrator_id=0, shifter=np.array([0.5,0.5,0.5]))
+proj1 = ct.spherical_projection(64, d, 0, 20, integrator_id=1)
 
 proj0_high = ct.spherical_projection(256, d_high, 0, 30, integrator_id=0, shifter=np.array([0.5,0.5,0.5]))