Line data    Source code

       1             : #include "crpropa/module/AdiabaticCooling.h"
       2             : 
       3             : namespace crpropa {
       4             : 
       5           2 : AdiabaticCooling::AdiabaticCooling(ref_ptr<AdvectionField> advectionField) :
       6           2 :         advectionField(advectionField) {
       7           2 :         setLimit(0.1);
       8           2 : }
       9             : 
      10           1 : AdiabaticCooling::AdiabaticCooling(ref_ptr<AdvectionField> advectionField, double limit) :
      11           1 :         advectionField(advectionField) {
      12           1 :         setLimit(limit);
      13           1 : }
      14             : 
      15           2 : void AdiabaticCooling::process(Candidate *c) const {
      16             : 
      17           2 :         Vector3d pos = c->current.getPosition();
      18           2 :         double E = c->current.getEnergy(); // Note we use E=p/c (relativistic limit)
      19             :         
      20             :         double Div = 0.;        
      21             :         try {
      22           2 :                 Div +=  advectionField->getDivergence(pos);
      23             :         } 
      24           0 :         catch (std::exception &e) {
      25           0 :                 KISS_LOG_ERROR  << "AdiabaticCooling: Exception in getDivergence.\n" 
      26           0 :                                 << e.what();
      27           0 :         }
      28             :         
      29           2 :         double dEdt = -E / 3. * Div;    // cooling due to advection -p/3 * div(V_wind)
      30             :                                         // (see e.g. Kopp et al. Computer Physics Communication 183
      31             :                                         // (2012) 530-542)
      32           2 :         double dt = c->getCurrentStep() / c_light;
      33           2 :         double dE = dEdt * dt;
      34             :         
      35           2 :         c->current.setEnergy(E + dE);
      36           2 :         if (dEdt==0) {
      37             :                 return;
      38             :         }       
      39           1 :         c->limitNextStep(limit * E / fabs(dEdt) *c_light);
      40             : }
      41             : 
      42           3 : void AdiabaticCooling::setLimit(double l) {
      43           3 :         limit = l;
      44           3 : }
      45             : 
      46           1 : double AdiabaticCooling::getLimit() const {
      47           1 :         return limit;
      48             : }
      49             :         
      50             :         
      51             : 
      52             : 
      53             : 
      54             : } // end namespace crpropa