diff --git a/python/cosmotool/cic.py b/python/cosmotool/cic.py index 8fb4554..06553f3 100644 --- a/python/cosmotool/cic.py +++ b/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)) diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py index 4c279db..f3a51ce 100644 --- a/python_sample/icgen/borgicgen.py +++ b/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 diff --git a/python_sample/icgen/cosmogrowth.py b/python_sample/icgen/cosmogrowth.py index 65279c9..bee7fd3 100644 --- a/python_sample/icgen/cosmogrowth.py +++ b/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 diff --git a/python_sample/icgen/gen_ic_from_borg.py b/python_sample/icgen/gen_ic_from_borg.py index ab45742..5f9d210 100644 --- a/python_sample/icgen/gen_ic_from_borg.py +++ b/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)