diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py
index 79ccb0c..8a627ea 100644
--- a/python_sample/icgen/borgicgen.py
+++ b/python_sample/icgen/borgicgen.py
@@ -13,7 +13,7 @@ def gen_posgrid(N, L):
     y = ix[None,:,None].repeat(N, axis=0).repeat(N, axis=2)
     z = ix[None,None,:].repeat(N, axis=0).repeat(N, axis=1)
 
-    return x.flatten(), y.flatten(), z.flatten()
+    return x.reshape((x.size,)), y.reshape((y.size,)), z.reshape((z.size,))
 
 def bin_power(P, L, bins=20, range=(0,1.)):
 
@@ -73,15 +73,17 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, do_lpt2=True,
     D1_0 = D1/cgrowth.D(a_borg)
     velmul = cgrowth.compute_velmul(a_ic)
     
-    D2 = -3./7 * D1_0**2
+    D2 = 3./7 * D1_0**2
 
     for j in xrange(3):
         # Generate psi_j (displacement along j)
         print("LPT1 axis=%d" % j)
-        psi = D1_0*lpt.lpt1(j).flatten()
+        psi = D1_0*lpt.lpt1(j)
+        psi = psi.reshape((psi.size,))
         if do_lpt2:
           print("LPT2 axis=%d" % j)
-          psi += D2 * lpt.lpt2(j).flatten()
+          psi2 = lpt.lpt2(j)
+          psi += D2 * psi2.reshape((psi2.size,))
         # Generate posx
         posx.append(((posq[j] + psi)%L).astype(np.float32))
         # Generate vel
@@ -99,7 +101,9 @@ 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)).transpose(), L, a_ic, **cosmo)
+        ct.writeGrafic("ic_velc%s" % c, vel[i].reshape((N,N,N)), L, a_ic, **cosmo)
+# This used to be necessary. However this has been fixed in writeGrafic now
+#        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 9030a99..2fb1d9b 100644
--- a/python_sample/icgen/gen_ic_from_borg.py
+++ b/python_sample/icgen/gen_ic_from_borg.py
@@ -7,6 +7,7 @@ cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
 cosmo['omega_k_0'] = 0
 cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
+cosmo['ns']=0.9624
 
 zstart=50
 astart=1/(1.+zstart)
diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py
index cfedb09..7f069a8 100644
--- a/python_sample/icgen/test_ic_from_borg.py
+++ b/python_sample/icgen/test_ic_from_borg.py
@@ -1,6 +1,7 @@
 import numpy as np
 import cosmotool as ct
 import borgicgen as bic
+import cosmogrowth as cg
 import sys
 
 cosmo={'omega_M_0':0.3175, 'h':0.6711}
@@ -8,12 +9,25 @@ cosmo['omega_lambda_0']=1-cosmo['omega_M_0']
 cosmo['omega_k_0'] = 0
 cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
-N0=128
+cosmo['ns']=0.9624
+N0=256
+
+doSimulation=True
 
 snap_id=int(sys.argv[1])
+astart=1/100.
 
-s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim/", snap_id, doublePrecision=True)
-astart=s.getTime()
+if doSimulation:
+  s = ct.loadRamsesAll("/nethome/lavaux/remote2/borgsim/", snap_id, doublePrecision=True)
+  astart=s.getTime()
+  L = s.getBoxsize()
+
+  dsim = ct.cicParticles(s.getPositions(), L, N0)
+  dsim /= np.average(np.average(np.average(dsim, axis=0), axis=0), axis=0)
+  dsim -= 1
+
+  dsim_hat = np.fft.rfftn(dsim)*(L/N0)**3
+  Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, range=(0,1.), bins=150)
 
 pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=2, do_lpt2=True)
 
@@ -21,16 +35,20 @@ dcic = ct.cicParticles(pos, L, N0)
 dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
 dcic -= 1
 
-dsim = ct.cicParticles(s.getPositions(), L, N0)
-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
+dcic_hat = np.fft.rfftn(dcic)*(L/N0)**3
+dens_hat = np.fft.rfftn(density)*(L/N0)**3
 
 Pcic, bcic = bic.bin_power(np.abs(dcic_hat)**2/L**3, L, range=(0,1.), bins=150)
-Psim, bsim = bic.bin_power(np.abs(dsim_hat)**2/L**3, L, range=(0,1.), bins=150)
+Pdens, bdens = bic.bin_power(np.abs(dens_hat)**2/L**3, L, range=(0,1.), bins=150)
+
+cgrowth = cg.CosmoGrowth(**cosmo)
+D1 = cgrowth.D(astart)
+D1_0 = D1/cgrowth.D(1)#0.001)
+
+Pref, bref = bic.compute_ref_power(L, N0, cosmo, range=(0,1.), bins=150)
+
+Pcic /= D1_0**2
+Pdens /= D1_0**2
 
 borg_evolved = ct.read_borg_vol("final_density_1380.dat")