diff --git a/examples/gravity.cpp b/examples/gravity.cpp index abf0e1df..f6e7fa9f 100644 --- a/examples/gravity.cpp +++ b/examples/gravity.cpp @@ -6,6 +6,11 @@ #include #include #include +#include +#include +#include +#include +#include #include @@ -47,44 +52,197 @@ struct particle geom::point pos; geom::vector vel; + float angle; + float angle_vel; + float radius; float mass; + float density; + + geom::vector delta_pos{0.f, 0.f}; + geom::vector delta_vel{0.f, 0.f}; + + geom::vector acc{0.f, 0.f}; float T = 0.f; }; +float const G = 50.f; +float const GG = 0*1000.f; +float const GC = 0*100000.f; + +float const dt = 0.01f; + +float const world_size = 10000.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_; +}; + struct myapp : app::app { myapp() - : app("Test app") + : app("Test app", 4) , rng_{std::random_device{}()} { gl::ClearColor(1.f, 1.f, 1.f, 1.f); + vsync(true); - std::uniform_real_distribution d{-200.f, 200.f}; + 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}; - for (int i = 0; i < 1000; ++i) + float min_R = 0.f; + float max_R = 50.f; + std::uniform_real_distribution rR{0.f, 1.f}; + + bool star = false; + + 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}); + +// float planet_R[] = {200.f, 300.f, 400.f}; +// float planet_a[] = {0.f, geom::rad(120.f), geom::rad(240.f)}; + +// if(false) + for (int i = 0; i < 500; ++i) { - float r = rr(rng_); - float m = r * r * r; - - float v = 0.1f; + geom::vector v{0.f, 0.f}; geom::point p; + float m; + float r; + float den; while (true) { - p = {((i % 2) ? -300.f : 300.f) + d(rng_), d(rng_)}; +// p = {((i % 2) ? -200.f : 200.f) + d(rng_), d(rng_)}; + + r = rr(rng_); + den = rden(rng_); + m = geom::pi * r * r * den; + + auto a = ra(rng_); + float R = std::sqrt(rR(rng_)) * (max_R - min_R) + min_R; +// if (std::uniform_int_distribution(0, 1)(rng_) == 0) +// R += 200.f; + +// auto pl = std::uniform_int_distribution(0, std::size(planet_R) - 1)(rng_); +// a = planet_a[pl]; +// R = planet_R[pl]; + +// R += std::uniform_real_distribution(-10.f, 10.f)(rng_); +// a += geom::rad(std::uniform_real_distribution(-2.f, 2.f)(rng_)); + + p = + { + R * std::cos(a), + R * std::sin(a), + }; + p += geom::vector{((i % 2) ? -1000.f : 1000.f), 0.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; } - particles_.push_back({p, {v * d(rng_), v * d(rng_)}, r, m}); -// particles_.push_back({{ 400.f + d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); + particles_.push_back({p, v, 0.f, 0.f, r, m, den}); + } -// particles_.push_back({{d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); -// particles_.push_back({{d(rng), d(rng)}, {0.f, 0.f}, 1.f, 1.f}); + float total_M = 0.f; + if (star) + total_M += particles_[0].mass; + for (std::size_t i = star ? 1 : 0; i < particles_.size(); ++i) + total_M += particles_[i].mass / 2.f; + + if(false) + for (std::size_t i = star ? 1 : 0; i < particles_.size(); ++i) + { + auto r = particles_[i].pos - geom::point{0.f, 0.f}; + auto R = geom::length(r); + float V = std::sqrt(G * total_M / R); + particles_[i].vel = geom::ort(r / R) * V; + (void)V; } } @@ -112,15 +270,27 @@ struct myapp : app::app geom::point const screen_center { window_size_[0] / 2, window_size_[1] / 2 }; auto target = camera_center_ + flip_y(geom::cast(*mouse_ - screen_center) * scale); + force_target_ = target; - geom::point pos{100.f, 0.f}; +// geom::point pos{100.f, 0.f}; - geom::vector vel = geom::normalized(target - pos) * 500.f; +// geom::vector vel = geom::normalized(target - pos) * 40.f; - particles_.push_back({pos, vel, 1.f, 1.f}); +// particles_.push_back({pos, vel, 1.f, 1.f}); + +// for (auto & p : particles_) +// { +// auto r = p.pos - target; +// p.vel += 4000.f * r / geom::length_sqr(r) / p.mass; +// } } } + void on_left_button_up() override + { + force_target_ = std::nullopt; + } + void on_right_button_down() override { camera_drag_ = mouse_; @@ -144,48 +314,147 @@ struct myapp : app::app camera_center_ += flip_y(geom::cast(*camera_drag_ - *mouse_) * scale); camera_drag_ = mouse_; } + + if (force_target_) + { + geom::scale const flip_y({1.f, -1.f}); + float const scale = camera_size_ / window_size_[1]; + geom::point const screen_center { window_size_[0] / 2, window_size_[1] / 2 }; + auto target = camera_center_ + flip_y(geom::cast(*mouse_ - screen_center) * scale); + force_target_ = target; + } + } + + void on_key_down(SDL_Keycode key) override + { + app::app::on_key_down(key); + + 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}); + } } void update() override { -// std::shuffle(particles_.begin(), particles_.end(), rng_); - - float const dt = 0.01f; - float const G = 100.f; - for (std::size_t step = 0; step < 1; ++step) { { util::clock<> clock; + for (auto & p : particles_) + p.acc = {0.f, 0.f}; + + for (auto & p : particles_) + p.acc += geom::vector{0.f, -GG}; + + for (auto & p : particles_) + { + auto r = (p.pos - p.pos.zero()); + p.acc -= GC * r / std::pow(1.f + geom::length(r), 3.f); + } + +// for (std::size_t i = 0; i < particles_.size(); ++i) +// { +// log::info() << "Start: #" << i << " pos = " << std::setprecision(10) << particles_[i].pos << ", vel = " << particles_[i].vel +// << ", |vel| = " << geom::length(particles_[i].vel); +// } + 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 = std::max(particles_[i].radius + particles_[j].radius, length(r)); float const l = length(r); auto const f = G * particles_[i].mass * particles_[j].mass * r / std::pow(l, 3.f); - particles_[i].vel -= f * dt / particles_[i].mass; - particles_[j].vel += f * dt / particles_[j].mass; +// log::info() << "Force: #" << i << "," << j << " = " << std::setprecision(10) << f; + + particles_[i].acc -= f / particles_[i].mass; + particles_[j].acc += f / particles_[j].mass; } } total_forces_ += clock.count(); } +// for (std::size_t i = 0; i < particles_.size(); ++i) +// { +// log::info() << "Force: #" << i << " = " << std::setprecision(10) << particles_[i].acc; +// } + + if (force_target_ && false) + { + for (auto & p : particles_) + { + auto r = p.pos - *force_target_; + p.acc += 100000.f * r / geom::length_sqr(r) / p.mass; + } + } + +// for (auto & p : particles_) +// { +// auto old = p.pos; +// p.pos += (p.pos - p.old_pos) + p.acc * dt * dt; +// p.old_pos = old; +// } + +// for (std::size_t i = 0; i < particles_.size(); ++i) +// { +// log::info() << "Integrate: #" << i << " new pos = " << std::setprecision(10) << particles_[i].pos; +// } + + if (false) + { + float s = 100.f; + world_center[0] += dt * std::uniform_real_distribution(-s, s)(rng_); + world_center[1] += dt * std::uniform_real_distribution(-s, s)(rng_); + } + + if (force_target_) + world_center = *force_target_; + for (auto & p : particles_) { + p.vel += p.acc * dt; p.pos += p.vel * dt; + p.angle += p.angle_vel * dt; } { util::clock<> clock; - for (std::size_t iteration = 0; iteration < 10; ++iteration) + for (std::size_t iteration = 0; iteration < 1; ++iteration) { + for (auto & p : particles_) + { + auto r = p.pos - world_center; + auto l = geom::length(r); + + if (l + p.radius > world_size) + { + auto n = r / l; + + auto t = geom::ort(n); + + auto vn = n * geom::dot(p.vel, n); + + auto vt = t * geom::dot(p.vel + t * p.angle_vel * p.radius, t); + + auto pt = - vt * (1.f - std::exp(- 10.f * dt)); + float I = 0.5f * p.mass * geom::sqr(p.radius); + + p.pos -= n * (l + p.radius - world_size); + p.vel -= 1.75f * vn; + + p.vel += pt / p.mass; + p.angle_vel += geom::det(n * p.radius, pt) / I; + } + } + for (std::size_t i = 0; i < particles_.size(); ++i) { for (std::size_t j = i + 1; j < particles_.size(); ++j) @@ -225,24 +494,111 @@ struct myapp : app::app auto const vij = particles_[i].vel - particles_[j].vel; - if (l < R && dot(vij, n) < 0.f) + // merging collision + if (l < R && false) { - auto const d = n * (R - l) / 2.f * 0.5f;// * (particles_[i].mass + particles_[j].mass); - particles_[i].pos += d; - particles_[j].pos -= d; + particle p; - auto dp = - n * dot(n, vij) * 2.f / (1.f / particles_[i].mass + 1.f / particles_[j].mass); + auto M = (particles_[i].mass + particles_[j].mass); + p.pos = particles_[i].pos + (particles_[j].pos - particles_[i].pos) * particles_[j].mass / M; + p.vel = (particles_[i].vel * particles_[i].mass + particles_[j].vel * particles_[j].mass) / M; + p.mass = M; + p.radius = std::sqrt(geom::sqr(particles_[i].radius) + geom::sqr(particles_[j].radius)); + p.density = p.mass / (geom::pi * geom::sqr(p.radius)); - particles_[i].vel += dp / particles_[i].mass; - particles_[j].vel -= dp / particles_[j].mass; + particles_[i] = p; + particles_.erase(particles_.begin() + j); + break; + } + + // inelastic collision + // if (l < R && dot(vij, n) < 0.f) + if (l < R && false) + { + float D = R - l; + + float A = 0.5f * (1.f / particles_[i].mass + 1.f / particles_[j].mass); + float B = dot(n, vij); + auto np = n * (- B / A); +// np *= 0.5f + 0.5f * std::exp(- 0.f * dt); + + auto p = particles_[j].pos + (particles_[j].radius - D / 2.f) * n; + + auto ri = p - particles_[i].pos; + auto rj = p - particles_[j].pos; + + auto t = geom::ort(n); + + auto vri = geom::dot(t, particles_[i].vel + geom::ort(ri) * particles_[i].angle_vel); + auto vrj = geom::dot(t, particles_[j].vel + geom::ort(rj) * particles_[j].angle_vel); + + auto vrij = vri - vrj; + + auto tp = - t * vrij * (1.f - std::exp(- 10.f * dt));; + + particles_[i].vel += (np + tp) / particles_[i].mass; + particles_[j].vel -= (np + tp) / particles_[j].mass; + + float Ii = 0.5f * particles_[i].mass * geom::sqr(particles_[i].radius); + float Ij = 0.5f * particles_[j].mass * geom::sqr(particles_[j].radius); + + auto ti = geom::det(ri, tp); + auto tj = geom::det(rj, -tp); + + particles_[i].angle_vel += ti / Ii; + particles_[j].angle_vel += tj / Ij; + + 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) + { + auto f = 10000.f * n * (R - l); + particles_[i].vel += dt * f / particles_[i].mass; + particles_[j].vel -= dt * f / particles_[j].mass; + + auto vn = geom::dot(vij, n) * n; +// vn *= 1.f - std::exp(- (1.f - l / R)); + 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) + { + float D = (R - l) * 1.f; + float ki = D * particles_[j].mass / (particles_[i].mass + particles_[j].mass); + float kj = - D * particles_[i].mass / (particles_[i].mass + particles_[j].mass); + + auto di = ki * n; + auto dj = kj * n; + + particles_[i].pos += di; + particles_[j].pos += dj; } } } + + for (auto & p : particles_) + { + p.vel += p.delta_vel; + p.pos += p.delta_pos; + p.delta_pos = {0.f, 0.f}; + p.delta_vel = {0.f, 0.f}; + } } total_collisions_ += clock.count(); } + float Ep = 0.f; float Ek = 0.f; @@ -256,20 +612,41 @@ struct myapp : app::app } } - log::info() << "Energy: " << Ek << " - " << (-Ep) << " = " << (Ek + Ep); + float omega = 0.f; + for (std::size_t i = 0; i < particles_.size(); ++i) + { + omega += geom::det(particles_[i].mass * particles_[i].vel, particles_[i].pos - geom::point{0.f, 0.f}); + } + +// log::info() << "Angular velocity: " << omega; + +// log::info() << "Energy: " << Ek << " - " << (-Ep) << " = " << (Ek + Ep); ++frame_count_; + +// std::reverse(particles_.begin(), particles_.end()); } } void present() override { + gl::ClearColor(0.f, 0.f, 0.1f, 0.f); gl::Clear(gl::COLOR_BUFFER_BIT); + painter_.circle(world_center, world_size, {255, 255, 255, 255}, 72); + for (auto & p : particles_) { - float c = 2.f / (std::exp(-p.T / 100000.f) + 1.f) - 1.f; - painter_.circle(p.pos, p.radius, {static_cast(c * 255.f), 0, 0, 255}); +// float c = 2.f / (std::exp(-p.T / 100000.f) + 1.f) - 1.f; + float c = 1.f - p.density; + auto x = static_cast(c * 255.f); + + 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 * 0.75f, {255, 255, 255, 255}); } geom::orthographic_camera camera; @@ -287,6 +664,8 @@ struct myapp : app::app { log::info() << "Avg forces time: " << (total_forces_ / frame_count_); log::info() << "Avg collision time: " << (total_collisions_ / frame_count_); + + prof::dump(); } private: @@ -294,6 +673,8 @@ private: geom::vector window_size_; + std::optional> force_target_; + geom::point camera_center_ { 0.f, 0.f }; float camera_size_ = 500.f; float camera_ratio_ = 1.f; @@ -308,6 +689,8 @@ private: int frame_count_ = 0; float total_forces_ = 0.f; float total_collisions_ = 0.f; + + audio::engine audio_; }; int main()