void PhysicsEngine::applyForce(Body *b1, Body *b2)
{
    double EPS = 10;      // softening parameter (just to avoid infinities)
    glm::dvec3 p1 = b1->getPos();
    glm::dvec3 p2 = b2->getPos();

    double dist = glm::distance(p1, p2);

    double fmag = (G * b1->getMass() * b2->getMass()) / (dist * dist + EPS * EPS);

    glm::dvec3 f1 = b1->getForce();
    glm::dvec3 f2 = b2->getForce();

    glm::dvec3 f = (p2 - p1) / dist;
    f *= fmag;

    f1 += f;
    f2 -= f;

    b1->updateForce(f1);
    b2->updateForce(f2);
}