Line data    Source code

       1              : #include "crpropa/module/MomentumDiffusion.h"
       2              : 
       3              : using namespace crpropa;
       4              : 
       5            1 : ConstantMomentumDiffusion::ConstantMomentumDiffusion(double Dpp) {
       6            1 :         setLimit(0.1);
       7            1 :         setDpp(Dpp);
       8            1 : }
       9              : 
      10            0 : ConstantMomentumDiffusion::ConstantMomentumDiffusion(double Dpp, double limit) {
      11            0 :         setLimit(limit);
      12            0 :         setDpp(Dpp);
      13            0 : }
      14              : 
      15            0 : void ConstantMomentumDiffusion::process(Candidate *c) const {
      16            0 :         double rig = c->current.getRigidity();
      17            0 :         if (std::isinf(rig)) {
      18              :                 return; // Only charged particles
      19              :         }
      20              :         
      21            0 :         double p = c->current.getEnergy() / c_light; // Note we use E=p/c (relativistic limit)
      22            0 :         double dt = c->getCurrentStep() / c_light;
      23              :         
      24            0 :         double eta =  Random::instance().randNorm();
      25            0 :         double domega = eta * sqrt(dt);
      26              :         
      27            0 :         double AScal = calculateAScalar(p);
      28            0 :         double BScal = calculateBScalar();
      29              : 
      30            0 :         double dp = AScal * dt + BScal * domega;
      31            0 :         c->current.setEnergy((p + dp) * c_light);
      32              :         
      33            0 :         c->limitNextStep(limit * p / AScal * c_light);
      34              : }
      35              : 
      36            0 : double ConstantMomentumDiffusion::calculateAScalar(double p) const {
      37            0 :         double a = + 2. / p * Dpp;
      38            0 :         return a; 
      39              : }
      40              : 
      41            0 : double ConstantMomentumDiffusion::calculateBScalar() const {
      42            0 :         double b = sqrt(2 * Dpp);
      43            0 :         return b;
      44              : }
      45              : 
      46            1 : void ConstantMomentumDiffusion::setDpp(double d) {
      47            1 :         if (d < 0 )
      48              :                 throw std::runtime_error(
      49            0 :                                 "ConstantMomentumDiffusion: Dpp must be non-negative");
      50            1 :         Dpp = d;
      51            1 : }
      52              : 
      53            2 : void ConstantMomentumDiffusion::setLimit(double l) {
      54            2 :         limit = l;
      55            2 : }
      56              : 
      57            0 : double ConstantMomentumDiffusion::getDpp() const {
      58            0 :         return Dpp;
      59              : }
      60              : 
      61            0 : double ConstantMomentumDiffusion::getLimit() const {
      62            0 :         return limit;
      63              : }
      64              : 
      65            0 : std::string ConstantMomentumDiffusion::getDescription() const {
      66            0 :         std::stringstream s;
      67            0 :         s << "limit: " << limit << "\n";
      68            0 :         s << "Dpp: " << Dpp / (meter * meter / second) << " m^2/s";
      69              : 
      70            0 :         return s.str();
      71            0 : }