//
// Created by HP on 2021/12/26.
//

#include "sph_math.h"

double Vec2::module() const {
    double t_x = this->x;
    double t_y = this->y;
    double xs = t_x * t_x;
    double ys = t_y * t_y;
    double m = std::sqrt(xs + ys);
    return m;
}

double Vec2::operator*(const Vec2 &other) const {
    double xs = this->x * other.x;
    double ys = this->y * other.y;
    return xs + ys;
}

Vec2 Vec2::operator*(double scale) const {
    double new_x = this->x * scale;
    double new_y = this->y * scale;
    Vec2 new_v = Vec2 {.x = new_x, .y = new_y};
    return new_v;
}

Vec2 Vec2::operator-(const Vec2 &other) const {
    double new_x = this->x - other.x;
    double new_y = this->y - other.y;
    Vec2 new_v = Vec2 {.x = new_x, .y = new_y};
    return new_v;
}

Vec2 Vec2::operator+(const Vec2 &other) const {
    double new_x = this->x + other.x;
    double new_y = this->y + other.y;
    Vec2 new_v = Vec2 {.x = new_x, .y = new_y};
    return new_v;
}

bool Vec2::operator==(const Vec2 &other) const {
    return this->x == other.x && this->y == other.y;
}

double density(const Vec2& i, const Vec2& j, double h) {
    // (h^2 - |i - j|)^3
    Vec2 new_v = i - j;

    double vxs = new_v.x * new_v.x;
    double vys = new_v.y * new_v.y;
    double vms = vxs + vys;

    double hs = h * h;

    double difference = hs - vms;

    double result = difference * difference * difference;

    return result;
}

Vec2 a_pressure(const Vec2& i, const Vec2& j, double h,
                double dens_i, double dens_j, double dens_0, double K) {
    // 计算单个粒子产生的压力
    double p_i = K * (dens_i - dens_0);
    double p_j = K * (dens_j - dens_0);
    double p_avg = (p_i + p_j) / 2.0;

    // 计算 r = |i - j|
    Vec2 new_v = i - j;
    double r = new_v.module();
    // 将 new_v 标准化
    new_v = new_v * (1.0 / r);

    // 计算结果
    double hrs = (h - r) * (h - r);
    double scale = (p_avg * hrs) / (dens_i * dens_j);
    Vec2   result = new_v * scale;
    return result;
}

Vec2 a_viscosity(const Vec2& i, const Vec2& j, const Vec2& vi, const Vec2& vj,
                 double dens_i, double dens_j, double h) {
    Vec2 new_v = vj - vi;

    // 计算 r = |i - j|
    Vec2 new_vec = i - j;
    double r = new_vec.module();

    // 计算结果
    double difference = h - r;
    double scale = difference / (dens_i * dens_j);
    Vec2 result = new_v * scale;
    return result;
}