Fixed cic. Fixed lpt1 evaluation in cosmogrowth

python/cosmotool/cic.py

@@ -3,18 +3,25 @@ import numpy as np
 
 def cicParticles(particles, L, N):
 
-    for d in xrange(3):
-        q = particles[d]*N/L
-        i.append(np.floor(q).astype(int))
-        r.append(q-i[-1])
-    
-    density  = np.bincount(shifted(i, (0,0,0)), weights=   r[0]*    r[1]*   r[2], minlength=N*N*N)
-    density += np.bincount(shifted(i, (0,0,1)), weights=   r[0]*    r[1]*(1-r[2]), minlength=N*N*N)
-    density += np.bincount(shifted(i, (0,1,0)), weights=   r[0]*(1-r[1])*   r[2], minlength=N*N*N)
-    density += np.bincount(shifted(i, (0,1,1)), weights=   r[0]*(1-r[1])*(1-r[2]), minlength=N*N*N)
-    density += np.bincount(shifted(i, (1,0,0)), weights=(1-r[0])*    r[1]*   r[2], minlength=N*N*N)
-    density += np.bincount(shifted(i, (1,0,1)), weights=(1-r[0])*    r[1]*(1-r[2]), minlength=N*N*N)
-    density += np.bincount(shifted(i, (1,1,0)), weights=(1-r[0])*(1-r[1])*   r[2], minlength=N*N*N)
-    density += np.bincount(shifted(i, (1,1,1)), weights=(1-r[0])*(1-r[1])*(1-r[2]), minlength=N*N*N)
 
-    return density
+    def shifted(i, t):
+      return (i[2]+t[2])%N + N*((i[1]+t[1])%N + N*((i[0]+t[0])%N))
+
+    i =[]
+    r = []
+    for d in xrange(3):
+        q = (particles[d]%L)*N/L
+        o = np.floor(q).astype(int)
+        i.append(o)
+        r.append(q-o)
+    
+    density  = np.bincount(shifted(i, (1,1,1)), weights=   r[0] *   r[1] *   r[2], minlength=N*N*N)
+    density += np.bincount(shifted(i, (1,1,0)), weights=   r[0] *   r[1] *(1-r[2]), minlength=N*N*N)
+    density += np.bincount(shifted(i, (1,0,1)), weights=   r[0] *(1-r[1])*   r[2], minlength=N*N*N)
+    density += np.bincount(shifted(i, (1,0,0)), weights=   r[0] *(1-r[1])*(1-r[2]), minlength=N*N*N)
+    density += np.bincount(shifted(i, (0,1,1)), weights=(1-r[0])*   r[1] *   r[2], minlength=N*N*N)
+    density += np.bincount(shifted(i, (0,1,0)), weights=(1-r[0])*   r[1] *(1-r[2]), minlength=N*N*N)
+    density += np.bincount(shifted(i, (0,0,1)), weights=(1-r[0])*(1-r[1])*   r[2], minlength=N*N*N)
+    density += np.bincount(shifted(i, (0,0,0)), weights=(1-r[0])*(1-r[1])*(1-r[2]), minlength=N*N*N)
+
+    return density.reshape((N,N,N))

python_sample/icgen/borgicgen.py

@@ -34,6 +34,7 @@ def compute_power_from_borg(input_borg, a_borg, cosmo, bins=10, range=(0,1)):
   cgrowth = CosmoGrowth(**cosmo)
   D1 = cgrowth.D(1)
   D1_0 = D1/cgrowth.D(a_borg)
+  print("D1_0=%lg" % D1_0)
 
   density_hat, L = ba.half_pixel_shift(borg_vol)
 
@@ -76,7 +77,7 @@ def run_generation(input_borg, a_borg, a_ic, **cosmo):
 
     for j in xrange(3):
         # Generate psi_j (displacement along j)
-        psi = D1_0*lpt.lpt1(j).flatten()*(N/L)**3
+        psi = D1_0*lpt.lpt1(j).flatten()
         # Generate posx
         posx.append(((posq[j] + psi)%L).astype(np.float32))
         # Generate vel

python_sample/icgen/cosmogrowth.py

@@ -57,7 +57,7 @@ class LagrangianPerturbation(object):
         k2 = self._get_k2()
         k2[0,0,0] = 1
     
-        return self._gradient(-self.dhat/k2, direction)
+        return self._gradient(self.dhat/k2, direction)
         
     def new_shape(self,direction, q=3, half=False):
         N0 = (self.N/2+1) if half else self.N

python_sample/icgen/gen_ic_from_borg.py

@@ -10,5 +10,9 @@ cosmo['SIGMA8']=0.8344
 zstart=0
 astart=1/(1.+zstart)
 
+pos,_,density,N,L,_ = bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo)
+
+dcic = ct.cicParticles(pos, L, N) 
+
 #if __name__=="__main__":
 #    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, **cosmo), **cosmo)