This commit is contained in:
Guilhem Lavaux 2015-02-26 10:38:39 +01:00
commit ae4187809e
4 changed files with 35 additions and 13 deletions

View File

@ -1,5 +1,6 @@
from libcpp cimport bool
from libcpp cimport string as cppstring
from libcpp cimport vector as cppvector
import numpy as np
cimport numpy as np
from cpython cimport PyObject, Py_INCREF
@ -29,7 +30,7 @@ cdef extern from "cic.hpp" namespace "CosmoTool":
@cython.cdivision(True)
@cython.wraparound(False)
def leanCic(float[:,:] particles, float L, int Resolution):
cdef CICParticles p
cdef cppvector.vector[CICParticles] *p
cdef CICFilter *cic
cdef np.uint64_t i
cdef CICType *field
@ -40,25 +41,34 @@ def leanCic(float[:,:] particles, float L, int Resolution):
cdef np.uint64_t j
cic = new CICFilter(Resolution, L)
print("Reset mesh")
cic.resetMesh()
if particles.shape[1] != 3:
raise ValueError("Particles must be Nx3 array")
p.mass = 1
for i in xrange(particles.shape[0]):
p.coords[0] = particles[i,0]
p.coords[1] = particles[i,1]
p.coords[2] = particles[i,2]
cic.putParticles(&p, 1)
print("Inserting particles")
# p = new cppvector.vector[CICParticles](particles.shape[0])
# for i in xrange(particles.shape[0]):
# *p[i].mass = 1
# *p[i].coords[0] = particles[i,0]
# *p[i].coords[1] = particles[i,1]
# *p[i].coords[2] = particles[i,2]
# cic.putParticles(&p[0], particles.shape[0])
del p
print("Done")
field = <CICType*>0
dummyRes = 0
cic.getDensityField(field, dummyRes)
print("Got to allocate a numpy %dx%dx%d" % (dummyRes, dummyRes,dummyRes))
out_field = np.empty((dummyRes, dummyRes, dummyRes), dtype=np.float64)
out_field0 = out_field.reshape(out_field.size)
out_field_buf = out_field
print("Copy")
for j in xrange(out_field_buf.size):
out_field_buf[j] = field[j]

View File

@ -75,6 +75,9 @@ cdef class CosmologyPower:
self.power.updateCosmology()
def setNormalization(self,A):
self.power.setNormalization(A)
def normalize(self,s8,k_min=-1,k_max=-1):
"""normalize(self, sigma8)

View File

@ -186,6 +186,11 @@ void CICFilter::putParticles(CICParticles *particles, uint32_t N)
for (long p = tUsedMin[t]; p < tUsedMax[t]; p++)
densityGrid[p] += threadedDensity[t][p];
}
for (int t = 0; t < threadUsed; t++)
{
delete[] threadedDensity;
}
}
void CICFilter::getDensityField(CICType*& field, uint32_t& res)

View File

@ -90,7 +90,8 @@ static double powC(double q, double alpha_c)
static double T_tilde_0(double q, double alpha_c, double beta_c)
{
double a = log(M_E + 1.8 * beta_c * q);
static const double c_E = 2.718282; //M_E;
double a = log(c_E + 1.8 * beta_c * q);
return a / ( a + powC(q, alpha_c) * q * q);
}
@ -133,7 +134,7 @@ void CosmoPower::updateHuWigglesConsts()
double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674));
double b2_zd = 0.238 * pow(OmegaEff, 0.223);
double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd));
double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OMEGA_B*h*h, b2_zd));
double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d;
double Req = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_eq;
@ -183,7 +184,10 @@ double CosmoPower::powerHuWiggles(double k)
double q = k / (13.41 * k_eq);
double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c);
double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde);
double T_b = (
T_tilde_0(q, 1, 1) / (1 + pow(xx / 5.2, 2)) +
alpha_b / (1 + pow(beta_b / xx, 3)) * exp(-pow(k/k_silk, 1.4))
) * j_0(k * s_tilde);
double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c;
@ -286,7 +290,7 @@ void CosmoPower::normalize(double k_min, double k_max)
ofstream ff("PP_k.txt");
for (int i = 0; i < 100; i++)
{
double k = pow(10.0, 4.0*i/100.-2);
double k = pow(10.0, 8.0*i/100.-4);
ff << k << " " << power(k) << endl;
}
@ -376,5 +380,5 @@ void CosmoPower::setFunction(CosmoFunction f)
void CosmoPower::setNormalization(double A_K)
{
normPower = A_K/power(0.002);
normPower = A_K;///power(0.002);
}