core/wave.hh

#ifndef VTESTBED_CORE_WAVE_HH
#define VTESTBED_CORE_WAVE_HH

#include <cmath>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/grid.hh>
#include <vtestbed/core/units.hh>

namespace vtb {

    namespace core {

        namespace bits {

            template <class T, int N>
            struct Wave_number;

            template <class T>
            struct Wave_number<T,1> {
                static inline blitz::TinyVector<T,1>
                wave_number(Length<T> length, Direction<T> angle) {
                    return 2*T{M_PI}/length;
                }
            };

            template <class T>
            struct Wave_number<T,2> {
                static inline blitz::TinyVector<T,2>
                wave_number(Length<T> length, Direction<T> angle) {
                    const T a = 2*T{M_PI}/length;
                    return {a*std::cos(angle), a*std::sin(angle)};
                }
            };

        }

        template <class T, int N>
        inline blitz::TinyVector<T,N>
        wave_number(Length<T> length, Direction<T> angle) {
            return bits::Wave_number<T,N>::wave_number(length, angle);
        }

        template <class T, int N>
        class alignas(16) Wave_base {

        public:
            typedef T value_type;
            typedef blitz::TinyVector<T,N-1> vec;

        protected:
            T _amplitude{1};
            vec _k{T{0.05}};
            T _angfrequency{1};

        public:

            Wave_base() = default;
            Wave_base(const Wave_base&) = default;
            Wave_base& operator=(const Wave_base&) = default;

            inline explicit
            Wave_base(Amplitude<T> amplitude, const vec& k, Frequency<T> freq):
            _amplitude{amplitude},
            _k{k},
            _angfrequency{freq} {}

            inline explicit
            Wave_base(
                Amplitude<T> amplitude,
                Length<T> length,
                Direction<T> direction,
                Frequency<T> velocity
            ):
            _amplitude{amplitude},
            _k{::vtb::core::wave_number<T,N-1>(length, direction)},
            _angfrequency{velocity} {}

            inline T
            amplitude() const noexcept {
                return this->_amplitude;
            }

            inline const vec&
            wave_number() const noexcept {
                return this->_k;
            }

            inline T
            angular_frequency() const noexcept {
                return this->_angfrequency;
            }

            inline vec
            length() const noexcept {
                return T{2}*T{M_PI} / this->wave_number();
            }

            inline T
            scalar_length() const noexcept {
                return T{2}*T{M_PI} / blitz::length(this->wave_number());
            }

            inline T
            direction() const noexcept {
                return std::atan2(this->_k(1), this->_k(0));
            }

            inline T
            period() const noexcept {
                return T{2}*T{M_PI} / this->angular_frequency();
            }

            inline vec
            phase_velocity() const noexcept {
                return this->angular_frequency() / this->wave_number();
            }

            inline constexpr static int
            dimensions() noexcept {
                return N;
            }

        };

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

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

        public:

            No_wave() = default;
            No_wave(const No_wave&) = default;
            No_wave& operator=(const No_wave&) = default;

            inline T
            operator()(const vec&, T) const noexcept {
                return T{0};
            }

        };

        template <class T, int N, class Wave>
        inline T
        courant_number_sum(const Wave& wave, const Grid<T,N>& grid) {
            static_assert(N == 3, "bad N");
            const auto& dx = grid.select(1,2).delta();
            return sum(wave.max_velocity() / dx);
        }

        template <class T, int N, class Wave>
        inline T
        courant_number(const Wave& wave, const Grid<T,N>& grid) {
            static_assert(N == 3, "bad N");
            const T& dt = grid.delta(0);
            return dt * courant_number_sum<T,N,Wave>(wave, grid);
        }

    }

}

#endif // vim:filetype=cpp