From b0bd606727e0f011df75489b65f2120a9159517b Mon Sep 17 00:00:00 2001 From: lisyarus Date: Thu, 3 Nov 2022 22:15:40 +0300 Subject: [PATCH] More gravity experiments --- examples/gravity.cpp | 768 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 675 insertions(+), 93 deletions(-) diff --git a/examples/gravity.cpp b/examples/gravity.cpp index f6e7fa9f..e47f6165 100644 --- a/examples/gravity.cpp +++ b/examples/gravity.cpp @@ -6,11 +6,13 @@ #include #include #include +#include #include #include #include #include #include +#include #include @@ -62,6 +64,8 @@ struct particle geom::vector delta_pos{0.f, 0.f}; geom::vector delta_vel{0.f, 0.f}; + geom::point old_pos{0.f, 0.f}; + geom::vector acc{0.f, 0.f}; float T = 0.f; @@ -70,88 +74,16 @@ struct particle float const G = 50.f; float const GG = 0*1000.f; float const GC = 0*100000.f; +float const K = 10000.f; +float const FR = 0.99f; float const dt = 0.01f; -float const world_size = 10000.f; +float const world_size = 10000000.f; geom::point world_center{0.f, 0.f}; -struct sound_stream - : audio::stream -{ - sound_stream(std::vector const & particles) - : particles_(particles) - {} - - std::optional length() const override - { - return std::nullopt; - } - - std::size_t read(float * data, std::size_t sample_count) override - { - std::size_t time = played_.load() / 2; - - if (velocity_.size() < particles_.size()) - { - velocity_.resize(particles_.size(), 0.f); - amplitude_.resize(particles_.size(), 0.f); - oscillator_.resize(particles_.size(), 0.f); - } - - float mean_v = 0.f; - for (std::size_t i = 0; i < particles_.size(); ++i) - { - float v = geom::length(particles_[i].vel); - velocity_[i] = geom::lerp(velocity_[i], v, 0.1f); - mean_v += velocity_[i]; - } - mean_v /= particles_.size(); - - for (std::size_t i = 0; i < particles_.size(); ++i) - { -// oscillator_[i].frequency(500.f * velocity_[i] / mean_v); -// oscillator_[i].frequency(geom::lerp(100.f, 10000.f, velocity_[i] / 1000.f)); - oscillator_[i].frequency(10.f * velocity_[i]); -// oscillator_[i].frequency(0.1f * geom::sqr(velocity_[i])); - } - - - for (std::size_t i = 0; i < sample_count; i += 2) - { - float v = 0.f; - for (std::size_t j = 0; j < particles_.size(); ++j) - v += oscillator_[j].next().imag(); - - if (i >= 2) - { - v = (v + data[i - 1]) / 2.f; - } - - data[i + 0] = v; - data[i + 1] = v; - - time += 1; - } - - played_.fetch_add(sample_count); - - return sample_count; - } - - // The number of samples already played from this stream - std::size_t played() const override - { - return played_.load(); - } - -private: - std::vector const & particles_; - std::vector velocity_; - std::vector amplitude_; - std::vector oscillator_; - std::atomic played_; -}; +int const SOLVE_ITERATIONS = 16; +float const BIAS = 1.f / 8.f; struct myapp : app::app { @@ -163,15 +95,13 @@ struct myapp : app::app vsync(true); - audio_.output()->stream(audio::compressor(std::make_shared(particles_), audio::from_db(-4.f), 0.95f)); - std::uniform_real_distribution d{-50.f, 50.f}; std::uniform_real_distribution rr{0.5f, 2.f}; std::uniform_real_distribution rden{0.25f, 1.f}; std::uniform_real_distribution ra{0.f, 2.f * geom::pi}; float min_R = 0.f; - float max_R = 50.f; + float max_R = 100.f; std::uniform_real_distribution rR{0.f, 1.f}; bool star = false; @@ -179,8 +109,10 @@ struct myapp : app::app if (star) particles_.push_back({{0.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 10.f, geom::pi * 10000.f, 1.f}); -// particles_.push_back({{-1.f, 0.f}, {0.f, -1.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); -// particles_.push_back({{ 1.f, 0.f}, {0.f, 1.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); +// particles_.push_back({{-2.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); +// particles_.push_back({{ 2.f, 0.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); +// particles_.push_back({{ 0.f, 3.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); +// particles_.push_back({{ 0.f, -3.f}, {0.f, 0.f}, 0.f, 0.f, 1.f, 1.f, 1.f}); // float planet_R[] = {200.f, 300.f, 400.f}; // float planet_a[] = {0.f, geom::rad(120.f), geom::rad(240.f)}; @@ -219,9 +151,9 @@ struct myapp : app::app R * std::cos(a), R * std::sin(a), }; - p += geom::vector{((i % 2) ? -1000.f : 1000.f), 0.f}; +// p += geom::vector{((i % 2) ? -1000.f : 1000.f), 0.f}; - v[0] = {(i % 2) ? 100.f : -100.f}; +// v[0] = {(i % 2) ? 100.f : -100.f}; if (std::all_of(particles_.begin(), particles_.end(), [&](particle const & q){ return geom::distance(q.pos, p) > q.radius + r; })) break; @@ -244,10 +176,15 @@ struct myapp : app::app particles_[i].vel = geom::ort(r / R) * V; (void)V; } + + for (auto & p : particles_) + p.old_pos = p.pos; } void on_resize(int width, int height) override { + app::on_resize(width, height); + gl::Viewport(0, 0, width, height); window_size_ = {width, height}; @@ -331,12 +268,21 @@ struct myapp : app::app if (key == SDLK_SPACE) { - particles_.push_back({{200.f, 20.f}, {-1000.f, 0.f}, 0.f, 0.f, 1.f, 100.f, 1.f}); + paused_ = !paused_; + } + + if (key == SDLK_c) + { + particles_.push_back({{200.f, 0.f}, {-5000.f, 0.f}, 0.f, 0.f, 1.f, 100.f, 1.f}); + particles_.back().old_pos = particles_.back().pos; } } void update() override { + if (paused_) + return; + for (std::size_t step = 0; step < 1; ++step) { { @@ -378,6 +324,39 @@ struct myapp : app::app } } + // Force-based collisions +// if (false) + for (std::size_t i = 0; i < particles_.size(); ++i) + { + for (std::size_t j = i + 1; j < particles_.size(); ++j) + { + auto const r = particles_[i].pos - particles_[j].pos; + float const l = length(r); + + float R = particles_[i].radius + particles_[j].radius; + auto n = r / l; + + if (l < R) + { + auto f = K * n * (R - l); + (void)n; + +// auto const f = -2.f * G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f); +// f -= 2.f * G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f); + + particles_[i].acc += f / particles_[i].mass; + particles_[j].acc -= f / particles_[j].mass; + } + +// auto vij = particles_[i].vel - particles_[j].vel; + +// auto vn = geom::dot(vij, n) * n; + +// particles_[i].acc -= vn * particles_[j].mass / (particles_[i].mass + particles_[j].mass); +// particles_[j].acc += vn * particles_[i].mass / (particles_[i].mass + particles_[j].mass); + } + } + total_forces_ += clock.count(); } @@ -386,6 +365,12 @@ struct myapp : app::app // log::info() << "Force: #" << i << " = " << std::setprecision(10) << particles_[i].acc; // } + if (is_key_down(SDLK_m)) + { + for (auto & p : particles_) + p.vel *= std::exp(-100.f*dt); + } + if (force_target_ && false) { for (auto & p : particles_) @@ -417,13 +402,466 @@ struct myapp : app::app if (force_target_) world_center = *force_target_; +// if (false) for (auto & p : particles_) { p.vel += p.acc * dt; + p.vel *= FR; p.pos += p.vel * dt; - p.angle += p.angle_vel * dt; } + // Verlet + if (false) + for (auto & p : particles_) + { + auto old = p.pos; + p.pos += (p.pos - p.old_pos) + p.acc * dt * dt; + p.old_pos = old; + } + + // New iterative algorithm + if (false) + { + std::vector> old_pos(particles_.size()); + for (std::size_t i = 0; i < particles_.size(); ++i) + old_pos[i] = particles_[i].pos; + + struct collision + { + std::size_t i, j; + geom::vector n; // i -> j + }; + + std::vector collisions; + + for (std::size_t i = 0; i < particles_.size(); ++i) + { + for (std::size_t j = i + 1; j < particles_.size(); ++j) + { + auto const r = particles_[i].pos - particles_[j].pos; + + float const R = particles_[i].radius + particles_[j].radius; + float const l = length(r); + + if (l < R) + collisions.push_back({i, j, - r / l}); + } + } + + for (std::size_t iteration = 0; iteration < 16; ++iteration) + { + for (auto const & c : collisions) + { + auto const r = particles_[c.j].pos - particles_[c.i].pos; + auto const R = particles_[c.j].radius + particles_[c.i].radius; + auto const l = geom::dot(r, c.n); + + if (l < R) + { + auto const M = particles_[c.j].mass + particles_[c.i].mass; + + auto const d = (l - R) * c.n / M; + + particles_[c.i].pos += d * particles_[c.j].mass; + particles_[c.j].pos -= d * particles_[c.i].mass; + } + } + } + + for (std::size_t i = 0; i < particles_.size(); ++i) + particles_[i].vel += (particles_[i].pos - old_pos[i]) / dt; + } + + // First-event algorithm + if (false) + { +// std::vector particle_time(particles_.size(), 0.f); + + for (std::size_t i = 0; i < particles_.size(); ++i) + { + bool collided = false; + + for (std::size_t j = 0; j < particles_.size(); ++j) + { + if (i == j) continue; + +// auto dp = (particles_[i].pos - particles_[i].vel * particle_time[i]) - (particles_[j].pos - particles_[j].vel * particle_time[j]); + auto dp = particles_[i].pos - particles_[j].pos; + auto dv = particles_[i].vel - particles_[j].vel; + auto R = particles_[i].radius + particles_[j].radius; + + float A = geom::dot(dv, dv); + float B = geom::dot(dp, dv); + float C = geom::dot(dp, dp) - geom::sqr(R); + + if (B >= 0.f) return; + + float collision_time; + + if (C <= 0.f) + { + collision_time = 0.f; + } + else + { + float D = B * B - A * C; + + if (D <= 0.f) + continue; + + collision_time = C / (- B + std::sqrt(D)); + } + + if (collision_time >= dt) + continue; + + collided = true; + particles_[i].pos += particles_[i].vel * collision_time; + +// particles_[i].pos += particles_[i].vel * (e.t - particle_time[i]); +// particles_[j].pos += particles_[j].vel * (e.t - particle_time[j]); + +// particle_time[i] = e.t; +// particle_time[j] = e.t; + + auto n = geom::normalized(particles_[i].pos - particles_[j].pos); + + float S = 0.5f * (1.f / particles_[i].mass + 1.f / particles_[j].mass); + float T = dot(n, particles_[i].vel - particles_[j].vel); + + auto np = n * (- T / S); + + np *= std::exp(-10.f * dt); + + particles_[i].vel += np / particles_[i].mass; + particles_[j].vel -= np / particles_[j].mass; + + break; + } + + if (!collided) + particles_[i].pos += particles_[i].vel * dt; + } + } + + // Event-based algorithm + if (false) + { +// log::info() << "===== ITERATION ====="; + collisions_ = 0; + struct collision_event + { + std::uint32_t i0, i1; + float t; + bool erased = false; + }; + + struct comparator + { + bool operator()(collision_event const & e1, collision_event const & e2) const + { + return std::tie(e1.t, e1.i0, e1.i1) < std::tie(e2.t, e2.i0, e2.i1); + } + }; + + using events_container = std::set; + events_container events; + std::vector erased_events; + std::vector> event_list(particles_.size()); + + std::vector particle_time(particles_.size(), 0.f); + + std::unordered_map, std::vector> cells; + + float time = 0.f; + + auto check_collision = [&](std::size_t i, std::size_t j) + { + if (i > j) + std::swap(i, j); + + // (pi + vi * (t - ti) - pj - vj * (t - tj))^2 = (ri+rj)^2 + + auto dp = (particles_[i].pos - particles_[i].vel * particle_time[i]) - (particles_[j].pos - particles_[j].vel * particle_time[j]); + auto dv = particles_[i].vel - particles_[j].vel; + auto R = particles_[i].radius + particles_[j].radius; + + float A = geom::dot(dv, dv); + float B = geom::dot(dp, dv); + float C = geom::dot(dp, dp) - geom::sqr(R); + +// log::info() << "DV = " << dv; +// log::info() << "DP = " << dp; +// log::info() << "B = " << B; + +// if (C < 0.f) +// log::info() << "DEFINITELY A COLLISION!"; + +// if (B > 0.f) return; + +// auto result = geom::solve_quadratic(A, B, C); +// if (!result) return; + +// log::info() << result->first << " " << result->second; + +// float min_time = 0.f; +// min_time = std::max(particle_time[i], particle_time[j]); + +// std::optional collision_time; + +// if (result->first >= min_time && result->first <= dt) +// collision_time = result->first; +// else if (result->second >= min_time && result->second <= dt) +// collision_time = result->second; + +// if (result->first >= 0.f && result->first <= dt) +// if (result->first <= dt) +// collision_time = result->first; +// collision_time = std::max(min_time, result->first); +// (void)min_time; + +// if (!collision_time) return; + + if (B >= -0.01f * geom::length(dv) * geom::length(dp)) return; + + float collision_time; + + if (C <= 0.f) + { + collision_time = time; + } + else + { + float D = B * B - A * C; + + if (D <= 0.f) + return; + + collision_time = C / (- B + std::sqrt(D)); + } + + if (collision_time >= dt) return; + + auto dpn = (particles_[i].pos + particles_[i].vel * (collision_time - particle_time[i])) - (particles_[j].pos + particles_[j].vel * (collision_time - particle_time[j])); + if (geom::dot(dv, dpn) >= -0.01f * geom::length(dv) * geom::length(dpn)) return; + + auto insert_result = events.insert(collision_event{i, j, collision_time}); + if (insert_result.second) + { + event_list[i].push_back(insert_result.first); + event_list[j].push_back(insert_result.first); + } + }; + + auto erase = [&](events_container::iterator it) + { + if (!it->erased) + { + auto n = events.extract(it); + n.value().erased = true; + erased_events.push_back(std::move(n)); + } + }; + + auto clear_collisions = [&](std::size_t i) + { + for (auto it : event_list[i]) + erase(it); + event_list[i].clear(); + }; + + auto foreach_cells = [&](std::size_t i, auto && callback) + { + geom::box bbox; + + bbox |= particles_[i].pos - geom::vector{1.f, 1.f} * particles_[i].radius; + bbox |= particles_[i].pos + geom::vector{1.f, 1.f} * particles_[i].radius; + + auto npos = particles_[i].pos + particles_[i].vel * dt;//(dt - particle_time[i]); + + bbox |= npos - geom::vector{1.f, 1.f} * particles_[i].radius; + bbox |= npos + geom::vector{1.f, 1.f} * particles_[i].radius; + + int xmin = std::floor(bbox[0].min); + int xmax = std::floor(bbox[0].max); + int ymin = std::floor(bbox[1].min); + int ymax = std::floor(bbox[1].max); + + for (int x = xmin; x <= xmax; ++x) + { + for (int y = ymin; y <= ymax; ++y) + { + callback(cells[{x, y}]); + } + } + }; + + auto clear_cells = [&](std::size_t i) + { + foreach_cells(i, [i](auto & cell){ + auto it = std::find(cell.begin(), cell.end(), i); + if (it != cell.end()) + cell.erase(it); + }); + }; + + auto fill_cells = [&](std::size_t i) + { + foreach_cells(i, [i](auto & cell){ cell.push_back(i); }); + }; + + auto find_collisions = [&](std::size_t i) + { + foreach_cells(i, [&, i](auto & cell){ + for (auto j : cell) + if (i != j) check_collision(i, j); + }); + }; + + auto handle_collision = [&](collision_event const & e) + { + clear_cells(e.i0); + clear_cells(e.i1); + + particles_[e.i0].pos += particles_[e.i0].vel * (e.t - particle_time[e.i0]); + particles_[e.i1].pos += particles_[e.i1].vel * (e.t - particle_time[e.i1]); + +// log::info() << "Collision " << e.i0 << " " << e.i1 << " " << geom::distance(particles_[e.i0].pos, particles_[e.i1].pos); + + particle_time[e.i0] = e.t; + particle_time[e.i1] = e.t; + + auto n = geom::normalized(particles_[e.i0].pos - particles_[e.i1].pos); + + float A = 0.5f * (1.f / particles_[e.i0].mass + 1.f / particles_[e.i1].mass); + float B = dot(n, particles_[e.i0].vel - particles_[e.i1].vel); + + auto BB = B / geom::length(particles_[e.i0].vel - particles_[e.i1].vel); + (void)BB; + auto np = n * (- B / A); + + np *= std::exp(-10.f * dt); + + particles_[e.i0].vel += np / particles_[e.i0].mass; + particles_[e.i1].vel -= np / particles_[e.i1].mass; + + fill_cells(e.i0); + fill_cells(e.i1); + }; + + // 0 - naive + // 1 - optimized + int algorithm = 1; + + if (algorithm == 0) + { + collisions_ = 0; + + while (time < dt) + { + particle_time.assign(particles_.size(), 0.f); + + collision_event next{0, 0, std::numeric_limits::infinity()}; + + for (std::size_t i = 0; i < particles_.size(); ++i) + { + for (std::size_t j = i + 1; j < particles_.size(); ++j) + { + auto rij = particles_[i].pos - particles_[j].pos; + auto vij = particles_[i].vel - particles_[j].vel; + auto R = particles_[i].radius + particles_[j].radius; + + float B = geom::dot(rij, vij); + + if (B >= 0.f) continue; + + std::optional e; + + float C = geom::dot(rij, rij) - geom::sqr(R); + + if (C <= 0.f) + { + e = collision_event{i, j, 0.f}; + } + else + { + float A = geom::dot(vij, vij); + + float D = B * B - A * C; + + if (D <= 0.f) + continue; + + float t = C / (- B + std::sqrt(D)); + e = collision_event{i, j, t}; + } + + if (e && e->t < next.t) + next = *e; + } + } + + if (time + next.t < dt) + { + handle_collision(next); + ++collisions_; + + for (std::size_t i = 0; i < particles_.size(); ++i) + particles_[i].pos += particles_[i].vel * (next.t - particle_time[i]); + + time += next.t; + } + else + break; + } + + for (std::size_t i = 0; i < particles_.size(); ++i) + particles_[i].pos += particles_[i].vel * (dt - time); + } + else if (algorithm == 1) + { + for (std::size_t i = 0; i < particles_.size(); ++i) + fill_cells(i); + + for (std::size_t i = 0; i < particles_.size(); ++i) + find_collisions(i); + + while (!events.empty()) + { + collision_event e = *events.begin(); + time = e.t; + +// if (e.t < time) +// { +// erase(events.begin()); +// continue; +// } + +// time = e.t; + +// if (e.t < particle_time[e.i0]) +// log::info() << e.t << " < " << particle_time[e.i0]; +// if (e.t < particle_time[e.i1]) +// log::info() << e.t << " < " << particle_time[e.i1]; + + handle_collision(e); + clear_collisions(e.i0); + clear_collisions(e.i1); + + find_collisions(e.i0); + find_collisions(e.i1); + + ++collisions_; + } + + for (std::size_t i = 0; i < particles_.size(); ++i) + { + particles_[i].pos += particles_[i].vel * (dt - particle_time[i]); + } + } + } + + if (false) { util::clock<> clock; @@ -512,7 +950,7 @@ struct myapp : app::app } // inelastic collision - // if (l < R && dot(vij, n) < 0.f) +// if (l < R && dot(vij, n) < 0.f) if (l < R && false) { float D = R - l; @@ -555,8 +993,30 @@ struct myapp : app::app particles_[j].pos += kj * n; } + // inelastic collision without rotation + if (l < R && dot(vij, n) < 0.f && false) +// if (l < R) + { + float D = R - l; + + float elasticity = 0.25f; + + float A = (1.f / particles_[i].mass + 1.f / particles_[j].mass); + float B = dot(n, vij); + auto np = n * (- B / A) * (1.f + elasticity); + + particles_[i].vel += np / particles_[i].mass; + particles_[j].vel -= np / particles_[j].mass; + + float ki = D * particles_[j].mass / (particles_[i].mass + particles_[j].mass); + float kj = - D * particles_[i].mass / (particles_[i].mass + particles_[j].mass); + + particles_[i].pos += ki * n; + particles_[j].pos += kj * n; + } + // collision replaced by force - if (l < R) + if (l < R && false) { auto f = 10000.f * n * (R - l); particles_[i].vel += dt * f / particles_[i].mass; @@ -564,14 +1024,14 @@ struct myapp : app::app auto vn = geom::dot(vij, n) * n; // vn *= 1.f - std::exp(- (1.f - l / R)); - vn *= 0.5f; +// vn *= 0.5f; particles_[i].vel -= vn * particles_[j].mass / (particles_[i].mass + particles_[j].mass); particles_[j].vel += vn * particles_[i].mass / (particles_[i].mass + particles_[j].mass); } // Verlet collision - if (l < R && false) + if (l < R) { float D = (R - l) * 1.f; float ki = D * particles_[j].mass / (particles_[i].mass + particles_[j].mass); @@ -583,6 +1043,19 @@ struct myapp : app::app particles_[i].pos += di; particles_[j].pos += dj; } + + // impulse-based collision + if (l < R && dot(n, particles_[i].vel - particles_[j].vel) < 0.f && false) + { + float bias = 0.5f; + float constraint = l - R; + float reduced_mass = 1.f / (1.f / particles_[i].mass + 1.f / particles_[j].mass); + float lambda = (- geom::dot(n, particles_[i].vel - particles_[j].vel) - bias / dt * constraint) * reduced_mass; + auto impulse = lambda * n; + + particles_[i].vel += impulse / particles_[i].mass; + particles_[j].vel -= impulse / particles_[j].mass; + } } } @@ -598,6 +1071,70 @@ struct myapp : app::app total_collisions_ += clock.count(); } + // Verlet finalization + if(false) + for (auto & p : particles_) + p.vel = (p.pos - p.old_pos) / dt; + + // Iterative impulse-based collision + if (false) + { + struct collision + { + std::size_t i, j; + geom::vector normal; + float value; + float impulse = 0.f; + }; + + std::vector collisions; + + for (std::size_t i = 0; i < particles_.size(); ++i) + { + for (std::size_t j = i + 1; j < particles_.size(); ++j) + { + auto const r = particles_[i].pos - particles_[j].pos; + + float const l = length(r); + + float const R = particles_[i].radius + particles_[j].radius; + + auto const n = r / l; + + if (l < R) + { + collisions.push_back({i, j, n, l - R}); + } + } + } + + for (std::size_t iteration = 0; iteration < SOLVE_ITERATIONS; ++iteration) + { + for (auto & c : collisions) + { + float bias = BIAS; + + float reduced_mass = 1.f / (1.f / particles_[c.i].mass + 1.f / particles_[c.j].mass); + + auto dv = geom::dot(c.normal, particles_[c.i].vel - particles_[c.j].vel); + float elasticity = 0.75f; + + float lambda = (- dv - bias / dt * c.value - dv * elasticity) * reduced_mass; + float new_impulse = std::max(0.f, c.impulse + lambda); + auto delta = new_impulse - c.impulse; + c.impulse = new_impulse; + + particles_[c.i].vel += delta * c.normal / particles_[c.i].mass; + particles_[c.j].vel -= delta * c.normal / particles_[c.j].mass; + } + } + } + + if (false) + for (auto & p : particles_) + { + p.pos += p.vel * dt; + } float Ep = 0.f; float Ek = 0.f; @@ -618,6 +1155,8 @@ struct myapp : app::app omega += geom::det(particles_[i].mass * particles_[i].vel, particles_[i].pos - geom::point{0.f, 0.f}); } + energy_ = Ek + Ep; + // log::info() << "Angular velocity: " << omega; // log::info() << "Energy: " << Ek << " - " << (-Ep) << " = " << (Ek + Ep); @@ -635,6 +1174,14 @@ struct myapp : app::app painter_.circle(world_center, world_size, {255, 255, 255, 255}, 72); +// for (auto const & c : cells_) +// { +// float a = 0.25f; +// if (!c.second.empty()) +// a = 0.5f; +// painter_.rect({{{c.first[0], c.first[0] + 1.f}, {c.first[1], c.first[1] + 1.f}}}, gfx::to_coloru8(gfx::color_4f{1.f, 0.f, 0.f, a})); +// } + for (auto & p : particles_) { // float c = 2.f / (std::exp(-p.T / 100000.f) + 1.f) - 1.f; @@ -644,8 +1191,8 @@ struct myapp : app::app float s = window_size_[1] / camera_size_; float r = std::max(p.radius * s, 1.5f) / s; - painter_.circle(p.pos, r, {x, x, x, 255}); - painter_.line(p.pos, p.pos + r * geom::direction(p.angle), r / 16.f, {255, 0, 0, 255}, false); + painter_.circle(p.pos, r, {x, x, x, 191}); +// painter_.line(p.pos, p.pos + r * geom::direction(p.angle), r / 16.f, {255, 0, 0, 255}, false); // painter_.circle(p.pos, r * 0.75f, {255, 255, 255, 255}); } @@ -658,6 +1205,35 @@ struct myapp : app::app camera.box[2].max = 1.f; painter_.render(camera.transform()); + + float rotation = 0.f; + for (auto const & p : particles_) + rotation += geom::det(p.pos - p.pos.zero(), p.vel) * p.mass; + rotation_.push(rotation); + + { + gfx::painter::text_options opts; + opts.scale = 3.f; + opts.c = {0, 0, 0, 255}; + opts.x = gfx::painter::x_align::center; + opts.y = gfx::painter::y_align::top; + + auto put = [&](geom::point const & pos, std::string const & str) + { + return; + painter_.text(pos, str, opts); + painter_.text(pos + geom::vector{1.f, 0.f}, str, opts); + }; + + put({width() / 2.f, 30.f}, util::to_string("ITERATIONS: ", SOLVE_ITERATIONS)); + put({width() / 2.f, 60.f}, BIAS == 1.f ? "BIAS: 1" : util::to_string("BIAS: 1/", std::round(1.f / BIAS))); + +// painter_.text({10.f, 10.f}, util::to_string("Angular velocity: ", rotation_.average()), opts); +// painter_.text({10.f, 30.f}, util::to_string("Collisions: ", collisions_, " = ", (collisions_ * 1.f / particles_.size() / particles_.size()), " N^2"), opts); +// painter_.text({10.f, 50.f}, util::to_string("Energy: ", energy_), opts); + } + + painter_.render(geom::window_camera{width(), height()}.transform()); } ~myapp() @@ -676,7 +1252,7 @@ private: std::optional> force_target_; geom::point camera_center_ { 0.f, 0.f }; - float camera_size_ = 500.f; + float camera_size_ = 150.f; float camera_ratio_ = 1.f; std::optional> mouse_; @@ -691,6 +1267,12 @@ private: float total_collisions_ = 0.f; audio::engine audio_; + + util::moving_average rotation_{300}; + int collisions_ = 0; + float energy_ = 0.f; + + bool paused_ = true; }; int main()