diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py
index d76ed0c..df6d151 100644
--- a/python_sample/icgen/borgicgen.py
+++ b/python_sample/icgen/borgicgen.py
@@ -71,7 +71,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True)
     # Compute LPT scaling coefficient
     D1 = cgrowth.D(a_ic)
     D1_0 = D1/cgrowth.D(a_borg)
-    velmul = cgrowth.compute_velmul(a_ic)*D1_0
+    velmul = cgrowth.compute_velmul(a_ic)
     
     D2 = -3./7 * D1_0**2
 
@@ -87,6 +87,8 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True)
         # Generate vel
         vel.append((psi*velmul).astype(np.float32)) 
 
+    print("velmul=%lg" % (cosmo['h']*velmul))
+
     density = np.fft.irfftn(density_hat*D1_0)*(N/L)**3
 
     return posx,vel,density,N,L,a_ic
@@ -97,7 +99,7 @@ def write_icfiles(*generated_ic, **cosmo):
 
     ct.simpleWriteGadget("borg.gad", posx, velocities=vel, boxsize=L, Hubble=cosmo['h'], Omega_M=cosmo['omega_M_0'], time=a_ic) 
     for i,c in enumerate(['x','y','z']):
-        ct.writeGrafic("ic_velc%s" % c, vel[i].reshape((N,N,N)), L, a_ic, **cosmo)
+        ct.writeGrafic("ic_velc%s" % c, vel[i].reshape((N,N,N)).transpose(), L, a_ic, **cosmo)
 
     ct.writeGrafic("ic_deltab", density, L, a_ic, **cosmo)
     ct.writeWhitePhase("white.dat", density)
diff --git a/python_sample/icgen/gen_ic_from_borg.py b/python_sample/icgen/gen_ic_from_borg.py
index 81d5d35..de73e4f 100644
--- a/python_sample/icgen/gen_ic_from_borg.py
+++ b/python_sample/icgen/gen_ic_from_borg.py
@@ -8,23 +8,9 @@ cosmo['omega_k_0'] = 0
 cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
 
-TestCase=True
-zstart=10
+zstart=200
 astart=1/(1.+zstart)
-
-if TestCase:
-  pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=1, do_lpt2=True)
-
-  dcic = ct.cicParticles(pos, L, N)
-  dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
-  dcic -= 1
-
-  dcic_hat = np.fft.rfftn(dcic)*(L/N)**3
-
-  Pcic, bcic = bic.bin_power(np.abs(dcic_hat)**2/L**3, L, bins=50)
-
-  borg_evolved = ct.read_borg_vol("final_density_1380.dat")
+snapid=1
 
 if __name__=="__main__":
-  if not TestCase:
-    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, do_lpt2=True), **cosmo)
+    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, do_lpt2=False), **cosmo)
diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py
new file mode 100644
index 0000000..dcef2a6
--- /dev/null
+++ b/python_sample/icgen/test_ic_from_borg.py
@@ -0,0 +1,34 @@
+import numpy as np
+import cosmotool as ct
+import borgicgen as bic
+
+cosmo={'omega_M_0':0.3175, 'h':0.6711}
+cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
+cosmo['omega_k_0'] = 0
+cosmo['omega_B_0']=0.049
+cosmo['SIGMA8']=0.8344
+
+snap_id=11
+
+s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim/", snap_id, doublePrecision=True)
+astart=s.getTime()
+
+pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=1, do_lpt2=True)
+
+dcic = ct.cicParticles(pos, L, N)
+dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
+dcic -= 1
+
+dsim = ct.cicParticles(s.getPositions(), L, N)
+dsim /= np.average(np.average(np.average(dsim, axis=0), axis=0), axis=0)
+dsim -= 1
+
+
+dcic_hat = np.fft.rfftn(dcic)*(L/N)**3
+dsim_hat = np.fft.rfftn(dsim)*(L/N)**3
+
+Pcic, bcic = bic.bin_power(np.abs(dcic_hat)**2/L**3, L, bins=50)
+Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, bins=50)
+
+borg_evolved = ct.read_borg_vol("final_density_1380.dat")
+