Line data    Source code

       1             : #include "crpropa/module/Acceleration.h"
       2             : #include <crpropa/Common.h>
       3             : #include <crpropa/Random.h>
       4             : #include <cmath>
       5             : 
       6             : namespace crpropa {
       7             : 
       8           0 : AbstractAccelerationModule::AbstractAccelerationModule(double _stepLength)
       9           0 :         : crpropa::Module(), stepLength(_stepLength) {}
      10             : 
      11             : 
      12           0 : void AbstractAccelerationModule::add(StepLengthModifier *modifier) {
      13           0 :         modifiers.push_back(modifier);
      14           0 : }
      15             : 
      16             : 
      17           0 : void AbstractAccelerationModule::scatter(
      18             :         crpropa::Candidate *candidate,
      19             :         const crpropa::Vector3d &scatter_center_velocity) const {
      20             :         // particle momentum in lab frame
      21           0 :         const double E = candidate->current.getEnergy();
      22           0 :         const crpropa::Vector3d p = candidate->current.getMomentum();
      23             : 
      24             :         // transform to rest frame of scatter center (p: prime)
      25           0 :         const double beta = scatter_center_velocity.getR() / crpropa::c_light;
      26           0 :         const double gamma = 1. / sqrt(1 - beta * beta);
      27           0 :         const double Ep = gamma * (E - scatter_center_velocity.dot(p));
      28             :         const crpropa::Vector3d pp = (p - scatter_center_velocity* E /
      29             :                 (crpropa::c_light * crpropa::c_light)) * gamma;
      30             : 
      31             :         // scatter into random direction
      32           0 :         const crpropa::Vector3d pp_new = crpropa::Random::instance().randVector() * pp.getR();
      33             : 
      34             :         // transform back
      35           0 :         const double E_new = gamma * (Ep + scatter_center_velocity.dot(pp_new));
      36             :         const crpropa::Vector3d p_new = (pp_new + scatter_center_velocity * Ep /
      37             :                 (crpropa::c_light * crpropa::c_light)) * gamma;
      38             : 
      39             :         // update candidate properties
      40           0 :         candidate->current.setEnergy(E_new);
      41           0 :         candidate->current.setDirection(p_new / p_new.getR());
      42           0 : }
      43             : 
      44             : 
      45           0 : void AbstractAccelerationModule::process(crpropa::Candidate *candidate) const {
      46           0 :         double currentStepLength = stepLength;
      47           0 :         for (auto m : modifiers) {
      48           0 :                 currentStepLength = m->modify(currentStepLength, candidate);
      49             :         }
      50             : 
      51           0 :         double step = candidate->getCurrentStep();
      52           0 :         while (step > 0) {
      53           0 :                 double randDistance = -1. * log(crpropa::Random::instance().rand()) * currentStepLength;
      54             : 
      55           0 :                 if (step < randDistance) {
      56           0 :                         candidate->limitNextStep(0.1 * currentStepLength);
      57           0 :                         return;
      58             :                 }
      59           0 :                 scatter(candidate, scatterCenterVelocity(candidate));
      60           0 :                 step -= randDistance;
      61             :         }
      62             : }
      63             : 
      64             : 
      65           0 : SecondOrderFermi::SecondOrderFermi(double scatterVelocity, double stepLength,
      66           0 :                                                                    unsigned int sizeOfPitchangleTable)
      67             :         : AbstractAccelerationModule(stepLength),
      68           0 :           scatterVelocity(scatterVelocity) {
      69           0 :         setDescription("SecondOrderFermi Acceleration");
      70           0 :         angle.resize(sizeOfPitchangleTable);
      71           0 :         angleCDF.resize(sizeOfPitchangleTable);
      72             : 
      73             :         // have a discretized table of beamed pitch angles
      74           0 :         for (unsigned int i =0; i < sizeOfPitchangleTable; i++) {
      75           0 :                 angle[i] = i * M_PI / (sizeOfPitchangleTable-1);
      76           0 :                 angleCDF[i] = (angle[i] +scatterVelocity / crpropa::c_light * sin(angle[i])) / M_PI;
      77             :         }
      78           0 : }
      79             : 
      80             : 
      81           0 : crpropa::Vector3d SecondOrderFermi::scatterCenterVelocity(crpropa::Candidate *candidate) const
      82             : {
      83           0 :         size_t idx = crpropa::closestIndex(crpropa::Random::instance().rand(), angleCDF);
      84           0 :         crpropa::Vector3d rv = crpropa::Random::instance().randVector();
      85           0 :         crpropa::Vector3d rotationAxis = candidate->current.getDirection().cross(rv);
      86             : 
      87           0 :         rv = candidate->current.getDirection().getRotated(rotationAxis, M_PI - angle[idx]);
      88           0 :         return rv * scatterVelocity;
      89             : }
      90             : 
      91             : 
      92           0 : DirectedFlowOfScatterCenters::DirectedFlowOfScatterCenters(
      93           0 :         const Vector3d &scatterCenterVelocity)
      94           0 :         : __scatterVelocity(scatterCenterVelocity) {}
      95             : 
      96             : 
      97           0 : double DirectedFlowOfScatterCenters::modify(double steplength, Candidate* candidate)
      98             : {
      99           0 :         double directionModifier = (-1. * __scatterVelocity.dot(candidate->current.getDirection()) + c_light) / c_light;
     100           0 :         return steplength / directionModifier;
     101             : }
     102             : 
     103             : 
     104           0 : DirectedFlowScattering::DirectedFlowScattering(
     105           0 :         crpropa::Vector3d scatterCenterVelocity, double stepLength)
     106           0 :         : __scatterVelocity(scatterCenterVelocity),
     107           0 :           AbstractAccelerationModule(stepLength) {
     108             : 
     109             :         // In a directed field of scatter centers, the probability to encounter a
     110             :         // scatter center depends on the direction of the candidate.
     111           0 :         StepLengthModifier *mod = new DirectedFlowOfScatterCenters(__scatterVelocity);
     112           0 :         this->add(mod);
     113           0 : }
     114             : 
     115             : 
     116           0 : crpropa::Vector3d DirectedFlowScattering::scatterCenterVelocity(
     117             :         crpropa::Candidate *candidate) const { // does not depend on candidate here.
     118           0 :         return __scatterVelocity;
     119             : }
     120             : 
     121             : 
     122           0 : QuasiLinearTheory::QuasiLinearTheory(double referenecEnergy,
     123             :                                                                          double turbulenceIndex,
     124           0 :                                                                          double minimumRigidity)
     125           0 :         : __referenceEnergy(referenecEnergy), __turbulenceIndex(turbulenceIndex),
     126           0 :           __minimumRigidity(minimumRigidity) {}
     127             : 
     128             : 
     129           0 : double QuasiLinearTheory::modify(double steplength, Candidate* candidate)
     130             : {
     131           0 :         if (candidate->current.getRigidity() < __minimumRigidity)
     132             :         {
     133           0 :                 return steplength * std::pow(__minimumRigidity /
     134           0 :                         (__referenceEnergy / eV), 2. - __turbulenceIndex);
     135             :         }
     136             :         else
     137             :         {
     138           0 :                 return steplength * std::pow(candidate->current.getRigidity() /
     139           0 :                         (__referenceEnergy / eV), 2. - __turbulenceIndex);
     140             :         }
     141             : }
     142             : 
     143             : 
     144           0 : ParticleSplitting::ParticleSplitting(Surface *surface, int      crossingThreshold, 
     145           0 :         int numberSplits, double minWeight, std::string counterid)
     146           0 :         : surface(surface), crossingThreshold(crossingThreshold),
     147           0 :           numberSplits(numberSplits), minWeight(minWeight), counterid(counterid){};
     148             : 
     149           0 : void ParticleSplitting::process(Candidate *candidate) const {
     150             :         const double currentDistance =
     151           0 :                 surface->distance(candidate->current.getPosition());
     152             :         const double previousDistance =
     153           0 :                 surface->distance(candidate->previous.getPosition());
     154             : 
     155           0 :         if (currentDistance * previousDistance > 0)
     156             :                 // candidate remains on the same side
     157             :                 return;
     158             : 
     159           0 :         if (candidate->getWeight() < minWeight)
     160             :                 return;
     161             : 
     162             :         int num_crossings = 1;
     163           0 :         if (candidate->hasProperty(counterid))
     164           0 :                 num_crossings = candidate->getProperty(counterid).toInt32() + 1;
     165           0 :         candidate->setProperty(counterid, num_crossings);
     166             : 
     167           0 :         if (num_crossings % crossingThreshold != 0)
     168             :                 return;
     169             : 
     170           0 :         candidate->updateWeight(1. / numberSplits);
     171             : 
     172           0 :         for (size_t i = 1; i < numberSplits; i++) {
     173             :                 // No recursive split as the weights of the secondaries created
     174             :                 // before the split are not affected
     175           0 :                 ref_ptr<Candidate> new_candidate = candidate->clone(false);
     176           0 :                 new_candidate->parent = candidate;
     177           0 :                 uint64_t snr = Candidate::getNextSerialNumber();
     178           0 :                 Candidate::setNextSerialNumber(snr + 1);
     179           0 :                 new_candidate->setSerialNumber(snr);
     180             :                 candidate->addSecondary(new_candidate);
     181             :         }
     182             : };
     183             : 
     184             : } // namespace crpropa