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            1 : ConstantMomentumDiffusion::ConstantMomentumDiffusion(double Dpp, double limit) {
      11            1 :         setLimit(limit);
      12            1 :         setDpp(Dpp);
      13            1 : }
      14              : 
      15            1 : void ConstantMomentumDiffusion::process(Candidate *c) const {
      16            1 :         double rig = c->current.getRigidity();
      17            1 :         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            1 : double ConstantMomentumDiffusion::calculateAScalar(double p) const {
      37            1 :         double a = + 2. / p * Dpp;
      38            1 :         return a; 
      39              : }
      40              : 
      41            1 : double ConstantMomentumDiffusion::calculateBScalar() const {
      42            1 :         double b = sqrt(2 * Dpp);
      43            1 :         return b;
      44              : }
      45              : 
      46            2 : void ConstantMomentumDiffusion::setDpp(double d) {
      47            2 :         if (d < 0 )
      48              :                 throw std::runtime_error(
      49            0 :                                 "ConstantMomentumDiffusion: Dpp must be non-negative");
      50            2 :         Dpp = d;
      51            2 : }
      52              : 
      53            3 : void ConstantMomentumDiffusion::setLimit(double l) {
      54            3 :         limit = l;
      55            3 : }
      56              : 
      57            2 : double ConstantMomentumDiffusion::getDpp() const {
      58            2 :         return Dpp;
      59              : }
      60              : 
      61            2 : double ConstantMomentumDiffusion::getLimit() const {
      62            2 :         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 : }