#include <cmath>
#include <matplot/matplot.h>
#include <vector>

std::vector<double> get_x();
std::vector<double> get_y();
std::vector<double> get_z();

int main() {
    using namespace matplot;

    // Positions
    vector_1d x = get_x();
    vector_1d y = get_y();
    vector_1d z = get_z();

    // Vectors
    vector_1d u, v, w, m;
    for (size_t i = 0; i < x.size(); ++i) {
        double t = std::exp(-std::pow(x[i], 2) - std::pow(y[i], 2) -
                            std::pow(z[i], 2));
        u.emplace_back(x[i] * t);
        v.emplace_back(y[i] * t);
        w.emplace_back(z[i] * t);
        m.emplace_back(std::sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]));
    }

    quiver3(x, y, z, u, v, w, m, 0.3)->normalize(true).line_width(2);

    show();
    return 0;
}

std::vector<double> get_x() {
    return {-1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,
            -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,   -1,
            -1,   -1,   -1,   -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
            -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
            -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, 0,    0,    0,    0,    0,
            0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,
            0,    0,    0,    0,    0,    0,    0,    0,    0,    0.5,  0.5,
            0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,
            0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,
            0.5,  1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
            1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,
            1,    1,    1,    1};
}

std::vector<double> get_y() {
    return {-1,   -1,   -1, -1,  -1,  -0.5, -0.5, -0.5, -0.5, -0.5, 0,    0,
            0,    0,    0,  0.5, 0.5, 0.5,  0.5,  0.5,  1,    1,    1,    1,
            1,    -1,   -1, -1,  -1,  -1,   -0.5, -0.5, -0.5, -0.5, -0.5, 0,
            0,    0,    0,  0,   0.5, 0.5,  0.5,  0.5,  0.5,  1,    1,    1,
            1,    1,    -1, -1,  -1,  -1,   -1,   -0.5, -0.5, -0.5, -0.5, -0.5,
            0,    0,    0,  0,   0,   0.5,  0.5,  0.5,  0.5,  0.5,  1,    1,
            1,    1,    1,  -1,  -1,  -1,   -1,   -1,   -0.5, -0.5, -0.5, -0.5,
            -0.5, 0,    0,  0,   0,   0,    0.5,  0.5,  0.5,  0.5,  0.5,  1,
            1,    1,    1,  1,   -1,  -1,   -1,   -1,   -1,   -0.5, -0.5, -0.5,
            -0.5, -0.5, 0,  0,   0,   0,    0,    0.5,  0.5,  0.5,  0.5,  0.5,
            1,    1,    1,  1,   1};
}

std::vector<double> get_z() {
    return {-1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1, -1, -0.5, 0, 0.5, 1,
            -1, -0.5, 0, 0.5, 1};
}