Line data    Source code

       1             : #ifndef CRPROPA_SIMPLEGRIDTURBULENCE_H
       2             : #define CRPROPA_SIMPLEGRIDTURBULENCE_H
       3             : 
       4             : #ifdef CRPROPA_HAVE_FFTW3F
       5             : 
       6             : #include "crpropa/magneticField/turbulentField/GridTurbulence.h"
       7             : 
       8             : #include "kiss/logger.h"
       9             : #include "kiss/string.h"
      10             : 
      11             : namespace crpropa {
      12             : /**
      13             :  * \addtogroup MagneticFields
      14             :  * @{
      15             :  */
      16             : 
      17             : /**
      18             :  @class SimpleTurbulenceSpectrum
      19             :  @brief Defines the energy spectrum of simple power-law turbulence
      20             :  */
      21             : class SimpleTurbulenceSpectrum : public TurbulenceSpectrum {
      22             :         const int constScaleBendover = 1000;    // define the bandover scale as 1000 * lMax to ensure k * lBendover >> 1. The bendover scale is necessary for the implementation of PlaneWaveTurbulence. 
      23             :   public:
      24             :         /**
      25             :          * The bend-over scale is set to 1000 times lMax to ensure to be in the inertial range. This should not be changed.
      26             :          @param Brms            Root mean square field strength for generated field
      27             :          @param lMin            Minimum physical scale of the turbulence
      28             :          @param lMax            Maximum physical scale of the turbulence
      29             :          @param sIndex          Spectral index of the energy spectrum in the inertial range
      30             :         */
      31             :         SimpleTurbulenceSpectrum(double Brms, double lMin, double lMax,
      32             :                                  double sIndex = 5. / 3)
      33           2 :             : TurbulenceSpectrum(Brms, lMin, lMax, constScaleBendover * lMax, sIndex, 0) {}
      34           0 :         ~SimpleTurbulenceSpectrum() {}
      35             : 
      36             :         /**
      37             :         General energy spectrum for synthetic turbulence models
      38             :         */
      39      206961 :         double energySpectrum(double k) const {
      40      206961 :                 return std::pow(k, -getSindex() - 2);
      41             :         }
      42             : 
      43             :         /**
      44             :             @brief       compute the magnetic field coherence length according to
      45             :            the formula in  Harari et al. JHEP03(2002)045
      46             :             @return Lc   coherence length of the magnetic field
      47             :         */
      48           0 :         double getCorrelationLength() const {
      49           0 :                 return turbulentCorrelationLength(getLmin(), getLmax(),
      50           0 :                                                   getSindex());
      51             :         }
      52           1 :         static double turbulentCorrelationLength(double lMin, double lMax,
      53             :                                                  double s) {
      54           1 :                 double r = lMin / lMax;
      55           1 :                 return lMax / 2 * (s - 1) / s * (1 - pow(r, s)) / (1 - pow(r, s - 1));
      56             :         }
      57             : };
      58             : 
      59             : /**
      60             :  @class SimpleGridTurbulence
      61             :  @brief Turbulent grid-based magnetic field with a simple power-law spectrum
      62             :  */
      63           3 : class SimpleGridTurbulence : public GridTurbulence {
      64             :   public:
      65             :         /**
      66             :          Create a random initialization of a turbulent field.
      67             :          @param spectrum    TurbulenceSpectrum instance to define the spectrum of
      68             :          turbulence
      69             :          @param gridProp        GridProperties instance to define the underlying grid
      70             :          @param seed     Random seed
      71             :          */
      72             :         SimpleGridTurbulence(const SimpleTurbulenceSpectrum &spectrum,
      73             :                              const GridProperties &gridProp, unsigned int seed = 0);
      74             : 
      75             :         static void initTurbulence(ref_ptr<Grid3f> grid, double Brms, double lMin,
      76             :                                    double lMax, double alpha, int seed);
      77             : };
      78             : 
      79             : // Compatibility with old functions from GridTurbulence:
      80             : 
      81             : /** Analytically calculate the correlation length of the simple model turbulent
      82             :  * field */
      83           1 : inline double turbulentCorrelationLength(double lMin, double lMax,
      84             :                                          double alpha = -11 / 3.) {
      85           2 :         KISS_LOG_WARNING
      86             :             << "turbulentCorrelationLength is deprecated and will be "
      87             :                "removed in the future. Replace it with a more appropriate "
      88           1 :                "turbulent field model and call getCorrelationLength().";
      89           1 :         return SimpleTurbulenceSpectrum::turbulentCorrelationLength(lMin, lMax,
      90           1 :                                                                     -alpha - 2);
      91             : }
      92             : 
      93             : /**
      94             :  Create a random initialization of a turbulent field.
      95             :  @param grid    grid on which the turbulence is calculated
      96             :  @param lMin    Minimum wavelength of the turbulence
      97             :  @param lMax    Maximum wavelength of the turbulence
      98             :  @param alpha   Power law index of <B^2(k)> ~ k^alpha (alpha = -11/3 corresponds
      99             :  to a Kolmogorov spectrum)
     100             :  @param Brms    RMS field strength
     101             :  @param seed    Random seed
     102             :  */
     103           4 : inline void initTurbulence(ref_ptr<Grid3f> grid, double Brms, double lMin,
     104             :                            double lMax, double alpha = -11 / 3., int seed = 0) {
     105           8 :         KISS_LOG_WARNING
     106             :             << "initTurbulence is deprecated and will be removed in the future. "
     107           4 :                "Replace it with a more appropriate turbulent field model instance.";
     108           4 :         SimpleGridTurbulence::initTurbulence(grid, Brms, lMin, lMax, alpha, seed);
     109           1 : }
     110             : 
     111             : /** @}*/
     112             : } // namespace crpropa
     113             : 
     114             : #endif // CRPROPA_HAVE_FFTW3F
     115             : 
     116             : #endif // CRPROPA_SIMPLEGRIDTURBULENCE_H