diff --git a/python/_cosmotool.pyx b/python/_cosmotool.pyx index f2f2214..edf411a 100644 --- a/python/_cosmotool.pyx +++ b/python/_cosmotool.pyx @@ -42,25 +42,28 @@ cdef extern from "loadRamses.hpp" namespace "CosmoTool": class PySimulationBase(object): def getPositions(self): - raise NotImplemented("getPositions is not implemented") + raise NotImplementedError("getPositions is not implemented") def getVelocities(self): - raise NotImplemented("getVelocities is not implemented") + raise NotImplementedError("getVelocities is not implemented") def getIdentifiers(self): - raise NotImplemented("getIdentifiers is not implemented") + raise NotImplementedError("getIdentifiers is not implemented") def getOmega_M(self): - raise NotImplemented("getOmega_M is not implemented") + raise NotImplementedError("getOmega_M is not implemented") def getOmega_Lambda(self): - raise NotImplemented("getOmega_Lambda is not implemented") + raise NotImplementedError("getOmega_Lambda is not implemented") def getTime(self): - raise NotImplemented("getTime is not implemented") + raise NotImplementedError("getTime is not implemented") def getHubble(self): - raise NotImplemented("getHubble is not implemented") + raise NotImplementedError("getHubble is not implemented") + + def getBoxsize(self): + raise NotImplementedError("getBoxsize is not implemented") cdef class Simulation: @@ -86,6 +89,10 @@ cdef class Simulation: def __get__(Simulation self): return self.data.Omega_M + property Omega_Lambda: + def __get__(Simulation self): + return self.data.Omega_Lambda + property positions: def __get__(Simulation self): return self.positions @@ -93,6 +100,10 @@ cdef class Simulation: property velocities: def __get__(Simulation self): return self.velocities + + property identifiers: + def __get__(Simulation self): + return self.identifiers property numParticles: def __get__(Simulation self): @@ -126,7 +137,13 @@ class _PySimulationAdaptor(PySimulationBase): return self.simu.time def getHubble(self): - return self.simul.Hubble + return self.simu.Hubble + + def getOmega_M(self): + return self.simu.Omega_M + + def getOmega_Lambda(self): + return self.simu.Omega_Lambda cdef class ArrayWrapper: cdef void* data_ptr @@ -180,7 +197,7 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ): cdef object wrap_float_array(float *p, np.uint64_t s): - return wrap_array(p, s, np.NPY_FLOAT) + return wrap_array(p, s, np.NPY_FLOAT32) cdef object wrap_int64_array(np.int64_t* p, np.uint64_t s): return wrap_array(p, s, np.NPY_INT64) @@ -192,10 +209,18 @@ cdef object wrap_simudata(SimuData *data, int flags): simu.data = data if flags & NEED_POSITION: simu.positions = [wrap_float_array(data.Pos[i], data.NumPart) for i in xrange(3)] + else: + simu.positions = None + if flags & NEED_VELOCITY: simu.velocities = [wrap_float_array(data.Vel[i], data.NumPart) for i in xrange(3)] + else: + simu.velocities = None + if flags & NEED_GADGET_ID: simu.identifiers = wrap_int64_array(data.Id, data.NumPart) + else: + simu.identifiers = None return simu def loadGadget(str filename, int snapshot_id, bool loadPosition = True, bool loadVelocity = False, bool loadId = False): @@ -221,8 +246,10 @@ def loadGadget(str filename, int snapshot_id, bool loadPosition = True, bool loa def writeGadget(str filename, object simulation): cdef SimuData simdata - cdef np.ndarray[np.float_t, ndim=1] pos, vel + cdef np.ndarray[np.float32_t, ndim=1] pos, vel + cdef np.ndarray[np.int64_t, ndim=1] ids cdef np.int64_t NumPart + cdef int j if not isinstance(simulation,PySimulationBase): raise TypeError("Second argument must be of type SimulationBase") @@ -239,8 +266,10 @@ def writeGadget(str filename, object simulation): simdata.Pos[j] = pos.data simdata.Vel[j] = vel.data - - simdata.BoxSize = simulation.getBoxSize() + + ids = simulation.getIdentifiers() + simdata.Id = ids.data + simdata.BoxSize = simulation.getBoxsize() simdata.time = simulation.getTime() simdata.Hubble = simulation.getHubble() simdata.Omega_M = simulation.getOmega_M() diff --git a/python/cosmotool/__init__.py b/python/cosmotool/__init__.py index bf0df59..c4dcfcd 100644 --- a/python/cosmotool/__init__.py +++ b/python/cosmotool/__init__.py @@ -1,2 +1,57 @@ from _cosmotool import * from borg import read_borg_vol + + +class SimulationBare(PySimulationBase): + + def __init__(self): + pass + + def getPositions(self): + return self.positions + + def getVelocities(self): + return self.velocities + + def getIdentifiers(self): + return self.identifiers + + def getTime(self): + return self.time + + def getHubble(self): + return self.Hubble + + def getBoxsize(self): + return self.boxsize + + def getOmega_M(self): + return self.Omega_M + + def getOmega_Lambda(self): + return self.Omega_Lambda + + +def simpleWriteGadget(filename, positions, boxsize=1.0, Hubble=100, Omega_M=0.30, time=1, velocities=None, identifiers=None): + + s = SimulationBare() + + s.positions = positions + + if velocities: + s.velocities = velocities + else: + s.velocities = [np.zeros(positions[0].size,dtype=np.float32)]*3 + + if identifiers: + s.identifiers = identifiers + else: + s.identifiers = np.arange(positions[0].size, dtype=np.int64) + + s.Hubble = Hubble + s.time = time + s.Omega_M = Omega_M + s.Omega_Lambda = 1-Omega_M + s.boxsize = boxsize + + writeGadget(filename, s) diff --git a/python/cosmotool/borg.py b/python/cosmotool/borg.py index 5ac0f7b..df10fae 100644 --- a/python/cosmotool/borg.py +++ b/python/cosmotool/borg.py @@ -11,16 +11,8 @@ import glob class BorgVolume(object): def __init__(self, density, ranges): - self._density = density - self._ranges = ranges - - property density: - def __get__(self): - return self._density - - property ranges: - def __get__(self): - return self._ranges + self.density = density + self.ranges = ranges def build_filelist(fdir): #builds list of all borg density fields which may be distributed over several directories diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index be1bf14..e121457 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -274,12 +274,16 @@ void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int Gadget int npart[6], npartTotal[6]; float mass[6]; - if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0) + if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0) { + cerr << "Invalid simulation data array" << endl; return; + } f = new UnformattedWrite(fname); - if (f == 0) + if (f == 0) { + cerr << "Cannot create file" << endl; return; + } for (int i = 0; i < 6; i++) {