diff --git a/include/crpropa/Units.h b/include/crpropa/Units.h index c8add3e48..c85492b3b 100644 --- a/include/crpropa/Units.h +++ b/include/crpropa/Units.h @@ -121,16 +121,26 @@ static const double centimeter = 0.01 * meter; static const double km = kilometer; static const double cm = centimeter; -// second +// time static const double nanosecond = 1e-9 * second; static const double microsecond = 1e-6 * second; static const double millisecond = 1e-3 * second; static const double minute = 60 * second; static const double hour = 3600 * second; +static const double day = 24 * hour; +static const double year = 365.25 * day; +static const double kiloyear = 1e3 * year; +static const double megayear = 1e6 * year; +static const double gigayear = 1e9 * year; + static const double ns = nanosecond; static const double mus = microsecond; static const double ms = millisecond; static const double sec = second; +static const double yr = year; +static const double kyr = kiloyear; +static const double Myr = megayear; +static const double Gyr = gigayear; // volume static const double ccm = cm*cm*cm; diff --git a/include/crpropa/module/HDF5Output.h b/include/crpropa/module/HDF5Output.h index aab716987..96c963d4f 100644 --- a/include/crpropa/module/HDF5Output.h +++ b/include/crpropa/module/HDF5Output.h @@ -54,6 +54,7 @@ class HDF5Output: public Output { typedef struct OutputRow { double D; + double T; double z; uint64_t SN; int32_t ID; diff --git a/include/crpropa/module/Output.h b/include/crpropa/module/Output.h index 76b28fc6e..7c84b43d6 100644 --- a/include/crpropa/module/Output.h +++ b/include/crpropa/module/Output.h @@ -24,6 +24,7 @@ namespace crpropa { or as columns (in a HDF5 file). The right columns are the names of each column for internal access. . D TrajectoryLengthColumn + . T TimeColumn . SN SerialNumberColumn . ID CurrentIdColumn . E CurrentEnergyColumn @@ -53,7 +54,7 @@ namespace crpropa { */ class Output: public Module { protected: - double lengthScale, energyScale; + double lengthScale, energyScale, timeScale; std::bitset<64> fields; struct Property @@ -72,6 +73,7 @@ class Output: public Module { public: enum OutputColumn { TrajectoryLengthColumn, + TimeColumn, ColumnDensityColumn, RedshiftColumn, CurrentIdColumn, @@ -124,6 +126,10 @@ class Output: public Module { /** Set type of output. @param outputType type of output: Trajectory1D, Trajectory3D, Event1D, Event3D, Everything */ + /** Set time scale + @param scale time scale (scale = 1 corresponds to 1 sec) + */ + void setTimeScale(double scale); void setOutputType(OutputType outputType); /** Determines whether a given column will be displayed in the output. @param field name of the field to be added/removed from output diff --git a/src/module/HDF5Output.cpp b/src/module/HDF5Output.cpp index 5dcc5768f..1a4a95a62 100644 --- a/src/module/HDF5Output.cpp +++ b/src/module/HDF5Output.cpp @@ -106,6 +106,8 @@ void HDF5Output::open(const std::string& filename) { sid = H5Tcreate(H5T_COMPOUND, sizeof(OutputRow)); if (fields.test(TrajectoryLengthColumn)) H5Tinsert(sid, "D", HOFFSET(OutputRow, D), H5T_NATIVE_DOUBLE); + if (fields.test(TimeColumn)) + H5Tinsert(sid, "T", HOFFSET(OutputRow, T), H5T_NATIVE_DOUBLE); if (fields.test(RedshiftColumn)) H5Tinsert(sid, "z", HOFFSET(OutputRow, z), H5T_NATIVE_DOUBLE); if (fields.test(SerialNumberColumn)) @@ -205,6 +207,7 @@ void HDF5Output::open(const std::string& filename) { insertStringAttribute("Version", g_GIT_DESC); insertDoubleAttribute("LengthScale", this->lengthScale); insertDoubleAttribute("EnergyScale", this->energyScale); + insertDoubleAttribute("TimeScale", this->timeScale); // add ranom seeds std::vector< std::vector > seeds = Random::getSeedThreads(); @@ -258,6 +261,7 @@ void HDF5Output::process(Candidate* candidate) const { OutputRow r; r.D = candidate->getTrajectoryLength() / lengthScale; + r.T = candidate->getTime() / timeScale; r.z = candidate->getRedshift(); r.SN = candidate->getSerialNumber(); diff --git a/src/module/Output.cpp b/src/module/Output.cpp index dd3978869..29ff0c4ab 100644 --- a/src/module/Output.cpp +++ b/src/module/Output.cpp @@ -5,11 +5,11 @@ namespace crpropa { -Output::Output() : outputName(OutputTypeName(Everything)), lengthScale(Mpc), energyScale(EeV), oneDimensional(false), count(0) { +Output::Output() : outputName(OutputTypeName(Everything)), lengthScale(Mpc), energyScale(EeV), timeScale(Gyr), oneDimensional(false), count(0) { enableAll(); } -Output::Output(OutputType outputType) : outputName(OutputTypeName(outputType)), lengthScale(Mpc), energyScale(EeV), oneDimensional(false), count(0) { +Output::Output(OutputType outputType) : outputName(OutputTypeName(outputType)), lengthScale(Mpc), energyScale(EeV), timeScale(Gyr), oneDimensional(false), count(0) { setOutputType(outputType); } @@ -88,6 +88,11 @@ void Output::setLengthScale(double scale) { lengthScale = scale; } +void Output::setTimeScale(double scale) { + modify(); + timeScale = scale; +} + void Output::set1D(bool value) { modify(); oneDimensional = value; diff --git a/src/module/TextOutput.cpp b/src/module/TextOutput.cpp index 07bd5ea3c..c6a1eac39 100644 --- a/src/module/TextOutput.cpp +++ b/src/module/TextOutput.cpp @@ -54,6 +54,8 @@ void TextOutput::printHeader() const { *out << "#"; if (fields.test(TrajectoryLengthColumn)) *out << "\tD"; + if (fields.test(TimeColumn)) + *out << "\tT"; if (fields.test(RedshiftColumn)) *out << "\tz"; if (fields.test(SerialNumberColumn)) @@ -106,6 +108,8 @@ void TextOutput::printHeader() const { if (fields.test(TrajectoryLengthColumn)) *out << "# D Trajectory length [" << lengthScale / Mpc << " Mpc]\n"; + if (fields.test(TimeColumn)) + *out << "# T Time [" << timeScale / Gyr << " Gyr] \n"; if (fields.test(RedshiftColumn)) *out << "# z Redshift\n"; if (fields.test(SerialNumberColumn)) @@ -169,6 +173,10 @@ void TextOutput::process(Candidate *c) const { p += std::sprintf(buffer + p, "%8.5E\t", c->getTrajectoryLength() / lengthScale); + if (fields.test(TimeColumn)) + p += std::sprintf(buffer + p, "%8.5E\t", + c->getTime() / timeScale); + if (fields.test(RedshiftColumn)) p += std::sprintf(buffer + p, "%1.5E\t", c->getRedshift()); @@ -221,7 +229,6 @@ void TextOutput::process(Candidate *c) const { p += std::sprintf(buffer + p, "%8.5E\t%8.5E\t%8.5E\t", pos.x, pos.y, pos.z); } - } if (fields.test(SerialNumberColumn)) @@ -289,6 +296,7 @@ void TextOutput::load(const std::string &filename, ParticleCollector *collector) double lengthScale = Mpc; // default Mpc double energyScale = EeV; // default EeV + double timeScale = Gyr; // default Gyr if (!infile.good()) throw std::runtime_error("crpropa::TextOutput: could not open file " + filename); @@ -313,6 +321,8 @@ void TextOutput::load(const std::string &filename, ParticleCollector *collector) stream >> val_d; c->setTrajectoryLength(val_d*lengthScale); // D stream >> val_d; + c->setTime(val_d * timeScale); // T + stream >> val_d; c->setRedshift(val_d); // z stream >> val_i; c->setSerialNumber(val_i); // SN