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            2 :     double time = c->getTime();
      20              :         
      21              :         double Div = 0.;        
      22              :         try {
      23            2 :                 Div +=  advectionField->getDivergence(pos, time);
      24              :         } 
      25            0 :         catch (std::exception &e) {
      26            0 :                 KISS_LOG_ERROR  << "AdiabaticCooling: Exception in getDivergence.\n" 
      27            0 :                                 << e.what();
      28            0 :         }
      29              :         
      30            2 :         double dEdt = -E / 3. * Div;    // cooling due to advection -p/3 * div(V_wind)
      31              :                                         // (see e.g. Kopp et al. Computer Physics Communication 183
      32              :                                         // (2012) 530-542)
      33            2 :         double dt = c->getCurrentStep() / c_light;
      34            2 :         double dE = dEdt * dt;
      35              :         
      36            2 :         c->current.setEnergy(E + dE);
      37            2 :         if (dEdt==0) {
      38              :                 return;
      39              :         }       
      40            1 :         c->limitNextStep(limit * E / fabs(dEdt) *c_light);
      41              : }
      42              : 
      43            3 : void AdiabaticCooling::setLimit(double l) {
      44            3 :         limit = l;
      45            3 : }
      46              : 
      47            1 : double AdiabaticCooling::getLimit() const {
      48            1 :         return limit;
      49              : }
      50              :         
      51              :         
      52              : 
      53              : 
      54              : 
      55              : } // end namespace crpropa