diff --git a/CMakeLists.txt b/CMakeLists.txt index 02b728b52..3c7ae1d2c 100644 --- a/CMakeLists.txt +++ b/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} ) diff --git a/include/CRPropa.h b/include/CRPropa.h index 5861114e2..fd42a52cc 100644 --- a/include/CRPropa.h +++ b/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" diff --git a/include/crpropa/Units.h b/include/crpropa/Units.h index 0cdd7a259..51f01b152 100644 --- a/include/crpropa/Units.h +++ b/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 diff --git a/include/crpropa/advectionField/AdvectionField.h b/include/crpropa/advectionField/AdvectionField.h index 30d19398e..88de6bfee 100644 --- a/include/crpropa/advectionField/AdvectionField.h +++ b/include/crpropa/advectionField/AdvectionField.h @@ -1,7 +1,7 @@ #ifndef CRPROPA_ADVECTIONFIELD_H #define CRPROPA_ADVECTIONFIELD_H -#pragma once + #include #include #include @@ -44,7 +44,7 @@ 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; @@ -52,30 +52,39 @@ class UniformAdvectionField: public AdvectionField { 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,7 +145,15 @@ 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< #include #include @@ -18,13 +18,13 @@ 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 @@ -32,7 +32,13 @@ class ArchmedeanSpiralField: public MagneticField { 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 diff --git a/include/crpropa/module/AdiabaticCooling.h b/include/crpropa/module/AdiabaticCooling.h index 996482df2..ee3013181 100644 --- a/include/crpropa/module/AdiabaticCooling.h +++ b/include/crpropa/module/AdiabaticCooling.h @@ -1,4 +1,6 @@ -#pragma once +#ifndef CRPROPA_ADIABATICCOOLING_H +#define CRPROPA_ADIABATICCOOLING_H + #include #include #include @@ -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); AdiabaticCooling(ref_ptr advectionField, double limit); void process(Candidate *c) const; @@ -36,3 +44,4 @@ class AdiabaticCooling: public Module { }; // end namesspace crpropa +#endif // CRPROPA_ADIABATICCOOLING_H diff --git a/include/crpropa/module/DiffusionSDE.h b/include/crpropa/module/DiffusionSDE.h index 59f3794b6..3d63ffa4b 100644 --- a/include/crpropa/module/DiffusionSDE.h +++ b/include/crpropa/module/DiffusionSDE.h @@ -1,13 +1,20 @@ -#pragma once +#ifndef CRPROPA_DIFFUSIONSDE_H +#define CRPROPA_DIFFUSIONSDE_H + #include #include #include #include +#include +#include #include #include #include #include +#include + +#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 magneticField, double tolerance = 1e-4, double minStep=(10*pc), double maxStep=(1*kpc), double epsilon=0.1); DiffusionSDE(ref_ptr magneticField, ref_ptr 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 diff --git a/python/2_headers.i b/python/2_headers.i index 1a69d27dd..db2aed4e7 100644 --- a/python/2_headers.i +++ b/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" diff --git a/src/advectionField/AdvectionField.cpp b/src/advectionField/AdvectionField.cpp index f1ce4188d..45fd60d41 100644 --- a/src/advectionField/AdvectionField.cpp +++ b/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 diff --git a/src/magneticField/ArchmedeanSpiralField.cpp b/src/magneticField/ArchimedeanSpiralField.cpp similarity index 59% rename from src/magneticField/ArchmedeanSpiralField.cpp rename to src/magneticField/ArchimedeanSpiralField.cpp index 095c7dab8..23d3bebb8 100644 --- a/src/magneticField/ArchmedeanSpiralField.cpp +++ b/src/magneticField/ArchimedeanSpiralField.cpp @@ -1,15 +1,15 @@ -#include "crpropa/magneticField/ArchmedeanSpiralField.h" +#include "crpropa/magneticField/ArchimedeanSpiralField.h" namespace crpropa { -ArchmedeanSpiralField::ArchmedeanSpiralField(double B_0, double R_0, double Omega, double V_w) { +ArchimedeanSpiralField::ArchimedeanSpiralField(double B_0, double R_0, double Omega, double V_w) { setB0(B_0); setR0(R_0); setOmega(Omega); setVw(V_w); } -Vector3d ArchmedeanSpiralField::getField(const Vector3d &pos) const { +Vector3d ArchimedeanSpiralField::getField(const Vector3d &pos) const { double r = pos.getR(); double theta = pos.getTheta(); @@ -45,40 +45,40 @@ Vector3d ArchmedeanSpiralField::getField(const Vector3d &pos) const { return B; } -void ArchmedeanSpiralField::setR0(double R) { +void ArchimedeanSpiralField::setR0(double R) { R_0 = R; return; } -void ArchmedeanSpiralField::setB0(double B) { +void ArchimedeanSpiralField::setB0(double B) { B_0 = B; return; } -void ArchmedeanSpiralField::setOmega(double Om) { +void ArchimedeanSpiralField::setOmega(double Om) { Omega = Om; return; } -void ArchmedeanSpiralField::setVw(double v) { +void ArchimedeanSpiralField::setVw(double v) { V_w = v; return; } -double ArchmedeanSpiralField::getR0() const { +double ArchimedeanSpiralField::getR0() const { return R_0; } -double ArchmedeanSpiralField::getB0() const{ +double ArchimedeanSpiralField::getB0() const{ return B_0; } -double ArchmedeanSpiralField::getOmega() const{ +double ArchimedeanSpiralField::getOmega() const{ return Omega; } -double ArchmedeanSpiralField::getVw() const { +double ArchimedeanSpiralField::getVw() const { return V_w; } diff --git a/src/module/AdiabaticCooling.cpp b/src/module/AdiabaticCooling.cpp index ab1754022..66a7db1dc 100644 --- a/src/module/AdiabaticCooling.cpp +++ b/src/module/AdiabaticCooling.cpp @@ -22,8 +22,8 @@ void AdiabaticCooling::process(Candidate *c) const { Div += advectionField->getDivergence(pos); } catch (std::exception &e) { - std::cerr << "AdiabaticCooling: Exception in getDivergence." << std::endl; - std::cerr << e.what() << std::endl; + KISS_LOG_ERROR << "AdiabaticCooling: Exception in getDivergence.\n" + << e.what(); } double dEdt = -E / 3. * Div; // cooling due to advection -p/3 * div(V_wind) diff --git a/src/module/DiffusionSDE.cpp b/src/module/DiffusionSDE.cpp index c7278e460..0b540f396 100644 --- a/src/module/DiffusionSDE.cpp +++ b/src/module/DiffusionSDE.cpp @@ -1,14 +1,5 @@ #include "crpropa/module/DiffusionSDE.h" -#include - -#include -#include -#include -#include -#include - -#include "kiss/logger.h" using namespace crpropa; @@ -107,7 +98,7 @@ void DiffusionSDE::process(Candidate *candidate) const { Vector3d DirOut = Vector3d(0.); - double propTime = TStep * pow(h, 0.5) / c_light; + double propTime = TStep * sqrt(h) / c_light; size_t counter = 0; double r=42.; //arbitrary number larger than one @@ -128,7 +119,7 @@ void DiffusionSDE::process(Candidate *candidate) const { //std::cout << "counter = " << counter << "\n"; size_t stepNumber = pow(2, counter-1); - double allowedTime = TStep * pow(h, 0.5) / c_light / stepNumber; + double allowedTime = TStep * sqrt(h) / c_light / stepNumber; Vector3d Start = PosIn; Vector3d PosOut = Vector3d(0.); Vector3d PosErr = Vector3d(0.); @@ -182,7 +173,7 @@ void DiffusionSDE::process(Candidate *candidate) const { } // Integration of the SDE with a Mayorama-Euler-method - Vector3d PO = PosOut + LinProp + (NVec * NStep + BVec * BStep) * pow(h, 0.5) ; + Vector3d PO = PosOut + LinProp + (NVec * NStep + BVec * BStep) * sqrt(h) ; // Throw error message if something went wrong with propagation. // Deactivate candidate. @@ -228,13 +219,13 @@ void DiffusionSDE::process(Candidate *candidate) const { /* const std::string AL = "arcLength"; if (candidate->hasProperty(AL) == false){ - double arcLen = (TStep + NStep + BStep) * pow(h, 0.5); + double arcLen = (TStep + NStep + BStep) * sqrt(h); candidate->setProperty(AL, arcLen); return; } else { double arcLen = candidate->getProperty(AL); - arcLen += (TStep + NStep + BStep) * pow(h, 0.5); + arcLen += (TStep + NStep + BStep) * sqrt(h); candidate->setProperty(AL, arcLen); } */ @@ -259,8 +250,8 @@ void DiffusionSDE::tryStep(const Vector3d &PosIn, Vector3d &POut, Vector3d &PosE BField = magneticField->getField(y_n, z); } catch (std::exception &e) { - std::cerr << "DiffusionSDE: Exception in getField." << std::endl; - std::cerr << e.what() << std::endl; + KISS_LOG_ERROR << "DiffusionSDE: Exception in magneticField::getField.\n" + << e.what(); } k[i] = BField.getUnitVector() * c_light; @@ -277,8 +268,8 @@ void DiffusionSDE::driftStep(const Vector3d &Pos, Vector3d &LinProp, double h) c AdvField = advectionField->getField(Pos); } catch (std::exception &e) { - std::cerr << "DiffusionSDE: Exception in getField." << std::endl; - std::cerr << e.what() << std::endl; + KISS_LOG_ERROR << "DiffusionSDE: Exception in advectionField::getField.\n" + << e.what(); } LinProp += AdvField * h; return; diff --git a/test/testAdiabaticCooling.cpp b/test/testAdiabaticCooling.cpp index ac122825f..ba6b872bc 100644 --- a/test/testAdiabaticCooling.cpp +++ b/test/testAdiabaticCooling.cpp @@ -5,7 +5,7 @@ #include "crpropa/module/AdiabaticCooling.h" #include "gtest/gtest.h" -#include +//#include namespace crpropa {