fixed right amplitude calc

cpp/full_code/ornstein_uhlenbeck_process.cpp

@@ -43,6 +43,8 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
     this->ou_kmin_squ = pow(ou_kmin,2);
     this->ou_kmax_squ = pow(ou_kmax,2);
     this->ou_energy_amplitude = ou_energy_amplitude;
+    this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
+
     assert(this->kk->kM2 >= this->ou_kmax_squ); //I can avoid one if statement later in the step method
     this->initialize_B();
 
@@ -72,7 +74,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
                     double k2){
         if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ)
         {
-            double g = this->gamma(k2);
+            double g = this->gamma(sqrt(k2));
             double random[6] =
             {
                 this->d(this->gen),this->d(this->gen),this->d(this->gen),
@@ -121,9 +123,11 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
                 }
                 else
                 {
+                    //sigma =
+                     //   sqrt(4*this->gamma(kabs)*this->energy(kabs)
+                       // /this->kk->nshell[(int)(kabs/this->kk->dk)]);
                     sigma =
-                        this->ou_energy_amplitude *
-                        sqrt(4*this->gamma(ksqu)*pow(kabs,-5./3.)
+                        sqrt(4*this->gamma(kabs)*this->energy(kabs)
                         /this->kk->nshell[(int)(kabs/this->kk->dk)]);
                 }
 

cpp/full_code/ornstein_uhlenbeck_process.hpp

@@ -20,6 +20,8 @@ class ornstein_uhlenbeck_process{
         double ou_kmin_squ;
         double ou_kmax_squ;
         double ou_energy_amplitude;
+        double kolmogorov_constant = 2;
+        double epsilon;
 
         field<rnumber,be,THREE> *ou_field;
         field<rnumber,be,THREE> *ou_field_vort;
@@ -47,10 +49,19 @@ class ornstein_uhlenbeck_process{
         ~ornstein_uhlenbeck_process(void);
 
         void step(double dt);
-        inline double gamma(double ksqu)
+        inline double energy(double kabs)
         {
-            return this->ou_energy_amplitude * pow(ksqu,1./3.);
+            return this->ou_energy_amplitude * pow(kabs,-5./3);
         }
+
+        inline double gamma(double kabs)
+        {
+            //return 1.5*this->epsilon/this->energy(kabs);
+            //return 1./this->energy(kabs);
+            return kabs*sqrt(this->energy(kabs));
+        }
+        
+
         void initialize_B(void);
 
         void let_converge(void);