stokes_wave.hh

#ifndef VTESTBED_CORE_STOKES_WAVE_HH
#define VTESTBED_CORE_STOKES_WAVE_HH

#include <stdexcept>

#include <vtestbed/base/constants.hh>
#include <vtestbed/core/wave.hh>

namespace vtb {

    namespace core {

        template <class T, int N>
        inline Frequency<T>
        stokes_wave_angular_frequency(
            const blitz::TinyVector<T,N>& k,
            Amplitude<T> amplitude
        ) {
            using std::sqrt;
            constexpr auto g = base::constants<T>::g();
            const T k_length{length(k)};
            const T steepness = amplitude * k_length;
            return Frequency<T>{sqrt(g*k_length*(T{1} +steepness*steepness))};
        }

        template <class T, int N>
        class Propagating_stokes_wave: public Wave_base<T,N> {

        public:
            using typename Wave_base<T,N>::vec;

        public:

            Propagating_stokes_wave() = default;
            ~Propagating_stokes_wave() = default;
            Propagating_stokes_wave(const Propagating_stokes_wave&) = default;
            Propagating_stokes_wave& operator=(const Propagating_stokes_wave&) = default;

            inline
            Propagating_stokes_wave(Amplitude<T> amplitude, const vec& k):
            Wave_base<T,N>(
                amplitude,
                k,
                stokes_wave_angular_frequency(k, amplitude)
            ) {}

            inline
            Propagating_stokes_wave(
                Amplitude<T> amplitude,
                Length<T> length,
                Direction<T> direction
            ):
            Propagating_stokes_wave(
                amplitude,
                wave_number<T,N-1>(length, direction)
            ) {}


            inline T
            operator()(const vec& x, T t) const noexcept {
                using blitz::dot;
                using std::cos;
                const T steepness = this->amplitude() * length(this->wave_number());
                const T theta = dot(this->wave_number(), x) - this->angular_frequency()*t;
                return this->amplitude() * (
                    cos(theta)
                    + T(0.5) * steepness * cos(T(2)*theta)
                    + (T(3)/T(8)) * steepness * steepness * cos(T(3)*theta)
                );
            }

            inline Grid<T,3>
            acf_grid() const {
                throw std::runtime_error("not implemented");
            }

        };

    }

}

#endif // vim:filetype=cpp