I've been struggling with this for a while. I dunno if any of you have noticed my n-body sim code anywhere, but I've discovered the collision method I was using is more bullcrap than I thought it was (I knew it was innacurate but not that innacurate), so I've been struggling to find a decent collision algorithm that will preserve momentum and rotational momentum over many steps. I have the physics pretty good, so that the total momentum of the system only changes in the range of 10^-5 or 10^-6 per step... but the rotational velocity of bodies keeps increasing
Meaning, if you have, say, 20 bodies colliding, at first they behave like you would expect, and form an approximate circle... but then - seemingly out of nowhere - the 20 bodies start spinning in one direction, faster and faster until the centripetal force rips the bodies apart, all from nowhere.
Here's the source code:
#include "math.h"
#include <vector>
#include "DarkGDK.h"
#include "D3DFunc.h"
using namespace std;
class circPhy;
class vector2
{
public:
float x;
float y;
vector2();
vector2(const float& a,const float& b);
const vector2 operator+ (const vector2& b) const; //vector addition
vector2& operator+= (const vector2& b);
const vector2 operator- (const vector2& b) const; //vector subtraction
vector2& operator-= (const vector2& b);
const vector2 operator* (const float& b) const; //vector scalar multiplication
vector2& operator*= (const float& b);
const vector2 operator/ (const float& b) const; //vector scalar division
vector2& operator/= (const float& b);
const float operator* (const vector2& b) const; //dot product
void normalize();
void normalize(const float& dist);
vector2 normalized();
vector2 normalized(const float& dist);
float length();
};
static const vector2 operator* (const float& b, const vector2& a);
class verletPoint
{
public:
vector2 pos;
vector2 vel;
vector2 acc;
vector2 lAcc;
verletPoint();
verletPoint(const vector2& position,const vector2& velocity,const vector2& acceleration);
verletPoint(const vector2& position,const vector2& velocity);
verletPoint(const vector2& position);
void integrate(const float& timestep);
};
class circPhy
{
public:
static float G;
static float COR;
float M;
float R;
bool updateable;
verletPoint center;
circPhy(){}
circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable);
void update(const float& timestep);
void coll(circPhy* colObj);
void addGrav(const vector2& pos, const float& mass);
void zeroAcc();
};
void DarkGDK(void)
{
dbSetDisplayMode(640,480,32);
int number_of_bodies = 1;
int lastclick=0;
d3dInit();
dbSyncOn();
dbRandomize(1);
vector<circPhy> r;
int moveobj=0;
for(int n=0;n<number_of_bodies;n++)
r.push_back(circPhy(dbRnd(dbScreenWidth()),dbRnd(dbScreenHeight()),20, 10,true));
while ( LoopGDK ( ) )
{
dbCLS();
dbInk(dbRGB (0, 0, 255),dbRGB (0, 0, 255));
vector2 sum=vector2(0.0,0.0);
if(!dbSpaceKey())
{
for(int n=0;n<number_of_bodies;n++)
{
for(int i=0;i<number_of_bodies;i++)
if(i!=n)
r[n].addGrav(r[i].center.pos,r[i].M);
}
for(int n=0;n<number_of_bodies;n++)
{
for(int i=0;i<number_of_bodies;i++)
if(i!=n)
r[n].coll(&r[i]);
}
for(int n=0;n<number_of_bodies;n++)
{
sum+=r[n].center.vel*r[n].M;
r[n].update(2.0f);
}
}
if(dbMouseClick() && lastclick!=dbMouseClick())
{
if(dbMouseClick()==1)
{
r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(20)+10,dbRnd(10)+2,true));
number_of_bodies++;
}
if(dbMouseClick()==2)
{
r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,1,false));
number_of_bodies++;
}
if(dbMouseClick()==4)
{
r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(100)+5,100,false));
number_of_bodies++;
moveobj=r.size()-1;
}
}
if(dbMouseClick()==4)
{
r[moveobj].center.pos=vector2(dbMouseX(),dbMouseY());
}
if(lastclick==4 && dbMouseClick()!=4)
{
r.erase(r.begin()+moveobj);
number_of_bodies--;
}
lastclick=dbMouseClick();
for(int n=0;n<number_of_bodies;n++)
d3dCircle(r[n].center.pos.x,r[n].center.pos.y,r[n].R,0);
dbText(0,0,"number of bodies:");
dbText(0,16,dbStr(number_of_bodies));
dbText(0,32,"framerate (capped @ 60):");
dbText(0,48,dbStr(dbScreenFPS()));
dbText(0,64,dbStr((float)sum.x));
dbText(0,80,dbStr((float)sum.y));
dbSync();
}
return;
}
//class vector2
vector2::vector2()
{
x = y = 0.0f;
}
vector2::vector2(const float& a,const float& b)
{
x=a;
y=b;
}
const vector2 vector2::operator+ (const vector2& b) const
{
vector2 temp;
temp.x=x+b.x;
temp.y=y+b.y;
return temp;
}
vector2& vector2::operator+= (const vector2& b)
{
if(this!=&b)
{
x=x+b.x;
y=y+b.y;
}
return *this;
}
const vector2 vector2::operator-(const vector2& b) const
{
vector2 temp;
temp.x=x-b.x;
temp.y=y-b.y;
return temp;
}
vector2& vector2::operator-=(const vector2& b)
{
if(this!=&b)
{
x=x-b.x;
y=y-b.y;
}
return *this;
}
const float vector2::operator* (const vector2& b) const
{
return (x*b.x+y*b.y);
}
const vector2 vector2::operator* (const float& b) const
{
vector2 temp;
temp.x=x*b;
temp.y=y*b;
return temp;
}
vector2& vector2::operator*= (const float& b)
{
x=x*b;
y=y*b;
return *this;
}
const vector2 vector2::operator/ (const float& b) const
{
vector2 temp=vector2(0.0f,0.0f);
if(b>0)
{
temp.x=x/b;
temp.y=y/b;
}
return temp;
}
vector2& vector2::operator/= (const float& b)
{
if(b>0)
{
x=x/b;
y=y/b;
}
return *this;
}
void vector2::normalize()
{
float dist = sqrt(x*x+y*y);
x/=dist;
y/=dist;
}
void vector2::normalize(const float& dist)
{
x/=dist;
y/=dist;
}
vector2 vector2::normalized()
{
float dist = sqrt(x*x+y*y);
return vector2(x/dist,y/dist);
}
vector2 vector2::normalized(const float& dist)
{
return vector2(x/dist,y/dist);
}
float vector2::length()
{
float n=x*x+y*y;
return sqrt(n>0?n:0);
}
static const vector2 operator* (const float& b,const vector2& a)
{
return a*b;
}
verletPoint::verletPoint()
{
vector2 pos (0.0f,0.0f);
vector2 vel (0.0f,0.0f);
vector2 acc (0.0f,0.0f);
vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position, const vector2 & velocity, const vector2 & acceleration)
{
pos=position;
vel=velocity;
acc=acceleration;
vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position,const vector2 & velocity)
{
pos=position;
vel=velocity;
vector2 acc (0.0f,0.0f);
vector2 lAcc(0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position)
{
pos=position;
vector2 vel (0.0f,0.0f);
vector2 acc (0.0f,0.0f);
vector2 lAcc(0.0f,0.0f);
}
void verletPoint::integrate(const float& timestep)
{
// vel+=acc*timestep;
//pos+=vel*timestep;
pos+=vel*timestep+lAcc*timestep*timestep/2; //P(t+dt)~=P(t)+V(t)*dt+A(t)*((dt)^2)/2
vel+=(lAcc+acc)*timestep/2; //V(t+dt)~=V(t)+(A(t)+A(t+dt))*dt/2
lAcc=acc;
}
float circPhy::G=2.0f;
float circPhy::COR=0.0f;
circPhy::circPhy(const float& center_x,const float& center_y,const float& radius,const float& mass, const bool& moveable)
{
M=mass;
R=radius;
updateable=moveable;
center=verletPoint(vector2(center_x,center_y));
}
void circPhy::update(const float& timestep)
{
if(updateable)
center.integrate(timestep);
zeroAcc();
}
void circPhy::coll(circPhy* colObj)
{
float d=0.0f;
vector2 x=center.pos-colObj->center.pos;
d=x.length();
x.normalize(d);
float r1=R+colObj->R;
if(r1>d)
{
float v1x=center.vel*x;
float v2x=colObj->center.vel*x;
center.vel=center.vel-(v1x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
colObj->center.vel=colObj->center.vel-(v2x-((M*v1x+colObj->M*v2x)/(M+colObj->M)))*x;
center.pos+=(r1-d)/2*x;
colObj->center.pos-=(r1-d)/2*x;
}
}
void circPhy::addGrav(const vector2& pos, const float& mass)
{
vector2 d=pos-center.pos;
float n=d.length();
center.acc+=d*(G*mass/(n*n*n));
}
void circPhy::zeroAcc()
{
center.acc=vector2(0.0f,0.0f);
}
If you run it (d3d_circle is used, but you could replace that with whatever command you want to draw a dot, circle, or box; the effect is the same), you can see that the total momentum (x and y values are the second to last and last values printed on the top left of the screen) stays close to zero. If you have two bodies collide (click to create one), you can see them behave as expected, even with preserving momentum when they collide. But... as you add more and more, you can see that it all spins out of control!!!
What am i doing wrong? First I find the force on every point, then I update collision on every point, then I update velocity and move every point. I've tried switching the order of those to no avail, and the problem doesn't just occur with objects of differing radius/mass/density, though it happens faster with smaller/denser objects.
I could imagine that when two objects aren't moving/accelerating towards each other in exactly opposite directions, they would slide and this would add to rotational velocity... but i've no clue how to solve this.
Is't life, I ask, is't even prudence, to bore thyself and bore thy students?
