/*
* Ball Collision, simulates collision between several
* freely moving rigid balls in 3D (intermediate).
* Uses FreeGLUT for graphics.
* Press SPACE or ENTER to generate new circles, ESC or Q to quit.
*/
#include <cstdlib> // for rand/srand
#include "GL/freeglut.h" // normally should be <>
#include <ctime> // for clock
#include <iostream> // logging to console
#include <cmath> // for sqrt
//
// simple definitions
/// 3D vector
struct Vec3
{
double x, y, z;
explicit Vec3(double nx = 0.0, double ny = 0.0, double nz = 0.0)
: x(nx), y(ny), z(nz)
{
}
Vec3& operator+=(const Vec3 &other)
{
x += other.x;
y += other.y;
z += other.z;
return *this;
}
Vec3& operator-=(const Vec3 &other)
{
x -= other.x;
y -= other.y;
z -= other.z;
return *this;
}
Vec3& operator*=(const double coef)
{
x *= coef;
y *= coef;
z *= coef;
return *this;
}
double length() const
{
return std::sqrt(x * x + y * y + z * z);
}
};
inline Vec3 operator+(const Vec3 &a, const Vec3 &b)
{
return Vec3(a.x + b.x, a.y + b.y, a.z + b.z);
}
inline Vec3 operator-(const Vec3 &a, const Vec3 &b)
{
return Vec3(a.x - b.x, a.y - b.y, a.z - b.z);
}
inline Vec3 operator*(double coef, const Vec3 &v)
{
return Vec3(coef * v.x, coef * v.y, coef * v.z);
}
inline double operator*(const Vec3 &a, const Vec3 &b)
{
return a.x * b.x + a.y * b.y + a.z * b.z;
}
std::ostream& operator<<(std::ostream &os, const Vec3 &v)
{
return os << '(' << v.x << ", " << v.y << ", " << v.z << ')';
}
/// RGBA color
struct Color
{
float red, green, blue, alpha;
explicit Color(float r = 0.0f, float g = 0.0f, float b = 0.0f, float a = 1.0f)
: red(r), green(g), blue(b), alpha(a)
{
}
void select() const
{
glColor4f(red, green, blue, alpha);
}
};
/// rigid ball
struct Ball
{
Vec3 p; ///< position
Vec3 v; ///< velocity
Vec3 a; ///< acceleration
double m; ///< mass
double r; ///< radius
Color color;
Ball
(
const Vec3 &pos = Vec3(),
const Vec3 &vel = Vec3(),
double mass = 1.0,
double rad = 10.0,
const Color &col = Color(1.0f, 1.0f) // yellow by default
)
: p(pos), v(vel), m(mass), r(rad), color(col), a()
{
}
/// OpenGL code to bring Ball to the frame buffer is here
void paint() const
{
glLoadIdentity();
glTranslated(p.x, p.y, p.z);
glScaled(r, r, r);
color.select();
glCallList(1);
}
/// calculate dynamics
void move(double time_step)
{
const Vec3 ta = time_step * a;
p += time_step * (v + 0.5 * ta);
v += ta;
a = Vec3();
}
};
std::ostream& operator<<(std::ostream &os, const Ball &c)
{
return os << "position: " << c.p
<< "\nvelocity: " << c.v
<< "\nmass: " << c.m
<< "\nradius: " << c.r << std::endl;
}
//
// constants and current (global) state
//
/// amount of bodies
const size_t N = 64;
/// relaxation coefficient (should fall into [0, 1])
const double k = 0.8;
/// central force constant
const double G = 0.5;
/// friction constant
const double F = 0.0;
/// actual bodies
Ball C[N];
/// last frame time moment
std::clock_t t;
//
// aux functions
//
/// returns a random number from [minv, maxv)
inline double random(double minv, double maxv)
{
static const double anti_rand_max = 1.0 / double(RAND_MAX);
return minv + (maxv - minv) * std::rand() * anti_rand_max;
}
/// generates randomly moving circles going to collide somewhere
void newBalls()
{
std::cout << "\nGenerating new circles:\n";
// target collision point
const Vec3 cp(random(-300.0, 300.0), random(-300.0, 300.0), random(-300.0, 300.0));
std::cout << "Collision point: " << cp << "\n";
for (size_t i = 0; i < N; ++i)
{
Ball &b = C[i]; // setup C[i] fields
b.p = Vec3(random(-300.0, 300.0), random(-300.0, 300.0), random(-300.0, 300.0)); // position
b.v = 0.25 * (cp - b.p); // 4 seconds before it may reach the collision point
b.r = random(4.0, 20.0); // radius
b.m = random(0.25, 25.0) * b.r * b.r * b.r; // mass
b.color.red = static_cast<float>(random(0.0, 1.0));
b.color.green = static_cast<float>(random(0.0, 1.0));
b.color.blue = static_cast<float>(random(0.0, 1.0));
std::cout << i << ")\n" << b;
}
}
/// gravity
inline void gravity(Ball &b1, Ball &b2)
{
Vec3 r = b2.p - b1.p;
const double dist = r.length();
const double dist3 = dist * dist * dist;
if (dist < b1.r + b2.r)
{
const double dist4 = dist3 * 0.5 + std::numeric_limits<double>::epsilon();
r *= -G / dist4;
}
else
r *= G / dist3;
b1.a += b2.m * r;
b2.a -= b1.m * r;
}
/// viscous friction force
inline void friction(Ball &b)
{
b.a -= F * b.r * b.r / b.m * b.v;
}
//
// GLUT callbacks
//
/// redraw our scene
void repaint()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
// paint all balls
for (size_t i = 0; i < N; ++i)
C[i].paint();
// finish and present our frame to the user
glFinish();
glutSwapBuffers();
}
/// calculate the next frame
void frameMove()
{
static const double anti_clocks_per_sec = 1.0 / double(CLOCKS_PER_SEC);
// check time
std::clock_t t1 = std::clock();
if(t1 != t)
{
double h = anti_clocks_per_sec * (t1 - t);
t = t1;
// check for balls intersection
# pragma omp parallel for
for (int i = 0; i < N; ++i)
{
Ball &a = C[i];
for (size_t j = i + 1; j < N; ++j)
{
Ball &b = C[j];
const Vec3 dp = a.p - b.p;
const double dplen = dp.length();
const double delta = dplen - (a.r + b.r);
if (delta < 0.0 && (a.v - b.v) * dp < 0.0) // are overlapping and approaching?
{
const Vec3 r = 1.0 / dplen * dp;
// collision detected -- calculate new velocities
const Vec3 u(-r.y, r.x);
// old velocities in (r, u) basis
const double
v1r = a.v * r,
v2r = b.v * r,
v1u = a.v * u,
v2u = b.v * u,
// new velocities
v1n = k / (a.m + b.m) * ((a.m - b.m) * v1r + 2.0 * b.m * v2r),
v2n = k / (a.m + b.m) * (2.0 * a.m * v1r + (b.m - a.m) * v2r);
// calculate resulting velocities in original coordinates
a.v = v1n * r + v1u * u;
b.v = v2n * r + v2u * u;
}
}
}
h *= 0.25;
for (int round = 0; round < 4; ++round)
{
// forces
# pragma omp parallel for
for (int i = 0; i < N; ++i)
{
Ball &b = C[i];
//applyCentralForce(b);
for (size_t j = i + 1; j < N; ++j)
gravity(b, C[j]);
}
// move
# pragma omp parallel for
for (int i = 0; i < N; ++i)
C[i].move(h);
}
// say GLUT that our frame has changed
glutPostRedisplay();
}
}
/// work out a keypress
void keyPressed(unsigned char key, int x, int y)
{
switch (key)
{
case 27: // ESC ASCII code
case 'q':
case 'Q':
std::cout << "Exiting\n";
glutExit(); // quit
case ' ':
case '\r':
case '\n':
newBalls();
default:
; // ignore it, suppress compiler's possible warning
}
}
/// work out window resize to width of w pixels and height of h pixels
void reshape(int w, int h)
{
// setup projection
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
if (w < h)
{
const double aspect = double(w) / double(h);
glOrtho(-500.0 * aspect, 500.0 * aspect, -500.0, 500.0, -500.0, 500.0);
}
else
{
const double aspect = double(h) / double(w);
glOrtho(-500.0, 500.0, -500.0 * aspect, 500.0 * aspect, -500.0, 500.0);
}
glMatrixMode(GL_MODELVIEW);
// setup viewport
glViewport(0, 0, w, h);
}
//
// main
//
int main(int argc, char **argv)
{
// initialization
glutInit(&argc, argv);
std::srand(static_cast<unsigned>(std::time(nullptr)));
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_MULTISAMPLE);
glutInitWindowPosition(100, 100);
glutInitWindowSize(800, 600);
glutCreateWindow("Circle collision simulation");
// setup some OpenGL state
glClearColor(0.0, 0.0, 0.0, 1.0); // default one (black), just to show glClearColor
glEnable(GL_COLOR_MATERIAL);
glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
glEnable(GL_DEPTH_TEST);
glEnable(GL_NORMALIZE);
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
// create a disk image
glNewList(1, GL_COMPILE);
glutSolidSphere(1.0, 16, 8);
glEndList();
// register callbacks
glutDisplayFunc(repaint);
glutKeyboardFunc(keyPressed);
glutReshapeFunc(reshape);
glutIdleFunc(frameMove);
// run
std::cout << "Starting...\n" << "k = " << k << "\n";
newBalls();
t = std::clock();
glutMainLoop();
}