add time to output

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;

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;

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

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<uint32_t> > 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();

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;

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