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, 1000 * 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