#ifndef OU_VORTICITY_EQUATION_HPP
#define OU_VORTICITY_EQUATION_HPP

#include <cstdlib>
#include "vorticity_equation.hpp"
#include "ornstein_uhlenbeck_process.hpp"


// vorticity_equation<double, FFTW> *test;

template <typename rnumber, field_backend be>
class ou_vorticity_equation : public vorticity_equation<rnumber, be>
{
    public:

      std::string ou_forcing_type;

      ornstein_uhlenbeck_process<rnumber,be> *ou;
      ou_vorticity_equation(
        const char *NAME,
        int nx,
        int ny,
        int nz,
        double ou_kmin,
        double ou_kmax,
        double ou_energy_amplitude,
        double ou_gamma_factor,
        double DKX = 1.0,
        double DKY = 1.0,
        double DKZ = 1.0,
        unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE):
          vorticity_equation<rnumber,be>(
            NAME, nx, ny, nz,
            DKX,DKY,DKZ, FFTW_PLAN_RIGOR){
              this->ou = new ornstein_uhlenbeck_process<rnumber,be>(
                NAME,
                nx, ny, nz,
                ou_kmin, ou_kmax, ou_energy_amplitude,ou_gamma_factor,
                DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
              this->ou_forcing_type = "replace";
            }

        ~ou_vorticity_equation();

        void omega_nonlin(int src);
        void step(double dt);

        void add_ou_forcing_velocity(void);
        void add_ou_forcing_vorticity(void);

};


#endif