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