Merge pull request #525 from LeanderSchlegel/time_master

Explicit time in CRPropa

Author: rafaelab

CMakeLists.txt

@@ -405,6 +405,7 @@ add_library(crpropa SHARED
   src/magneticField/KST24Field.cpp
   src/magneticField/CMZField.cpp
   src/advectionField/AdvectionField.cpp
+  src/advectionField/TimeDependentAdvectionField.cpp
   src/massDistribution/ConstantDensity.cpp
   src/massDistribution/Cordes.cpp
   src/massDistribution/Ferriere.cpp

doc/pages/acceleration.rst

@@ -30,7 +30,12 @@ The exact nature of the dependency is subject of ongoing research. Here,
 so far only quasi-linear theory, a simplistic approach to diffusion of charged
 particles in magnetic fields is implemented in the corresponding modifier :cpp:class:`crpropa::QuasiLinearTheory`.
 
+In the diffusive picture, the diffusive picture as an interplay between diffusion, 
+advection and adiabatic cooling. Simple shock configurations are shown in the following 
+example notebook using the DiffusionSDE module:
 
+.. toctree::
+   example_notebooks/acceleration/diffusive_shock_acceleration.ipynb
 
 
 

doc/pages/example_notebooks/acceleration/DSA-moving-shock.ipynb

@@ -40,7 +40,7 @@
     "import pandas as pd\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
-    "marker= ['.', 's', '^', 'v', '<', '>', 'd', '8']"
+    "marker = ['.', 's', '^', 'v', '<', '>', 'd', '8']"
    ]
   },
   {

doc/pages/example_notebooks/acceleration/diffusive_shock_acceleration.ipynb

@@ -47,6 +47,12 @@ example notebook.
 .. toctree::
    example_notebooks/acceleration/diffusive_shock_acceleration.ipynb
 
+In the next example time-dependent advection fields are introduced and acceleration
+at a moving shock and stationary shock are compared.
+
+.. toctree::
+    example_notebooks/acceleration/DSA-moving-shock.ipynb
+
 Momentum Diffusion 
 ^^^^^^^^^^^^^^^^^^
 

doc/pages/galactic_cosmic_rays.rst

@@ -74,6 +74,7 @@
 #include "crpropa/magneticField/turbulentField/TurbulentField.h"
 
 #include "crpropa/advectionField/AdvectionField.h"
+#include "crpropa/advectionField/TimeDependentAdvectionField.h"
 
 #include "crpropa/massDistribution/Density.h"
 #include "crpropa/massDistribution/Nakanishi.h"

include/CRPropa.h

@@ -47,6 +47,7 @@ class Candidate: public Referenced {
 	double currentStep; /**< Size of the currently performed step in [m] comoving units */
 	double nextStep; /**< Proposed size of the next propagation step in [m] comoving units */
 	std::string tagOrigin; /**< Name of interaction/source process which created this candidate*/
+	double time; /**< Time [s] that has passed in the laboratory frame of reference */
 
 	static uint64_t nextSerialNumber;
 	uint64_t serialNumber;
@@ -73,6 +74,8 @@ class Candidate: public Referenced {
 
 	void setTrajectoryLength(double length);
 	double getTrajectoryLength() const;
+	
+	double getVelocity() const;
 
 	void setRedshift(double z);
 	double getRedshift() const;
@@ -105,6 +108,12 @@ class Candidate: public Referenced {
 	void setTagOrigin(std::string tagOrigin);
 	std::string getTagOrigin() const;
 
+	/**
+	 Sets the time of the candidate.
+	 */
+	void setTime(double t);
+	double getTime() const;
+
 	/**
 	 Make a bid for the next step size: the lowest wins.
 	 */

include/crpropa/Candidate.h

@@ -127,10 +127,19 @@ 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 * 24 * hour;
+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/Units.h

@@ -24,8 +24,8 @@ class AdvectionField: public Referenced {
 public:
 	virtual ~AdvectionField() {
 	}
-	virtual Vector3d getField(const Vector3d &position) const = 0;
-	virtual double getDivergence(const Vector3d &position) const = 0;
+	virtual Vector3d getField(const Vector3d &position, const double &time=0) const = 0;
+	virtual double getDivergence(const Vector3d &position, const double &time=0) const = 0;
 };
 
 
@@ -37,8 +37,8 @@ class AdvectionFieldList: public AdvectionField {
 	std::vector<ref_ptr<AdvectionField> > fields;
 public:
 	void addField(ref_ptr<AdvectionField> field);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 };
 
 
@@ -50,8 +50,8 @@ class UniformAdvectionField: public AdvectionField {
 	Vector3d value;
 public:
 	UniformAdvectionField(const Vector3d &value);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	std::string getDescription() const;
 };
@@ -73,8 +73,8 @@ class ConstantSphericalAdvectionField: public AdvectionField {
 */
 
 	ConstantSphericalAdvectionField(const Vector3d origin, double vWind);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	void setOrigin(const Vector3d origin);
 	void setVWind(double vMax);
@@ -108,8 +108,8 @@ class SphericalAdvectionField: public AdvectionField {
 	@param alpha	Tuning parameter
 */
 	SphericalAdvectionField(const Vector3d origin, double radius, double vMax, double tau, double alpha);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	double getV(const double &r) const;
 
@@ -144,8 +144,8 @@ class OneDimensionalCartesianShock: public AdvectionField {
 	@param lShock //shock width
 */
 	OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	void setComp(double compressionRatio);
 	void setVup(double vUp);
@@ -178,8 +178,8 @@ class OneDimensionalSphericalShock: public AdvectionField {
 	@param coolUpstream //flag for upstream cooling
 */
 	OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	void setComp(double compressionRatio);
 	void setVup(double vUp);
@@ -217,8 +217,8 @@ class ObliqueAdvectionShock: public AdvectionField {
 	
 */
 	ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock);
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	void setComp(double compressionRatio);
 	void setVup(double vXUp);
@@ -257,8 +257,8 @@ class SphericalAdvectionShock: public AdvectionField {
 */
 	SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double lambda);
 
-	Vector3d getField(const Vector3d &position) const;
-	double getDivergence(const Vector3d &position) const;
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
 
 	double g(double R) const;
 	double g_prime(double R) const;

include/crpropa/advectionField/AdvectionField.h

@@ -0,0 +1,88 @@
+#ifndef CRPROPA_TIMEDEPENDENTADVECTIONFIELD_H
+#define CRPROPA_TIMEDEPENDENTADVECTIONFIELD_H
+
+#include "crpropa/advectionField/AdvectionField.h"
+
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+#include <sstream>
+
+#include "crpropa/Vector3.h"
+#include "crpropa/Referenced.h"
+#include "crpropa/Units.h"
+
+namespace crpropa {
+
+/**
+ @class OneDimensionalTimeDependentShock
+ @brief Advection field in x-direction with shock at x = x0 that propagates with a constant speed vsh = v_up - v_down
+ 		and width x_sh approximated by a tanh() with variable compression ratio r_comp = v_up/v_down.
+        Pre- and postshock speeds as well as shock speed must be specified. 
+        The shock position at t=0 can be specified, as well as the time t0 the shock starts to propagate.
+ */
+class OneDimensionalTimeDependentShock: public AdvectionField {
+	double v_sh; // shock speed assuming xsh = vsh * t + xsh_0
+	double v1; // speed behind the shock in lab frame
+	double v0; // undisturbed speed in lab frame
+	double l_sh; // shock width
+    double x_sh0; // shock position at t = 0
+	double t_sh0; // time the shock starts to propagate, before: vsh=0
+public:/** Constructor
+	@param v_sh; // shock speed
+    @param v1; // speed behind the shock in lab frame
+    @param v0; // undisturbed speed in lab frame
+	@param l_sh; // shock width
+*/
+	OneDimensionalTimeDependentShock(double v_sh, double v1, double v0, double l_sh);
+
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
+
+	void setShockSpeed(double v_sh);
+	void setSpeeds(double v1, double v0);
+	void setShockWidth(double l_sh);
+	void setShockPosition(double x_sh0);
+	void setShockTime(double t_sh0);
+
+    double getVshock() const;
+    double getV1() const;
+    double getV0() const;
+    double getShockWidth() const;
+    double getShockPosition(double time) const;
+    double getShockTime() const;
+};
+
+/**
+ @class SedovTaylorBlastWave
+ @brief Spherical advection field with shock at R(t), velocity vsh(t) and width l_sh approximated by tanh() 
+		Wind solution is given by Kahn 1975, see also Drury 1983 for acceleration at the expanding shock
+ */
+class SedovTaylorBlastWave: public AdvectionField {
+	double E0; 	    // energy of the explosion
+	double rho0; 	// initial density
+	double l_sh; 	// shock width
+public:/** Constructor
+	@param E0 // energy of the explosion
+	@param rho0	// initial density
+	
+*/
+	SedovTaylorBlastWave(double E0, double rho0, double l_sh);
+	Vector3d getField(const Vector3d &position, const double &time=0) const;
+	double getDivergence(const Vector3d &position, const double &time=0) const;
+
+	void setShockWidth(double l_sh);
+	void setEnergy(double E0);
+	void setDensity(double rho0);
+
+	double getShockRadius(double time) const;
+    double getShockSpeed(double time) const;
+    double getShockWidth() const;
+    double getEnergy() const;
+    double getDensity() const;
+};
+
+} // namespace crpropa
+
+#endif // CRPROPA_TIMEDEPENDENTADVECTIONFIELD_H