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

#include "sph_system.h"

/**
 *
 * @param width 网格或者说摄像机的宽
 * @param height 网格或者说摄像机的高
 * @param u 液体粘度系数
 * @param m 粒子质量
 * @param K 理想气体方程常数
 * @param dens_0 液体静止密度
 * @param h 光滑核半径
 * @param grid_scale 网格单元大小
 */
SPHSystem::SPHSystem(int width, int height, double u, double m, double K, double dens_0, double h, double grid_scale) {
    this->K = K;
    this->h = h;
    this->u = u;
    this->dens_0 = dens_0;
    this->m = m;
    this->grid_scale = grid_scale;
    this->cache = new char[(width+1) * height];
    for (int i = 0; i < height; i++) {
        cache[width + i * (width + 1)] = '\n';
    }
    cache[width + (height - 1) * (width + 1)] = '\0';

    h2 = h*h;
    h4 = h2*h2;
    h5 = h4*h;
    h8 = h4*h4;
    K_poly6 =   m * 4.0 / (PI * h8);
    K_spiky =   m * (double)10.0 / (PI * h5);
    K_viscosity = m * u * 30 / (PI * h4);

    this->grid = {
            .grid = new std::vector<int> [width * height],
            .h = height,
            .w = width
    };
}


#ifdef DEBUG
int cnt = 0;
#endif

void SPHSystem::update() {
    int size = this->particles.size();
    // 计算密度
    for (int id = 0; id < size; id++) {
        Particle& particle = particles[id];
        particle.dens = calculate_density(particle, id);
    }

    // 计算速度与位移
    for (int id = 0; id < size; id++) {
        Particle& pi = particles[id];
        if (pi.type == NORMAL) {
            double t = 0.006;
            Vec2 a_p = calculate_a_pressure(pi, id);
            Vec2 a_v = calculate_a_viscosity(pi, id);
            Vec2 a = g + a_p + a_v;
            Vec2 v0 = pi.v;
            Vec2 v1 = v0 + a * t;
            Vec2 s = (v0 + v1) * 0.5;

            pi.v = v1;
            pi.pos = pi.pos + s;
        }
    }

    // 刷新 grid
    int g_width = grid.w;
    int g_height = grid.h;
    int g_size = g_width * g_height;
    for (int i = 0; i < g_size; i++) {
        grid.grid[i].clear();
    }

    // 更新 grid
    for (int id = 0; id < size; id++) {
        Particle& pi = particles[id];
        int x = (int) floor(pi.pos.x / grid_scale);
        int y = (int) floor(pi.pos.y / grid_scale);
        if (y >= g_height) {
            pi.v.y = -0.5 * pi.v.y;
            pi.pos.y = (g_height - 1) * grid_scale;
        } else if (y < 0) {
            pi.v.y = -0.5 * pi.v.y;
            pi.pos.y = 0;
        }
        if (x >= g_width) {
            pi.v.x = -0.5 * pi.v.x;
            pi.pos.x = (g_width - 1) * grid_scale;
        } else if(x < 0) {
            pi.v.x = -0.5 * pi.v.x;
            pi.pos.x = 0;
        }
        if (x < g_width && x > 0 && y < g_height && y > 0) {
            int pos = x + y * g_width;
            this->grid.grid[pos].push_back(id);
        }
    }
#ifdef DEBUG
    cnt++;

    for (int i = 0; i < size; i++) {
        Particle& pi = particles[i];
        printf("%d, %d: pos = (%.2lf, %.2lf), v = (%.2lf, %.2lf), dens = %.2lf\n",
                      cnt, i,   pi.pos.x, pi.pos.y,     pi.v.x, pi.v.y, pi.dens);
    }

    getchar();

    for (int i = 0; i <= size; i++) {
        printf("\033[1A");
    }
#endif

    // 渲染一帧
    render();

    std::this_thread::sleep_for(std::chrono::milliseconds(50));
}

void SPHSystem::render() const {
    // 使用 Marching Squares 渲染网格
    flush();
    int width = grid.w;
    int height = grid.h;
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            // a --- b --- o
            // |  a  |  b  |
            // d --- c --- o
            // |  d  |  c  |
            // o --- o --- o
            bool a, b, c, d;
            int a_p = x + y * width;
            a = !grid.grid[a_p].empty();
            b = false;
            c = false;
            d = false;
            if (x + 1 < width) {
                int b_p = a_p + 1;
                b = !grid.grid[b_p].empty();
                if (y + 1 < height) {
                    c = !grid.grid[b_p + width].empty();
                }
            }
            if (y + 1 < height) {
                d = !grid.grid[a_p + width].empty();
            }

            bool na = !a;
            bool nb = !b;
            bool nc = !c;
            bool nd = !d;
            if        (na && nb && nc && nd ||  a &&  b &&  c &&  d) {
                // case 1 && case 9
                if (a) {
                    cache[x + y * (width + 1)] = '#';
                } else {
                    cache[x + y * (width + 1)] = ' ';
                }
            } else if (na && nb && nc &&  d ||  a &&  b &&  c && nd) {
                // case 2 && case 10
                cache[x + y * (width + 1)] = '.';
            } else if ( a && nb && nc && nd || na &&  b &&  c &&  d) {
                // case 3 && case 11
                cache[x + y * (width + 1)] = '\'';
            } else if (na &&  b && nc && nd ||  a && nb &&  c &&  d) {
                // case 4 && case 12
                cache[x + y * (width + 1)] = '`';
            } else if (na && nb &&  c && nd ||  a &&  b && nc &&  d) {
                // case 5 && case 13
                cache[x + y * (width + 1)] = ',';
            } else if ( a && nb && nc &&  d || na &&  b &&  c && nd) {
                // case 6 && case 14
                cache[x + y * (width + 1)] = '|';
            } else if (na && nb &&  c &&  d ||  a && b && nc && nd) {
                // case 7 && case 15
                cache[x + y * (width + 1)] = '_';
            } else if (na &&  b && nc &&  d ||  a && nb &&  c && nd) {
                // case 8 && case 16
                cache[x + y * (width + 1)] = '/';
            }
        }
    }
    puts(cache);
}

double SPHSystem::calculate_density(const Particle &pi, int id) {
     double dens = 0.0;
     for (auto& pj: this->particles) {
         if (pi.pos == pj.pos) {
            dens += h4;
             continue;
         }
         Vec2 v = pi.pos - pj.pos;
         double distance = v.x * v.x + v.y * v.y;
         if (distance < h2) {
             double result = density(pi.pos, pj.pos, h);
             dens += result;
         }
     }
    double result = dens * K_poly6;
    return result;
}

Vec2 SPHSystem::calculate_a_pressure(const Particle &pi, int id) {
    Vec2 a = Vec2 {.x = 0.0, .y = 0.0};
    for (auto& pj: this->particles) {
        if (pi.pos == pj.pos) {
            continue;
        }
        Vec2 v = pi.pos - pj.pos;
        double distance = v.x * v.x + v.y * v.y;
        if (distance < h2) {
            Vec2 result = a_pressure(pi.pos, pj.pos, h, pi.dens, pj.dens, dens_0, K);
            a = a + result;
        }
    }

    return a * K_spiky;
}

Vec2 SPHSystem::calculate_a_viscosity(const Particle &pi, int id) {
    Vec2 a = Vec2 {.x = 0.0, .y = 0.0};

    for (auto& pj: this->particles) {
        if (pi.pos == pj.pos) {
            continue;
        }
        Vec2 v = pi.pos - pj.pos;
        double distance = v.x * v.x + v.y * v.y;
        if (distance < h2) {
            Vec2 result = a_viscosity(pi.pos, pj.pos, pi.v, pj.v, pi.dens, pj.dens, h);
            a = a + result;
        }
    }

    return a * K_viscosity;
}