Line data    Source code

       1             : #ifndef CRPROPA_TURBULENTFIELD_H
       2             : #define CRPROPA_TURBULENTFIELD_H
       3             : 
       4             : #include "crpropa/magneticField/MagneticField.h"
       5             : #include <cmath>
       6             : 
       7             : namespace crpropa {
       8             : /**
       9             :  * \addtogroup MagneticFields
      10             :  * @{
      11             :  */
      12             : 
      13             : /**
      14             :  @class TurbulenceSpectrum
      15             :  @brief Defines the energy spectrum of turbulence parametrizied by A(k) ~ k^q /(1 + k^2)^{(s + q)/2 + 1}
      16             :  */
      17             : class TurbulenceSpectrum : public Referenced {
      18             :   private:
      19             :         const double Brms; /**< Brms value of the turbulent field (normalization) */
      20             :         const double sIndex; /**< Spectral index for the inertial range, for example
      21             :                           s=5/3 for Kolmogorov spectrum; in some parts of the code this
      22             :                           parameter is referred by alpha which is the total 3D isotropic
      23             :                           spectrum with additional k^2 and the minus sign, e.g.,
      24             :                           for Kolmogorov: alpha = -(s + 2) */
      25             :         const double qIndex; /**< Spectral index for the injection range, for
      26             :                           example q=4 for 3D homogeneous turbulence */
      27             :         const double lBendover;  /**< the bend-over scale */
      28             :         const double lMin, lMax; /**< Min and Max scale of turbulence */
      29             : 
      30             :   protected:
      31             :         /**
      32             :         Normalization for the below defined Lc
      33             :         */
      34           1 :         double spectrumNormalization() const {
      35           1 :                 return std::tgamma((sIndex + qIndex) / 2.0) /
      36           1 :                        (2.0 * std::tgamma((sIndex - 1) / 2.0) *
      37           1 :                         std::tgamma((qIndex + 1) / 2.0));
      38             :         }
      39             : 
      40             :   public:
      41             :         /**
      42             :          * @param Brms         root mean square field strength for generated field
      43             :          * @param lMin   Minimum physical scale of the turbulence
      44             :          * @param lMax   Maximum physical scale of the turbulence
      45             :          * @param lBendover        the bend-over scale
      46             :          * @param sIndex         Spectral index of the energy spectrum in the inertial range
      47             :          * @param qIndex         Spectral index of the energy spectrum in the energy range
      48             :         */
      49           9 :         TurbulenceSpectrum(double Brms, double lMin, double lMax,
      50             :                            double lBendover = 1, double sIndex = (5. / 3.),
      51             :                            double qIndex = 4)
      52           9 :             : Brms(Brms), lMin(lMin), lMax(lMax), lBendover(lBendover),
      53           9 :               sIndex(sIndex), qIndex(qIndex) {
      54           9 :                 if (lMin > lMax) {
      55           0 :                         throw std::runtime_error("TurbulenceSpectrum: lMin > lMax");
      56             :                 }
      57           9 :                 if (lMin <= 0) {
      58           0 :                         throw std::runtime_error("TurbulenceSpectrum: lMin <= 0");
      59             :                 }
      60           9 :         }
      61             : 
      62           7 :         ~TurbulenceSpectrum() {}
      63             : 
      64          31 :         double getBrms() const { return Brms; }
      65          19 :         double getLmin() const { return lMin; }
      66          19 :         double getLmax() const { return lMax; }
      67          22 :         double getLbendover() const { return lBendover; }
      68           5 :         double getSindex() const { return sIndex; }
      69           2 :         double getQindex() const { return qIndex; }
      70             :         
      71             :         /**
      72             :         General energy spectrum for synthetic turbulence models (not normalized!)
      73             :         with normalized ^k = k*lBendover
      74             :         */
      75      208128 :         virtual double energySpectrum(double k) const {
      76      208128 :                 double kHat = k * lBendover;
      77      208128 :                 return std::pow(kHat, qIndex) /
      78      208128 :                                        std::pow(1.0 + kHat * kHat,
      79      208128 :                                         (sIndex + qIndex) / 2.0 + 1.0);
      80             :         }
      81             : 
      82             :         /**
      83             :   Computes the magnetic field coherence length
      84             :         Obtained from the definition of \f$l_c = 1/B_{\rm rms}^2 \int_0^\infty dr\langleB(0)B^*(r)\rangle \f$
      85             :         Approximates the true value correctly as long as lBendover <= lMax/8 (~5%
      86             :   error) (for the true value the above integral should go from lMin to lMax)
      87             :         */
      88           0 :         virtual double getCorrelationLength() const {
      89           1 :                 return 4 * M_PI / ((sIndex + 2.0) * sIndex) * spectrumNormalization() *
      90           1 :                        lBendover;
      91             :         }
      92             : };
      93             : 
      94             : /**
      95             :  @class TurbulentField
      96             :  @brief An abstract base class for different models of turbulent magnetic fields
      97             : 
      98             :  This module provides common methods for all turbulent (synthetic) magnetic
      99             :  fields. Does not actually implement any turbulent field.
     100             :  */
     101             : class TurbulentField : public MagneticField {
     102             :   protected:
     103             :         const TurbulenceSpectrum &spectrum;
     104             : 
     105             :   public:
     106           8 :         TurbulentField(const TurbulenceSpectrum &spectrum) : spectrum(spectrum) {}
     107           0 :         virtual ~TurbulentField() {}
     108             : 
     109           0 :         double getBrms() const { return spectrum.getBrms(); }
     110           0 :         virtual double getCorrelationLength() const {
     111           0 :                 return spectrum.getCorrelationLength();
     112             :         }
     113             : };
     114             : 
     115             : /** @}*/
     116             : 
     117             : } // namespace crpropa
     118             : 
     119             : #endif // CRPROPA_TURBULENTFIELD_H