Fix typos. Add include guards. Add description for advection fields.

CMakeLists.txt

@@ -392,7 +392,7 @@ add_library(crpropa SHARED
 	src/magneticField/MagneticField.cpp
 	src/magneticField/MagneticFieldGrid.cpp
 	src/magneticField/PshirkovField.cpp
-	src/magneticField/ArchmedeanSpiralField.cpp
+	src/magneticField/ArchimedeanSpiralField.cpp
 	src/advectionField/AdvectionField.cpp
 	${CRPROPA_EXTRA_SOURCES}
 )

include/CRPropa.h

@@ -59,7 +59,7 @@
 #include "crpropa/magneticField/MagneticFieldGrid.h"
 #include "crpropa/magneticField/PshirkovField.h"
 #include "crpropa/magneticField/QuimbyMagneticField.h"
-#include "crpropa/magneticField/ArchmedeanSpiralField.h"
+#include "crpropa/magneticField/ArchimedeanSpiralField.h"
 
 #include "crpropa/advectionField/AdvectionField.h"
 

include/crpropa/Units.h

@@ -89,16 +89,15 @@ static const double cm = centimeter;
 
 // second
 static const double nanosecond = 1e-9 * second;
-static const double mikrosecond = 1e-6 * 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 ns = nanosecond;
-static const double mus = mikrosecond;
+static const double mus = microsecond;
 static const double ms = millisecond;
 static const double sec = second;
 static const double min = minute;
-static const double h = hour;
 
 } // namespace crpropa
 

include/crpropa/advectionField/AdvectionField.h

@@ -1,7 +1,7 @@
 #ifndef CRPROPA_ADVECTIONFIELD_H
 #define CRPROPA_ADVECTIONFIELD_H
 
-#pragma once
+
 #include <string>
 #include <iostream>
 #include <cmath>
@@ -44,38 +44,47 @@ class AdvectionFieldList: public AdvectionField {
 
 /**
  @class UniformAdvectionField
- @brief Advection field with one B-field vector.
+ @brief Advection field with one velocity/advection-field vector.
  */
 class UniformAdvectionField: public AdvectionField {
 	Vector3d value;
 public:
 	UniformAdvectionField(const Vector3d &value);
 	Vector3d getField(const Vector3d &position) const;
 	double getDivergence(const Vector3d &position) const;
+
+	std::string getDescription() const;
 };
 
 
 /**
 @class ConstantSphericalAdvectionField
-@brief Spherical advection with a constant wind speed
-
+@brief Spherical advection field with a constant wind speed
 */
 
 class ConstantSphericalAdvectionField: public AdvectionField {
 	Vector3d origin; //origin of the advection sphere
-	double vWind; // maximum wind velocity
+	double vWind; // wind velocity
 public:
-	ConstantSphericalAdvectionField(Vector3d origin, double vWind);
+	/** Constructor
+	 @param origin	Origin of the advection field
+	 @param vWind	Constant wind velocity
+
+*/
+
+	ConstantSphericalAdvectionField(const Vector3d origin, double vWind);
 	Vector3d getField(const Vector3d &position) const;
 	double getDivergence(const Vector3d &position) const;
 
-	void setOrigin(Vector3d origin);
+	void setOrigin(const Vector3d origin);
 	void setVWind(double vMax);
 
 	Vector3d getOrigin() const;
 	double getVWind() const;
 
-	//Add description method
+	std::string getDescription() const;
+
+
 };
 
 /**
@@ -91,13 +100,20 @@ class SphericalAdvectionField: public AdvectionField {
 	double tau; // transition distance
 	double alpha; //tuning parameter
 public:
-	SphericalAdvectionField(Vector3d origin, double radius, double vMax, double tau, double alpha);
+	/** Constructor
+	@param origin 	Origin of the advection sphere
+	@param radius 	Radius of the advection sphere
+	@param vMax	Maximum wind velocity
+	@param tau	Transition distance
+	@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;
 
 	double getV(const double &r) const;
 
-	void setOrigin(Vector3d origin);
+	void setOrigin(const Vector3d origin);
 	void setRadius(double radius);
 	void setVMax(double vMax);
 	void setTau(double tau);
@@ -129,15 +145,23 @@ class SphericalAdvectionShock: public AdvectionField {
 	double v_phi; // rotation speed at r_rot
 
 public:
-	SphericalAdvectionShock(Vector3d origin, double r_0, double v_0, double lambda);
+	/** Constructor
+	@param origin 	Origin of the advection sphere
+	@param r_0 	Position of the shock
+	@param v_0 	Constant velocity (r<<r_o)
+	@param lambda 	Transition width / width of the shock
+	@param r_rot 	Normalization radius for rotation speed
+	@param v_phi 	Rotation speed at r_rot
+*/
+	SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double lambda);
 
 	Vector3d getField(const Vector3d &position) const;
 	double getDivergence(const Vector3d &position) const;
 
 	double g(double R) const;
 	double g_prime(double R) const;
 
-	void setOrigin(Vector3d Origin);
+	void setOrigin(const Vector3d Origin);
 	void setR0(double r);
 	void setV0(double v);
 	void setLambda(double l);
@@ -150,6 +174,8 @@ class SphericalAdvectionShock: public AdvectionField {
 	double getLambda() const;
 	double getRRot() const;
 	double getAzimuthalSpeed() const;
+
+	std::string getDescription() const;
 };
 
 } // namespace crpropa

include/crpropa/magneticField/ArchmedeanSpiralField.h → include/crpropa/magneticField/ArchimedeanSpiralField.h

@@ -1,9 +1,9 @@
-#ifndef CRPROPA_ARCHMEDEANSPIRALFIELD_H
-#define CRPROPA_ARCHMEDEANSPIRALFIELD_H
+#ifndef CRPROPA_ARCHIMEDEANSPIRALFIELD_H
+#define CRPROPA_ARCHIMEDEANSPIRALFIELD_H
 
 #include "crpropa/magneticField/MagneticField.h"
 
-#pragma once
+
 #include <string>
 #include <iostream>
 #include <cmath>
@@ -18,21 +18,27 @@ namespace crpropa {
 
 /**
 
-@class ArchmedeanSpiralField
-@brief Magnetic field model following a archmedean spiral
+@class ArchimedeanSpiralField
+@brief Magnetic field model following a Archimedean spiral
 
 See e.g. Jokipii, Levy & Hubbard 1977
 */
 
-class ArchmedeanSpiralField: public MagneticField {
+class ArchimedeanSpiralField: public MagneticField {
 private:
 	double B_0; // Magnetic field strength at reference level
 	double R_0; // Reference level
 	double Omega; // Angular velocity of the rotation
 	double V_w; // Asymptotic wind speed
 
 public:
-	ArchmedeanSpiralField(double B_0, double R_0, double Omega, double V_w);
+/** Constructor
+	@param B_0	Magnetic field strength at reference level
+	@param R_0	Reference level
+	@param Omega	Angular velocity of the rotation
+	@param V_w	Asymptotic wind speed
+*/
+	ArchimedeanSpiralField(double B_0, double R_0, double Omega, double V_w);
 
 	Vector3d getField(const Vector3d &pos) const;	
 
@@ -50,4 +56,4 @@ class ArchmedeanSpiralField: public MagneticField {
 
 } // end namespace crpropa
 
-#endif // CRPROPA_ACHMEDEANSPIRALFIELD_H
+#endif // CRPROPA_ACHIMEDEANSPIRALFIELD_H

include/crpropa/module/AdiabaticCooling.h

@@ -1,4 +1,6 @@
-#pragma once
+#ifndef CRPROPA_ADIABATICCOOLING_H
+#define CRPROPA_ADIABATICCOOLING_H
+
 #include <string>
 #include <iostream>
 #include <cmath>
@@ -9,6 +11,8 @@
 #include "crpropa/Module.h"
 #include "crpropa/Units.h"
 #include "crpropa/advectionField/AdvectionField.h"
+#include "kiss/logger.h"
+
 
 namespace crpropa {
 
@@ -23,6 +27,10 @@ class AdiabaticCooling: public Module {
 		double limit;
 
 	public:
+	/** Constructor
+	@param advectionField 	The advection field used for the adiabatic energy change
+	@param limit 		Maximum relative energy change allowed	
+*/
 		AdiabaticCooling(ref_ptr<AdvectionField> advectionField);
 		AdiabaticCooling(ref_ptr<AdvectionField> advectionField, double limit);		
 		void process(Candidate *c) const;
@@ -36,3 +44,4 @@ class AdiabaticCooling: public Module {
 
 
 }; // end namesspace crpropa
+#endif // CRPROPA_ADIABATICCOOLING_H

include/crpropa/module/DiffusionSDE.h

@@ -1,13 +1,20 @@
-#pragma once
+#ifndef CRPROPA_DIFFUSIONSDE_H
+#define CRPROPA_DIFFUSIONSDE_H
+
 #include <iostream>
 #include <vector>
 #include <cmath>
 #include <string>
+#include <cstdlib>
+#include <stdexcept>
 
 #include <crpropa/Module.h>
 #include <crpropa/magneticField/MagneticField.h>
 #include <crpropa/advectionField/AdvectionField.h>
 #include <crpropa/Units.h>
+#include <crpropa/Random.h>
+
+#include "kiss/logger.h"
 
 namespace crpropa {
 
@@ -35,6 +42,14 @@ class DiffusionSDE : public Module{
 
 
 public:
+/** Constructor
+	@param minStep		minStep/c_light is the minimum integration timestep
+	@param maxStep		maxStep/c_light is the maximum integration timestep
+	@param tolerance	Tolerance is criterion for step adjustment. Step adjustment takes place when the tangential vector of the magnetic field line is calculated.
+	@param epsilon		Ratio of parallel and perpendicular diffusion coefficient D_par = epsilon*D_perp
+	@param alpha 		Power law index of the energy dependent diffusion coefficient: D\propto E^alpha
+	@param scale 		Scaling factor for the diffusion coefficient D = scale*D_0
+*/
 	    DiffusionSDE(ref_ptr<crpropa::MagneticField> magneticField, double tolerance = 1e-4, double minStep=(10*pc), double maxStep=(1*kpc), double epsilon=0.1);
 
 	    DiffusionSDE(ref_ptr<crpropa::MagneticField> magneticField, ref_ptr<crpropa::AdvectionField> advectionField, double tolerance = 1e-4, double minStep=(10*pc), double maxStep=(1*kpc), double epsilon=0.1);
@@ -65,3 +80,5 @@ class DiffusionSDE : public Module{
 }; 
 
 } //namespace crpropa
+
+#endif // CRPROPA_DIFFUSIONSDE_H

python/2_headers.i

@@ -308,7 +308,7 @@
 %include "crpropa/magneticField/AMRMagneticField.h"
 %include "crpropa/magneticField/JF12Field.h"
 %include "crpropa/magneticField/PshirkovField.h"
-%include "crpropa/magneticField/ArchmedeanSpiralField.h"
+%include "crpropa/magneticField/ArchimedeanSpiralField.h"
 %include "crpropa/module/BreakCondition.h"
 %include "crpropa/module/Boundary.h"
 

src/advectionField/AdvectionField.cpp

@@ -38,9 +38,15 @@ double UniformAdvectionField::getDivergence(const Vector3d &position) const {
 	return 0.;
 	}
 
+std::string UniformAdvectionField::getDescription() const {
+	std::stringstream s;
+	s << "v0: " << value / km * sec << " km/s, ";
+	return s.str();
+}
+
 //----------------------------------------------------------------
 
-ConstantSphericalAdvectionField::ConstantSphericalAdvectionField(Vector3d origin, double vWind) {
+ConstantSphericalAdvectionField::ConstantSphericalAdvectionField(const Vector3d origin, double vWind) {
 	setOrigin(origin);
 	setVWind(vWind);
 }
@@ -55,7 +61,7 @@ double ConstantSphericalAdvectionField::getDivergence(const Vector3d &position)
 	return 2*vWind/R;
 }
 
-void ConstantSphericalAdvectionField::setOrigin(Vector3d o) {
+void ConstantSphericalAdvectionField::setOrigin(const Vector3d o) {
 	origin=o;
 	return;
 }
@@ -73,11 +79,18 @@ double ConstantSphericalAdvectionField::getVWind() const {
 	return vWind;
 }
 
+std::string ConstantSphericalAdvectionField::getDescription() const {
+	std::stringstream s;
+	s << "Origin: " << origin / kpc  << " kpc, ";
+	s << "v0: " << vWind / km * sec << " km/s, ";
+	return s.str();
+}
+
 
 
 //----------------------------------------------------------------
 
-SphericalAdvectionField::SphericalAdvectionField(Vector3d origin, double radius, double vMax, double tau, double alpha) {
+SphericalAdvectionField::SphericalAdvectionField(const Vector3d origin, double radius, double vMax, double tau, double alpha) {
 	setOrigin(origin);
 	setRadius(radius);
 	setVMax(vMax);
@@ -109,7 +122,7 @@ double SphericalAdvectionField::getV(const double &r) const {
 	return f;
 }
 
-void SphericalAdvectionField::setOrigin(Vector3d o) {
+void SphericalAdvectionField::setOrigin(const Vector3d o) {
 	origin = o;
 	return;
 }
@@ -167,7 +180,7 @@ std::string SphericalAdvectionField::getDescription() const {
 
 //-----------------------------------------------------------------
 
-SphericalAdvectionShock::SphericalAdvectionShock(Vector3d origin, double r_0, double v_0, double l) {
+SphericalAdvectionShock::SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double l) {
 	setOrigin(origin);
 	setR0(r_0);
 	setV0(v_0);
@@ -210,7 +223,7 @@ double SphericalAdvectionShock::g_prime(double r) const {
 	return 1. / (2*lambda*(1+cosh(-a)));
 }	
 
-void SphericalAdvectionShock::setOrigin(Vector3d o) {
+void SphericalAdvectionShock::setOrigin(const Vector3d o) {
 	origin = o;
 }
 
@@ -258,5 +271,15 @@ double SphericalAdvectionShock::getAzimuthalSpeed() const {
 	return v_phi;
 }
 
+std::string SphericalAdvectionShock::getDescription() const {
+	std::stringstream s;
+	s << "Origin: " << origin / kpc  << " kpc, ";
+	s << "r0 (shock radius): " << r_0 / kpc  << " kpc, ";
+	s << "r_rot (norm. azimuthal velocity): " << r_rot / kpc  << " kpc, ";
+	s << "v0 (maximum radial speed): " << v_0 / km * sec << " km/s, ";
+	s << "v_phi (azimuthal speed @ r_rot): " << v_phi / km * sec << " km/s, ";
+	s << "lambda: " << lambda / pc << " pc";
+	return s.str();
+}
 
 } // namespace crpropa