     double g_[3];
     double kz_;
 };
+
+/**
+ * @brief Makes a random bath in the RZ plane
+ */
+struct BathRZ{
+    /**
+    * @param Rm2 squared Number of Fourier modes in R direction
+    * @param Zm2 squared Number of Fourier modes in Z direction
+    * @param RZm \f[RZm= (Rm2+Zm2)^{1/2};\f]
+    * @param norm normalisation factor \f[norm= (2/(Rm* Zm))^{1/2}\f]
+    * @param tpi \f[tpi= 2 \pi\f]
+    * @param tpi2 \f[tpi= 4 \pi^2\f]
+    * @param k \f[k=(k_R^2+k_Z^2)^{1/2}\f]
+    * @param sqEkvec  \f[ E_k^{1/2}\f]
+    * @param unif1 uniform real random variable between [0,2 pi]
+    * @param unif2 uniform real random variable between [0,2 pi]
+    * @param normal1 normal random variable with mean=0 and standarddeviation=1
+    * @param normal2 normal random variable with mean=0 and standarddeviation=1
+    * @param normalamp \f[normalamp= (normal1^2+normal2^2)^{1/2}\f]
+    * @param normalphase \f[normalphase= arctan(normal2/normal1)\f]
+    */
+  double Rm2, Zm2, RZm;
+  double norm,tpi,tpi2;    
+  std::vector<double> kvec;
+  std::vector<double> sqEkvec;
+  std::vector<double> unif1, unif2, normal1,normal2,normalamp,normalphase;
+      /**
+     * @brief Functor returning a random field in the RZ-plane or in the first RZ-plane
+     *
+     * @param Rm Number of Fourier modes in R direction
+     * @param Zm Number of Fourier modes in Z direction
+     * @param Nz Number of planes in phi direction
+     * @param R_min Minimal R (in units of rho_s)
+     * @param Z_min Minimal Z (in units of rho_s)
+     * @param gamma exponent of the energy function \f[E_k=(k/(k+_k0)^2)^\gamma\f](typical around 30)
+     * @param eddysize \f[k_0=2*\pi*eddysize/XYm \f]
+     * @param amp Amplitude
+     */  
+  BathRZ( unsigned Rm, unsigned Zm, unsigned Nz, double R_min, double Z_min, double gamma, double eddysize, double amp) : Rm_(Rm), Zm_(Zm), Nz_(Nz),R_min_(R_min), Z_min_(Z_min), gamma_(gamma), eddysize_(eddysize) , amp_(amp) {
+    Rm2=(double)Rm_*(double)Rm_;
+    Zm2=(double)Zm_*(double)Zm_;
+    RZm= sqrt(Rm2+Zm2);
+    norm=sqrt(2./(double)Rm_/(double)Zm_); 
+    tpi=2.*M_PI; tpi2=tpi*tpi;
+    double k0= tpi*eddysize_/RZm;
+    double Rmh = Rm_/2.;
+    double Zmh = Zm_/2.;
+    
+    kvec.resize(Rm_*Zm_);
+    sqEkvec.resize(Rm_*Zm_);
+    unif1.resize(Rm_*Zm_);
+    unif2.resize(Rm_*Zm_);
+    normal1.resize(Rm_*Zm_);
+    normal2.resize(Rm_*Zm_);
+    normalamp.resize(Rm_*Zm_);
+    normalphase.resize(Rm_*Zm_);
+    
+    thrust::random::minstd_rand generator;
+    thrust::random::normal_distribution<float> ndistribution;
+    thrust::random::uniform_real_distribution<float> udistribution(0.0,tpi);
+    
+    for (unsigned j=1;j<=Zm_;j++)
+    {
+      double kZ2=tpi2*(j-Zmh)*(j-Zmh)/(Zm2);
+      for (unsigned i=1;i<=Rm_;i++)
+      {
+        double kR2=tpi2*(i-Rmh)*(i-Rmh)/(Rm2);
+        int z=(j-1)*(Rm_)+(i-1);
+        kvec[z]= sqrt(kR2 + kZ2);  //radial k number
+        sqEkvec[z]=pow(kvec[z]*4.*k0/(kvec[z]+k0)/(kvec[z]+k0),gamma_/2.); //Energie in k space with max at 1.
+        unif1[z]=cos(udistribution(generator));
+        unif2[z]=sin(udistribution(generator));
+        normal1[z]=ndistribution(generator);
+        normal2[z]=ndistribution(generator);
+        normalamp[z]=sqrt(normal1[z]*normal1[z]+normal2[z]*normal2[z]);
+        normalphase[z]=atan2(normal2[z],normal1[z]);
+      }
+    }
+    
+  }
+      /**
+     * @brief Return the value of the Bath
+     *
+     * @param R R - coordinate
+     * @param Z Z - coordinate
+     *
+     */
+  double operator()(double R, double Z)
+  { 
+    double f, RZphasecos;
+    double  RR, ZZ;
+    RR=R-R_min_;
+    ZZ=Z-Z_min_;
+    f=0;
+    for (unsigned j=0;j<Zm_;j++)
+      {
+      for (unsigned i=0;i<Rm_;i++)
+       {
+        int z=j*Rm_+i;
+        RZphasecos= RR*unif1[z]+ZZ*unif2[z];        
+        f+= sqEkvec[z]*normalamp[z]*cos(kvec[z]*RZphasecos+normalphase[z]); 
+      }      
+    }
+    return amp_*norm*f;
+    
+  }
+    /**
+     * @brief Return the value of the Bath for first phi plane
+     *
+     * @param R R - coordinate
+     * @param Z Z - coordinate
+     * @param phi phi - coordinate
+     *
+     */
+  double operator()(double R, double Z, double phi) { 
+    double f, RZphasecos;
+    double  RR, ZZ;
+    RR=R-R_min_;
+    ZZ=Z-Z_min_;
+    f=0;
+    if (phi== M_PI/Nz_)
+    {
+        for (unsigned j=0;j<Zm_;j++)
+        {
+            for (unsigned i=0;i<Rm_;i++)
+            {
+                int z=(j)*(Rm_)+(i);
+                RZphasecos= RR*unif1[z]+ZZ*unif2[z];        
+                f+= sqEkvec[z]*normalamp[z]*cos(kvec[z]*RZphasecos+normalphase[z]); 
+            }      
+        }
+    return amp_*norm*f;
+    }
+    else {
+    return 0.;
+    }
+  }
+  private:
+  unsigned Rm_,Zm_,Nz_;
+  double gamma_, eddysize_, R_min_, Z_min_;
+  double amp_;
+};
 /**
  * @brief Exponential
  *

#define _DG_FUNCTORS_CUH_

#include <cmath>

namespace dg
{
 

struct AbsMax
{
    /**
     * @brief Return the asbolute maximum
     *
     * @param x left value
     * @param y right value
