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           0 :                 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 : }