diff --git a/demo/Table.h b/demo/Table.h index cb54441..e87b1fc 100644 --- a/demo/Table.h +++ b/demo/Table.h @@ -6,13 +6,14 @@ #include #include #include +#include #include namespace phy { struct Particle { - /// @brief Kinematic properties (position: y0, velocity: y1). - dyn::Yoshida kin; + /// @brief Kinematic properties. + std::complex xy, v; /// @brief Mass and radius. float mass = 1.0f, radius = 1.0f; @@ -28,11 +29,11 @@ struct Particle { /// @param r Radius, positive. constexpr Particle(std::complex xy, std::complex v, float m = 1, float r = 1) - : kin{xy, v}, mass{m}, radius{r} {} + : xy{xy}, v{v}, mass{m}, radius{r} {} /// @brief Construct a circle that represents this particle. [[nodiscard]] constexpr dyn::Circle circle() const { - return {kin.y0, radius}; + return {xy, radius}; } }; @@ -68,7 +69,10 @@ class Table : public std::vector { }; // step(): Integrate. // (Calls accel() above a few times with slightly different xy.) - copy[i].kin.step(dt, accel); + auto step = dyn::Yoshida<>{p.xy, p.v}; + step.step(dt, accel); + auto &q = copy[i]; + q.xy = step.y0, q.v = step.y1; } *this = copy; } @@ -79,10 +83,25 @@ class Table : public std::vector { dyn::Kahan> cm; dyn::Kahan m; for (auto &&p : *this) - cm += double(p.mass) * C<>(p.kin.y0), m += double(p.mass); + cm += double(p.mass) * C<>(p.xy), m += double(p.mass); auto c = C(cm() / m()); for (auto &&p : *this) - p.kin.y0 -= c; + p.xy -= c; + } + + /// @brief Refresh the "disk" used for parts of the calculation. + void refresh_disk() noexcept { gr.refresh_disk(); } + + /// @brief Test whether the simulation is in "good state." + bool good() noexcept { + auto constexpr goodf = [](float f) { return std::isfinite(f); }; + auto constexpr goodc = [](std::complex f) { + return std::isfinite(f.real()) && std::isfinite(f.imag()); + }; + for (auto &&p : *this) + if (!goodc(p.xy) || !goodc(p.v)) + return false; + return true; } }; diff --git a/demo/main.cpp b/demo/main.cpp index ab72dfa..c1968fd 100644 --- a/demo/main.cpp +++ b/demo/main.cpp @@ -1,12 +1,12 @@ #include #include -#include #include #include #include #include +#include #include #include "Table.h" @@ -19,7 +19,7 @@ static int do_main() { auto constexpr RADIUS = 0.05f; auto constexpr MASS = 1.0f; auto constexpr too_far = [](Particle &p) { - return std::abs(p.kin.y0) > 5'000.0f; + return std::abs(p.xy) > 5'000.0f; }; // Quasi-random number generator. @@ -183,7 +183,7 @@ static int do_main() { auto o = std::complex(n.x, n.y); auto p = random_particle(); // Set location - p.kin.y0 = o; + p.xy = o; table.push_back(p); interactive.spawned_last_frame = true; @@ -195,19 +195,18 @@ static int do_main() { } simulate: - // Clear floating point exceptions. - std::feclearexcept(FE_ALL_EXCEPT); - // Do the simulation! if (!interactive.freeze) table.step(dt); + table.refresh_disk(); + // Don't be moving the particles while it's being manipulated. - if (!interactive.spawned_last_frame && !interactive.other_user_manip) - table.center(); + // if (!interactive.spawned_last_frame && !interactive.other_user_manip) + // table.center(); // Inspect floating point exceptions. - if (std::fetestexcept(FE_DIVBYZERO | FE_INVALID)) + if (!table.good()) // NaN or infinity somewhere. Reset the simulation. goto reset_sim; @@ -216,16 +215,21 @@ static int do_main() { // Draw all particles in the frame. BeginMode2D(cam); - for (auto &&p : table) { - auto cp = p.circle(); - Vector2 xy{cp.real(), cp.imag()}; - // Should I issue the draw call or not? Tell it here: - // co stores screen dimensions calculated earlier. - auto co = cam.offset; - // x, y, width, then height. - Rectangle screen{-co.x / 2, -co.y / 2, co.x, co.y}; - if (CheckCollisionCircleRec(xy, cp.radius, screen)) - DrawCircleV(xy, cp.radius, WHITE); + { + // Compute the screen's less-less (ll) and greater-greater (gg) corners in + // the world coordinate system so that draw calls would only be issued for + // certainly visible particles. + std::complex scwh{float(GetScreenWidth()), + float(GetScreenHeight())}; + std::complex scof{cam.offset.x, cam.offset.y}; + auto scsc = (scwh - scof) / cam.zoom; + std::complex sctg{cam.target.x, cam.target.y}; + auto ll = sctg - scof / cam.zoom, gg = sctg + scsc; + for (auto &&p : table) { + auto cp = p.circle(); + if (dyn::disk_arrect_isct(cp, ll, gg)) + DrawCircleV({cp.real(), cp.imag()}, cp.radius, WHITE); + } } EndMode2D(); diff --git a/dyn/CMakeLists.txt b/dyn/CMakeLists.txt index 8efc6e3..5882d92 100644 --- a/dyn/CMakeLists.txt +++ b/dyn/CMakeLists.txt @@ -5,5 +5,6 @@ add_library(dyn INTERFACE kahan.h halton.h newton.h - circle.h) + circle.h + verlet.h) target_include_directories(dyn INTERFACE .) diff --git a/dyn/circle.h b/dyn/circle.h index 132ceb4..5791ff6 100644 --- a/dyn/circle.h +++ b/dyn/circle.h @@ -5,6 +5,7 @@ /// @brief Represent a circle by the center and radius. #include +#include namespace dyn { @@ -22,6 +23,55 @@ template struct Circle : public std::complex { : std::complex(center), radius(radius) {} }; +/// @brief Decide whether at least one intersection (point) exists between the +/// area of the disk centered at the origin with the given radius and the area +/// of the given rectangle (degenerate cases are unspecified). +/// @param ll The less-less corner of the rectangle. +/// @param gg The greater-greater corner of the rectangle. +template +bool origin_disk_arrect_isct(F radius, std::complex ll, + std::complex gg) noexcept { + // sin(45 deg) [also cos(45 deg)]. + constexpr auto sc45 = std::numbers::sqrt2_v / F(2); + + // Decide whether the given rectangle and a square centered at (0,0) having a + // "radius" (half side length) of `radius` (thus a full side length of twice + // the `radius`). + auto arrec_arsq_isct = [gg, ll](auto radius) { + return ll.real() < radius && -radius < gg.real() && ll.imag() < radius && + -radius < gg.imag(); + }; + + // The rectangle does not touch the bounding square of the circle. + if (!arrec_arsq_isct(radius)) + return false; + + // The rectangle touches the interior square of the circle. + if (arrec_arsq_isct(radius * sc45)) + return true; + + // At least one of the four corners of the given rectangle is in the disk. + F cc[4] = {std::abs(ll), std::abs(gg), + std::abs(std::complex{ll.real(), gg.imag()}), + std::abs(std::complex{gg.real(), ll.imag()})}; + for (auto &&c : cc) + if (c < radius) + return true; + + // They aren't touching. + return false; +} + +/// @brief Decide whether at least one intersection (point) exists between a +/// circular disk and the area of a rectangle (degenerate cases are +/// unspecified). +template +bool disk_arrect_isct(Circle circ, std::complex ll, + std::complex gg) noexcept { + ll -= circ, gg -= circ; + return origin_disk_arrect_isct(circ.radius, ll, gg); +} + } // namespace dyn #endif // GRASS_CIRCLE_H diff --git a/dyn/newton.h b/dyn/newton.h index e579c3a..3a3dd43 100644 --- a/dyn/newton.h +++ b/dyn/newton.h @@ -23,6 +23,14 @@ template class Gravity { /// for Monte Carlo integration in the case of overlapping circles. std::array, N_MONTE> disk{}; + // Each Halton sequence (a kind of low-discrepancy sequence) creates an + // evenly spaced set of points on the unit interval (0, 1); unlike the + // uniform distribution, however, the points look "uniformly distributed" + // (number of points being mostly proportional to length of any subset) + // even for a finite sample of points. + Halton h2; + Halton h3; + public: /// @brief Create an instance with a quasi-random internal state. constexpr Gravity() { refresh_disk(); } @@ -54,14 +62,6 @@ template class Gravity { /// case of intersecting circles) with new evenly distributed points on the /// unit disk centered about the origin. Call time to time. constexpr void refresh_disk() { - // Each Halton sequence (a kind of low-discrepancy sequence) creates an - // evenly spaced set of points on the unit interval (0, 1); unlike the - // uniform distribution, however, the points look "uniformly distributed" - // (number of points being mostly proportional to length of any subset) - // even for a finite sample of points. - Halton h2; - Halton h3; - // Fill `disk` with random points on unit disk centered about origin // by rejection sampling. for (auto &&p : disk) { diff --git a/dyn/verlet.h b/dyn/verlet.h new file mode 100644 index 0000000..be7b1cc --- /dev/null +++ b/dyn/verlet.h @@ -0,0 +1,21 @@ +#ifndef GRASS_VERLET_H +#define GRASS_VERLET_H + +#include + +namespace dyn { + +template struct Verlet { + std::complex y0, y1; + + template void step(F h, A y2) { + auto a = y2(y0); + y0 += h * y1 + h * h * F(0.5) * a; + auto b = y2(y0); + y1 += h * F(0.5) * (a + b); + } +}; + +} // namespace dyn + +#endif // GRASS_VERLET_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index cb27347..a8f2b94 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -11,7 +11,8 @@ include(GoogleTest) # Now simply link against gtest or gtest_main as needed. Eg add_executable(units yoshida_test.cpp - newton_test.cpp) + newton_test.cpp + circle_test.cpp) target_precompile_headers(units INTERFACE "gtest/gtest.h") target_link_libraries(units gtest_main dyn) gtest_discover_tests(units) diff --git a/tests/CircleTest0.png b/tests/CircleTest0.png new file mode 100644 index 0000000..9490d77 Binary files /dev/null and b/tests/CircleTest0.png differ diff --git a/tests/CircleTest1.png b/tests/CircleTest1.png new file mode 100644 index 0000000..d98c362 Binary files /dev/null and b/tests/CircleTest1.png differ diff --git a/tests/CircleTest2.png b/tests/CircleTest2.png new file mode 100644 index 0000000..9dbc353 Binary files /dev/null and b/tests/CircleTest2.png differ diff --git a/tests/CircleTest3.png b/tests/CircleTest3.png new file mode 100644 index 0000000..d29be3d Binary files /dev/null and b/tests/CircleTest3.png differ diff --git a/tests/circle_test.cpp b/tests/circle_test.cpp new file mode 100644 index 0000000..03fea44 --- /dev/null +++ b/tests/circle_test.cpp @@ -0,0 +1,58 @@ +#include "gtest/gtest.h" + +#include + +#include + +/// @brief See `CircleTest0.png` for reference. +class CircleTest0 : public testing::Test { +protected: + // xy, r; rectangles (2 corners) + dyn::Circle const circle{0, 2.4f}; + std::complex const less_less{-4.0f, -4.0f}; + std::complex const greater_greater{-2.0f, -2.0f}; +}; + +TEST_F(CircleTest0, Disjoint0) { + ASSERT_FALSE(dyn::disk_arrect_isct(circle, less_less, greater_greater)); +} + +class CircleTest1 : public testing::Test { +protected: + dyn::Circle const circle{0, 2.0f}; + std::complex const less_less{-5.0f, 1.0f}, + greater_greater{-1.0f, 5.0f}; +}; + +TEST_F(CircleTest1, In0) { + ASSERT_TRUE(dyn::disk_arrect_isct(circle, less_less, greater_greater)); +} + +class CircleTest2 : public testing::Test { +protected: + dyn::Circle const circle{0, 3.0f}; + std::complex const less_less{0.6f, 2.7f}, greater_greater{1.8f, 3.6f}; +}; + +TEST_F(CircleTest2, In0) { + ASSERT_TRUE(dyn::disk_arrect_isct(circle, less_less, greater_greater)); +} + +class CircleTest3 : public testing::Test { +protected: + dyn::Circle const circle{0, 4.0f}; + std::complex const c{-2.0f, -2.0f}, e{0, -1.0f}, l{-2.0f, -5.0f}, + n{5.0f, 3.0f}, s{7.0f, 2.0f}, u{9.0f, 3.0f}; +}; + +TEST_F(CircleTest3, In0) { + ASSERT_TRUE(dyn::disk_arrect_isct(circle, c, e)); +} + +TEST_F(CircleTest3, In1) { + ASSERT_TRUE(dyn::disk_arrect_isct(circle, l, n)); +} + +TEST_F(CircleTest3, Out0) { + ASSERT_FALSE(dyn::disk_arrect_isct(circle, s, u)); +}