wave_numbers.cc

#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/math.hh>
#include <vtestbed/core/wave_numbers.hh>

namespace {

    template <class T>
    inline bool
    same_sign(T lhs, T rhs) {
        return (lhs > T(0) && rhs >= T(0)) || (lhs < T(0) && rhs <= T(0));
    }

    template <class T>
    inline void
    update(
        vtb::core::Wave_statistics<T>& stats,
        const std::vector<vtb::core::Wave_half<T>>& halves
    ) {
        for (const auto& half : halves) {
            stats.length.update(half.length);
            stats.number.update(half.number);
            stats.height.update(half.height);
        }
    }

}

template <class T>
std::ostream&
vtb::core::operator<<(std::ostream& out, const Wave_half<T>& rhs) {
    return out
        << "length=" << rhs.length
        << ",wavenumber=" << rhs.number
        << ",height=" << rhs.height;
}

template <class T>
std::ostream&
vtb::core::operator<<(std::ostream& out, const Wave_statistics<T>& rhs) {
    out << "length: " << rhs.length << '\n';
    out << "number: " << rhs.number << '\n';
    out << "height: " << rhs.height << '\n';
    return out;
}

template <class T>
std::vector<blitz::TinyVector<T,2>>
vtb::core::find_extrema(blitz::Array<T,1> func, const Grid<T, 1>& grid) {
    std::vector<blitz::TinyVector<T,2>> extrema;
    const int n = func.extent(0);
    if (n < 3) {
        return extrema;
    }
    const T dx = grid.delta(0);
    for (int i=1; i<n-1; ++i) {
        const T y0 = func(i-1);
        const T y1 = func(i);
        const T y2 = func(i+1);
        // If (x1,y1) is local extrema.
        if (same_sign(y1-y0, y1-y2)) {
            // Use parabola to find precise vertex coordinate.
            // Formula was computed in Mathematica.
            const T x1 = grid(i, 0);
            const T k = dx*(dx + T(2)*x1);
            const T a = (y0 - T(2)*y1 + y2) / (T(2)*dx*dx);
            const T b = (-a*k - y1 + y2) / dx;
            const T c = -a*(k + x1*x1) - b*(x1 + dx) + y2;
            const T vx = -b / (T(2)*a);
            const T vy = c - b*b / (T(4)*a);
            extrema.emplace_back(vx, vy);
        }
    }
    return extrema;
}

template <class T>
std::vector<vtb::core::Wave_half<T>>
vtb::core::find_wave_halves(const std::vector<blitz::TinyVector<T,2>>& extrema) {
    constexpr auto pi_times_2 = base::constants<T>::pi(2);
    using std::abs;
    std::vector<Wave_half<T>> halves;
    const int n = extrema.size();
    if (n <= 1) {
        return halves;
    }
    blitz::TinyVector<T,2> p;
    p(0) = extrema[0](0);
    p(1) = extrema[0](1);
    for (int i=1; i<n; ++i) {
        const T y0 = extrema[i-1](1);
        const T x1 = extrema[i](0);
        const T y1 = extrema[i](1);
        // We care only about big waves with trough and crest
        // located in different planes (positive and negative).
        if (same_sign(y0, y1)) {
            if (abs(p(1)) < abs(y1)) {
                p(0) = x1;
                p(1) = y1;
            }
        } else {
            const T length = T(2)*(x1 - p(0));
            const T number = pi_times_2 / length;
            const T height = abs(y1 - p(1));
            halves.emplace_back(length, number, height);
            p(0) = x1;
            p(1) = y1;
        }
    }
    return halves;
}

template <class T>
void
vtb::core::find_wave_number_range(
    const std::vector<Wave_half<T>>& halves,
    T& min_wave_number,
    T& max_wave_number
) {
    const int n = halves.size();
    for (int i=0; i<n; ++i) {
        const T wave_number = halves[i].number;
        if (wave_number < min_wave_number) {
            min_wave_number = wave_number;
        }
        if (max_wave_number < wave_number) {
            max_wave_number = wave_number;
        }
    }
}

template <class T, int N>
vtb::core::Grid<T,N>
vtb::core::wave_number_grid(const Wave_statistics_vector<T,N>& waves) {
    Grid<T,N> wngrid;
    for (int i=0; i<N; ++i) {
        wngrid.lbound(i) = waves(i).number.min();
        wngrid.ubound(i) = waves(i).number.max();
    }
    wngrid.begin() = 0;
    wngrid.end() = 2;
    for (int i=0; i<N; ++i) {
        if (!(wngrid.lbound(i) < wngrid.ubound(i))) {
            wngrid.lbound(i) = 0, wngrid.ubound(i) = 0;
        }
    }
    return wngrid;
}

template <class T>
vtb::core::Wave_statistics_vector<T,2>
vtb::core::wave_statistics(blitz::Array<T,2> surface, const Grid<T,2>& grid) {
    Wave_statistics_vector<T,2> waves;
    const auto& _ =  blitz::Range::all();
    const int n1 = surface.extent(0);
    const int n2 = surface.extent(1);
    // compute first dimension
    for (int i=0; i<n2; ++i) {
        const auto& extrema = find_extrema(surface(_,i), grid.select(0));
        update(waves(0), find_wave_halves(extrema));
    }
    // compute second dimension
    for (int i=0; i<n1; ++i) {
        const auto& extrema = find_extrema(surface(i,_), grid.select(1));
        update(waves(1), find_wave_halves(extrema));
    }
    return waves;
}

template <class T>
vtb::core::Wave_statistics<T>
vtb::core::temporal_wave_statistics(
    blitz::Array<T,1> surface,
    const Grid<T,1>& grid
) {
    Wave_statistics<T> waves;
    update(waves, find_wave_halves(find_extrema(surface, grid)));
    return waves;
}

template <class T>
vtb::core::Grid<T,2>
vtb::core::find_wave_number_range(blitz::Array<T,2> surface, const Grid<T,2>& grid) {
    using blitz::Range;
    Grid<T,2> wngrid{
        {std::numeric_limits<T>::max(), std::numeric_limits<T>::max()},
        {std::numeric_limits<T>::min(), std::numeric_limits<T>::min()},
        {2,2}
    };
    const int n1 = surface.extent(0);
    const int n2 = surface.extent(1);
    // compute first dimension
    for (int i=0; i<n2; ++i) {
        find_wave_number_range(
            find_wave_halves(find_extrema(surface(Range::all(), i), grid.select(0))),
            wngrid.lbound(0),
            wngrid.ubound(0)
        );
    }
    // compute second dimension
    for (int i=0; i<n1; ++i) {
        find_wave_number_range(
            find_wave_halves(find_extrema(surface(i, Range::all()), grid.select(1))),
            wngrid.lbound(1),
            wngrid.ubound(1)
        );
    }
    for (int i=0; i<2; ++i) {
        if (wngrid.lbound(i) > wngrid.ubound(i)) {
            wngrid.lbound(i) = T(0);
            wngrid.ubound(i) = T(0);
        }
    }
    return wngrid;
}

template std::vector<blitz::TinyVector<VTB_REAL_TYPE,2>>
vtb::core::find_extrema(
    blitz::Array<VTB_REAL_TYPE, 1> func,
    const Grid<VTB_REAL_TYPE, 1>& grid
);

template std::vector<vtb::core::Wave_half<VTB_REAL_TYPE>>
vtb::core::find_wave_halves(
    const std::vector<blitz::TinyVector<VTB_REAL_TYPE, 2>>& extrema
);

template void
vtb::core::find_wave_number_range(
    const std::vector<Wave_half<VTB_REAL_TYPE>>& halves,
    VTB_REAL_TYPE& min_wave_number,
    VTB_REAL_TYPE& max_wave_number
);


template vtb::core::Grid<VTB_REAL_TYPE,2>
vtb::core::find_wave_number_range(
    blitz::Array<VTB_REAL_TYPE, 2> surface,
    const Grid<VTB_REAL_TYPE, 2>& grid
);

template vtb::core::Grid<VTB_REAL_TYPE,2>
vtb::core::wave_number_grid(const Wave_statistics_vector<VTB_REAL_TYPE,2>&);

template vtb::core::Grid<VTB_REAL_TYPE,3>
vtb::core::wave_number_grid(const Wave_statistics_vector<VTB_REAL_TYPE,3>&);

template blitz::TinyVector<vtb::core::Wave_statistics<VTB_REAL_TYPE>,2>
vtb::core::wave_statistics(
    blitz::Array<VTB_REAL_TYPE,2> surface,
    const Grid<VTB_REAL_TYPE,2>& grid
);

template vtb::core::Wave_statistics<VTB_REAL_TYPE>
vtb::core::temporal_wave_statistics(
    blitz::Array<VTB_REAL_TYPE,1> surface,
    const Grid<VTB_REAL_TYPE,1>& grid
);

template std::ostream&
vtb::core::operator<<(std::ostream&, const Wave_half<VTB_REAL_TYPE>&);

template std::ostream&
vtb::core::operator<<(std::ostream&, const Wave_statistics<VTB_REAL_TYPE>&);