[GDK] O(n^2) n-body gravitational simulation

Author

Message

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 21st Mar 2010 04:15 Edited at: 21st Mar 2010 04:46

Link

Basically, I'm working on a large-scale capable n-body gravitational simulation. What I've worked up is more a cool bit of collision code along with a demonstration of verlet integration, wrapped up with gravitational forces.

What is inefficient about this is that I have a for loop that is basically:

for a = 1 to n
for b = 1 to n
- dostuff
next
next

I'm going to streamline this process and make it faster by using a barnes-hutt algorithm (which implements a quadtree). Here's the source code:

(this is my second project in c++, and my first time using a lot of features of c++. There might be errors in something-or-other that I did, but it's only ~350 lines.)

#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);
		vector2 operator+	(const vector2& b);	//vector addition
		vector2& operator+=	(const vector2& b);
		vector2 operator-	(const vector2& b);	//vector subtraction
		vector2& operator-=	(const vector2& b);
		vector2 operator*	(const float& b);	//vector scalar multiplication
		vector2& operator*=	(const float& b);
		vector2 operator/	(const float& b);	//vector scalar division
		vector2& operator/=	(const float& b);
		float operator*		(const vector2& b);	//dot product
		void normalize();
		void normalize(const float& dist);
		vector2 normalized();
		vector2 normalized(const float& dist);
		float length();
};
class verletPoint
{   
    public:   
        vector2 pos;
        vector2 lastpos;
        vector2 accel;
        verletPoint();
        verletPoint(const vector2& position,const vector2& last_position,const vector2& acceleration);
        verletPoint(const vector2& position,const vector2& last_position);
        verletPoint(const vector2& position);
        void updateVerlet(const float timestep);
};
struct collisionDataTemp
{
    float depth;
    vector2 normal;
    circPhy* bod;
};
class circPhy
{
    public:
        float M;
        float R;
		bool updateable;
        verletPoint center;
        collisionDataTemp coll;
        circPhy(){}
        circPhy(const float center_x,const float center_y,const float radius, bool moveable);
        void updateVerlet(const float timestep);
        bool checkColl(circPhy* colObj);
        bool checkColl(circPhy* colObj,const float dist);
        void processColl();
};

void DarkGDK(void)
{
	dbSetDisplayMode(640,480,32);
	int number_of_bodies = 5;
	int lastclick=0;
	d3dInit();
	dbSyncOn();
	dbSyncRate(60);
	dbRandomize(dbTimer());
	vector<circPhy> r;
	int moveobj=0;
	for(int n=0;n<number_of_bodies;n++)
		r.push_back(circPhy(dbRnd(dbScreenWidth()),dbRnd(dbScreenHeight()),dbRnd(5)+5,true));
	while ( LoopGDK ( ) )
	{
		dbCLS();
		dbInk(dbRGB (0, 0, 255),dbRGB (0, 0, 255));
		if(!dbSpaceKey())
			for(int n=0;n<number_of_bodies;n++)
			{
				if(dbControlKey()==1)
				{
					r[n].center.accel=vector2();
					r[n].center.lastpos=r[n].center.pos;
				}
				else
				{
					vector2 accel=vector2(0.0f,0.0f);
					for(int i=0;i<number_of_bodies;i++)
					{
						if(i!=n && r[i].updateable)
						{
							vector2 temp = r[i].center.pos-r[n].center.pos;
							float dist = temp.length();
							temp.normalize(dist);
							//now dist is the distance between the circle and the mouse, and accel is a unit vector
							dist=r[n].M*r[i].M*.001f/(dist*dist);
							//dist=(dist>.1f)?.1f:dist;
							temp*=dist;
							accel+=temp;
						}
					}
					r[n].center.accel=accel;
				}
				float* t=&r[n].center.pos.x;
				*t=(*t>=0)?  ((*t<=dbScreenWidth())?*t:dbScreenWidth())  :0;
				t=&r[n].center.pos.y;
				*t=(*t>=0)?  ((*t<=dbScreenHeight())?*t:dbScreenHeight())  :0;
				for(int i=0;i<number_of_bodies;i++)
				{
					if(r[n].checkColl(&r[i]))
					{r[n].processColl();}
				}
				r[n].updateVerlet(2);
			}

if(dbMouseClick() && lastclick!=dbMouseClick())
		{
			if(dbMouseClick()==1)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,true));
				number_of_bodies++;
			}
			if(dbMouseClick()==2)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,false));
				number_of_bodies++;
			}
			if(dbMouseClick()==4)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(100)+5,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()));
		dbCircle(365,351,63);
		dbSync();
	}
	return;
}
//class vector2
//{
	vector2::vector2()
	{
		x = y = 0.0f;
	}
	vector2::vector2(const float& a,const float& b)
	{
		x=a;
		y=b;
	}
	vector2 vector2::operator+ (const vector2& b)
	{
		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;
	}
	vector2 vector2::operator-(const vector2& b)
	{
		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;
	}
	float vector2::operator* (const vector2& b)
	{
		return (x*b.x+y*b.y);
	}
	vector2 vector2::operator* (const float& b)
	{
		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;
	}
	vector2 vector2::operator/ (const float& b)
	{
		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);
	}
//}

verletPoint::verletPoint()
{
	vector2 pos (0.0f,0.0f);
	vector2 lastpos (0.0f,0.0f);
	vector2 accel (0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position, const vector2 & last_position, const vector2 & acceleration)
{
	pos=position;
	lastpos=last_position;
	accel=acceleration;
}
verletPoint::verletPoint(const vector2 & position,const vector2 & last_position)
{
	pos=position;
	lastpos=last_position;
	vector2 accel (0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position)
{
	pos=position;
	lastpos=position;
	vector2 accel (0.0f,0.0f);
}
void verletPoint::updateVerlet(const float timestep)
{
	vector2 temp=pos;
	pos+=pos-lastpos+accel*timestep*timestep;
	lastpos=temp;
}

circPhy::circPhy(const float center_x,const float center_y,const float radius, bool movable)
{
	M=radius*radius;
	R=radius;
	updateable=movable;
	center=verletPoint(vector2(center_x,center_y),vector2(center_x,center_y));
}
void circPhy::updateVerlet(const float timestep)
{
	if(updateable)
		center.updateVerlet(timestep);
}
bool circPhy::checkColl(circPhy* colObj)
{
	if(this!=colObj)
	{
		vector2 n = colObj->center.pos - center.pos;
		float dist=n.length();
		if(dist>=R+colObj->R)
			return false;
		n.normalize();
		coll.normal = n;
		coll.depth = R+colObj->R-dist;
		coll.bod = colObj;
		return true;
	}
	else return false;
}
bool circPhy::checkColl(circPhy* colObj,const float dist)
{
	vector2 n = colObj->center.pos - center.pos;
	if(dist>=R+colObj->R)
		return false;
	n.normalize(dist);
	coll.normal=n;
	coll.depth=R+colObj->R-dist;
	coll.bod=colObj;
	return true;
}
void circPhy::processColl()
{
	if(updateable)
	{
		if(coll.bod->updateable)
		{
			vector2 movement = coll.normal*coll.depth*.5f;
			center.pos-=movement;
			coll.bod->center.pos+=movement;
		}
		else
		{
			vector2 movement = coll.normal*coll.depth;
			center.pos-=movement;
		}
	}
}

+ Code Snippet

#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);
		vector2 operator+	(const vector2& b);	//vector addition
		vector2& operator+=	(const vector2& b);
		vector2 operator-	(const vector2& b);	//vector subtraction
		vector2& operator-=	(const vector2& b);
		vector2 operator*	(const float& b);	//vector scalar multiplication
		vector2& operator*=	(const float& b);
		vector2 operator/	(const float& b);	//vector scalar division
		vector2& operator/=	(const float& b);
		float operator*		(const vector2& b);	//dot product
		void normalize();
		void normalize(const float& dist);
		vector2 normalized();
		vector2 normalized(const float& dist);
		float length();
};
class verletPoint
{   
    public:   
        vector2 pos;
        vector2 lastpos;
        vector2 accel;
        verletPoint();
        verletPoint(const vector2& position,const vector2& last_position,const vector2& acceleration);
        verletPoint(const vector2& position,const vector2& last_position);
        verletPoint(const vector2& position);
        void updateVerlet(const float timestep);
};
struct collisionDataTemp
{
    float depth;
    vector2 normal;
    circPhy* bod;
};
class circPhy
{
    public:
        float M;
        float R;
		bool updateable;
        verletPoint center;
        collisionDataTemp coll;
        circPhy(){}
        circPhy(const float center_x,const float center_y,const float radius, bool moveable);
        void updateVerlet(const float timestep);
        bool checkColl(circPhy* colObj);
        bool checkColl(circPhy* colObj,const float dist);
        void processColl();
};


void DarkGDK(void)
{
	dbSetDisplayMode(640,480,32);
	int number_of_bodies = 5;
	int lastclick=0;
	d3dInit();
	dbSyncOn();
	dbSyncRate(60);
	dbRandomize(dbTimer());
	vector<circPhy> r;
	int moveobj=0;
	for(int n=0;n<number_of_bodies;n++)
		r.push_back(circPhy(dbRnd(dbScreenWidth()),dbRnd(dbScreenHeight()),dbRnd(5)+5,true));
	while ( LoopGDK ( ) )
	{
		dbCLS();
		dbInk(dbRGB (0, 0, 255),dbRGB (0, 0, 255));
		if(!dbSpaceKey())
			for(int n=0;n<number_of_bodies;n++)
			{
				if(dbControlKey()==1)
				{
					r[n].center.accel=vector2();
					r[n].center.lastpos=r[n].center.pos;
				}
				else
				{
					vector2 accel=vector2(0.0f,0.0f);
					for(int i=0;i<number_of_bodies;i++)
					{
						if(i!=n && r[i].updateable)
						{
							vector2 temp = r[i].center.pos-r[n].center.pos;
							float dist = temp.length();
							temp.normalize(dist);
							//now dist is the distance between the circle and the mouse, and accel is a unit vector
							dist=r[n].M*r[i].M*.001f/(dist*dist);
							//dist=(dist>.1f)?.1f:dist;
							temp*=dist;
							accel+=temp;
						}
					}
					r[n].center.accel=accel;
				}
				float* t=&r[n].center.pos.x;
				*t=(*t>=0)?  ((*t<=dbScreenWidth())?*t:dbScreenWidth())  :0;
				t=&r[n].center.pos.y;
				*t=(*t>=0)?  ((*t<=dbScreenHeight())?*t:dbScreenHeight())  :0;
				for(int i=0;i<number_of_bodies;i++)
				{
					if(r[n].checkColl(&r[i]))
					{r[n].processColl();}
				}
				r[n].updateVerlet(2);
			}


		if(dbMouseClick() && lastclick!=dbMouseClick())
		{
			if(dbMouseClick()==1)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,true));
				number_of_bodies++;
			}
			if(dbMouseClick()==2)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(5)+5,false));
				number_of_bodies++;
			}
			if(dbMouseClick()==4)
			{
				r.push_back(circPhy(dbMouseX(),dbMouseY(),dbRnd(100)+5,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()));
		dbCircle(365,351,63);
		dbSync();
	}
	return;
}
//class vector2
//{
	vector2::vector2()
	{
		x = y = 0.0f;
	}
	vector2::vector2(const float& a,const float& b)
	{
		x=a;
		y=b;
	}
	vector2 vector2::operator+ (const vector2& b)
	{
		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;
	}
	vector2 vector2::operator-(const vector2& b)
	{
		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;
	}
	float vector2::operator* (const vector2& b)
	{
		return (x*b.x+y*b.y);
	}
	vector2 vector2::operator* (const float& b)
	{
		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;
	}
	vector2 vector2::operator/ (const float& b)
	{
		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);
	}
//}


	
verletPoint::verletPoint()
{
	vector2 pos (0.0f,0.0f);
	vector2 lastpos (0.0f,0.0f);
	vector2 accel (0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position, const vector2 & last_position, const vector2 & acceleration)
{
	pos=position;
	lastpos=last_position;
	accel=acceleration;
}
verletPoint::verletPoint(const vector2 & position,const vector2 & last_position)
{
	pos=position;
	lastpos=last_position;
	vector2 accel (0.0f,0.0f);
}
verletPoint::verletPoint(const vector2 & position)
{
	pos=position;
	lastpos=position;
	vector2 accel (0.0f,0.0f);
}
void verletPoint::updateVerlet(const float timestep)
{
	vector2 temp=pos;
	pos+=pos-lastpos+accel*timestep*timestep;
	lastpos=temp;
}


circPhy::circPhy(const float center_x,const float center_y,const float radius, bool movable)
{
	M=radius*radius;
	R=radius;
	updateable=movable;
	center=verletPoint(vector2(center_x,center_y),vector2(center_x,center_y));
}
void circPhy::updateVerlet(const float timestep)
{
	if(updateable)
		center.updateVerlet(timestep);
}
bool circPhy::checkColl(circPhy* colObj)
{
	if(this!=colObj)
	{
		vector2 n = colObj->center.pos - center.pos;
		float dist=n.length();
		if(dist>=R+colObj->R)
			return false;
		n.normalize();
		coll.normal = n;
		coll.depth = R+colObj->R-dist;
		coll.bod = colObj;
		return true;
	}
	else return false;
}
bool circPhy::checkColl(circPhy* colObj,const float dist)
{
	vector2 n = colObj->center.pos - center.pos;
	if(dist>=R+colObj->R)
		return false;
	n.normalize(dist);
	coll.normal=n;
	coll.depth=R+colObj->R-dist;
	coll.bod=colObj;
	return true;
}
void circPhy::processColl()
{
	if(updateable)
	{
		if(coll.bod->updateable)
		{
			vector2 movement = coll.normal*coll.depth*.5f;
			center.pos-=movement;
			coll.bod->center.pos+=movement;
		}
		else
		{
			vector2 movement = coll.normal*coll.depth;
			center.pos-=movement;
		}
	}
}

(also note it uses cloggy's d3d dll, because built in circle-and-stuff commands SUCK and don't work on vista)

left click to add a regular object
right click to add a nonmovable object with 0 mass (acts as a wall)
middle click to create a draggable collision object (0 mass)
spacebar to pause the simulation
control to zero all velocities
hhhhhh hhhhhh
hhhhh hhhhh
hhhh hhhh

Watch out while holding the space bar!!!!! since it's just a quick once-over I don't check for divide by zero! If you create two objects in the exact same spot while paused, then their positions will be undefined! If their positions are undefined, and you check the distance against that point, then it will return undefined, and every single object will have a position of undefined!

Screenie:

[edit]
and I'm posting this here because it's too small to be announced
[edit2]
right, forgot to upload the exe, here it is:
[href]http://www.filefactory.com/file/b0c214a/n/GDKtest.exe [/href]

Attachments

Login to view attachments

Back to top

Profile PM Email

Dr Tank

15

Years of Service

User Offline

Joined: 1st Apr 2009

Location: Southampton, UK

Posted: 21st Mar 2010 04:35

Link

Looks cool. Maybe i'll try to read it someday.

I'm not really fluent in C++ yet.

When I'm doing forces simulations, I bin the objects into a grid, and only do interactions within some cutoff range. This is far from perfect though. Maybe what you could do is, for stuff in grid squares beyond the cutoff range, treat the mass inside grid squares as all concentrated at the centre, or calculate the centre of mass for it, or something similar to that. I imagine if you're quoting names of authors then they'll have thought about all this stuff. Will be interesting to see what kind of speed increases you can get.

I don't have GDK, so maybe you could upload an exe for people to try.

SUPER WEBSITE 22G

Back to top

Profile PM

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 21st Mar 2010 04:48

Link

uploaded it, I'm also posting a video on youtube, might take a while if my computer is retarded.

Back to top

Profile PM Email

Dr Tank

15

Years of Service

User Offline

Joined: 1st Apr 2009

Location: Southampton, UK

Posted: 24th Mar 2010 18:55 Edited at: 24th Mar 2010 18:56

Link

Works well. I made a collision system for round objects, but when there were loads pushing up against eachother, it didn't work so well. Your system doesn't seem to suffer from this nearly as much. Hopefully i can learn something from your code. Thanks.

For about 200 objects, i get about 20fps. You can do better than order n^2, dude. Put in a grid system or something!

SUPER WEBSITE 22G

Back to top

Profile PM

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 26th Mar 2010 03:57

Link

I really want to work on this but... I'll be in japan for a week, and ma leaving 12 hours from now. Unless I program a glitch-free quadtree (what I'm planning to do) after I wake up before I leave... I won't be able to work on this for a week.
Also, this is the relevant code for the collision system:

+ Code Snippet

bool circPhy::checkColl(circPhy* colObj)
{
	if(this!=colObj)
	{
		vector2 n = colObj->center.pos - center.pos;
		float dist=n.length();
		if(dist>=R+colObj->R)
			return false;
		n.normalize();
		coll.normal = n;
		coll.depth = R+colObj->R-dist;
		coll.bod = colObj;
		return true;
	}
	else return false;
}
bool circPhy::checkColl(circPhy* colObj,const float dist)
{
	vector2 n = colObj->center.pos - center.pos;
	if(dist>=R+colObj->R)
		return false;
	n.normalize(dist);
	coll.normal=n;
	coll.depth=R+colObj->R-dist;
	coll.bod=colObj;
	return true;
}
void circPhy::processColl()
{
	if(updateable)
	{
		if(coll.bod->updateable)
		{
			vector2 movement = coll.normal*coll.depth*.5f;
			center.pos-=movement;
			coll.bod->center.pos+=movement;
		}
		else
		{
			vector2 movement = coll.normal*coll.depth;
			center.pos-=movement;
		}
	}
}

I think after I finish with this, I'll work on a 3d version. Should be awesome, and still really speedy

Back to top

Profile PM Email

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 12th Jul 2010 11:42

Link

took a look at this again and... Have you noticed how you can't really get a cluster to orbit a cluster? That's cuz clusters slow down after a while, and they don't obey conservation of momentum X_X

Fix'n right now. I think I gots me a good equation to determine velocity but... it's the AMs and I have to be up in 5 hrs.

If anyone wants to know, where P1,P2 are positions, V1,V2 are velocities, M1,M2 are masses, u is coefficient of restitution, and X is the normalized vector of P2-P1, or X=||P2-P1||, also . is vector dot product and * is scalar multiplication
V1'=V1-(V1.X)*X+((u*M1*(V1.X-V2.X)+M1*V1.X+M2*V2.X)/(M1+M2))*X
That should determine the new velocity of the verlet circle while conserving momentum.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 13th Jul 2010 23:38 Edited at: 14th Jul 2010 05:54

Link

Well... that vector method looks too computationally intense, and with further thought, I think the whole idea of having rigid circles colliding is bad.

Instead, I'm going to treat the particles not as circles, but as point "molecules" - attracting each other, until they get close, and then sharply repelling each other.

Normally, the equation for force over distance is newton's law of gravitation:

whose graph looks like this:

lets define F1(d) to be that function, in other words
F1(d)=G*m1*m2/d^2

Now... I'm not trying to make a true to life simulation, so I don't really need true to life equations. IMHO as long as my simulation obeys newtons laws, it's accurate enough. So, now we want an equation for a repelling force (in this case, a negative force repels, a positive force attracts), with an equally sharp drop off as F1(d) as the distance approaches 0. Basically:
F2(d)=C-M1*M2/d^2
where C is a constant. Now, this equation should also have the property that when d=r1+r2 (IE the circles are just touching), they are at an equilibrium, and the force between them is 0, IE
F2(r1+r2)=0
solving for C gives
C=M1*M2/(r1+r2)^2
so
F1(d)=G*m1*m2/d^2
F2(d)=M1*M2(1/(r1+r2)^2-1/d^2)
the graph of those two equations looks like this:

The lowest of the two values at any point seems to be a good graph, so: F(d)=min(F1(d),F2(d))=min(G*m1*m2/d^2,M1*M2(1/(r1+r2)^2-1/d^2))
graphed:

With two large radius/low density particles, they would probably bounce off each other repeatedly. With two extremely small radius/very high density particles, it would probably be vibrating or haywire collisions. Both of these situations are visually unappealing (circles slowly bouncing or going crazy, intruding into each others' boundaries), but I think with hundreds or thousands of such particles, better structures will form. There is also no need for collision detection, because everything works itself out. If structures don't form, and everything just repels everything else too much, velocity dampening could be added in the direction of the collision.

In vector notation, where P1 and P2 are vectors of the position of one such particle and another:
Force on particle 1 += (P1-P2)/( |P1-P2|)*F( |P1-P2|)
Force on particle 2 -= (P1-P2)/( |P1-P2|)*F( |P1-P2|)

I hafta go right now, but hopefully I'll have all this implemented soon.

[edit]
o rite, all graphs are force over separating distance of two particles.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Neuro Fuzzy

16

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 14th Jul 2010 06:59 Edited at: 14th Jul 2010 10:39

Link

Even better graph! Its pretty much as easy to make that function a continuous curve:
F(d)=G*m1*m2/d^2*(1-(r1+r2)^2/d^2)
Try it out! As the radii approach zero, or the larger x gets, the closer the function is to G*m1*m2/d^2

global pixels_per_unit as float
pixels_per_unit=100
sync on
r1 as float
r1=0
r2 as float
r2=0
m1 as float
m1=0
m2 as float
m2=0
G as float
G=0
set display mode 1280,800,32

do
   cls
   from#=-4
   ato#=8
  ` F1(d)=G*m1*m2/d^2
	`F2(d)=M1*M2(1/(r1+r2)^2-1/d^2)
	inc r1,(keystate(16)-keystate(30))/100.0
	inc r2,(keystate(17)-keystate(31))/100.0
	inc m1,(keystate(18)-keystate(32))/100.0
	inc m2,(keystate(19)-keystate(33))/100.0
	inc G,(keystate(20)-keystate(34))/100.0
	text 0,0,"r1="+str$(r1)
	text 0,16,"r2="+str$(r2)
	text 0,32,"m1="+str$(m1)
	text 0,48,"m2="+str$(m2)
	text 0,64,"G="+str$(G)
	text 0,80,"density1="+str$(m1/(3.14159265*r1*r1))
	text 0,96,"density2="+str$(m2/(3.14159265*r2*r2))
   for x#=from# to ato# step 0.005
      dot x#*pixels_per_unit+screen width()/2,-G*m1*m2/x#^2*pixels_per_unit+screen height()/2
   next x#
   for x#=from# to ato# step 0.005
      dot x#*pixels_per_unit+screen width()/2,-(G*m1*m2/x#^2)*(1-(r1+r2)^2/x#^2)*pixels_per_unit+screen height()/2
   next x#
   drawGraphPaper(pixels_per_unit,screen width()/2,screen height()/2,4)
   sync
loop

end
function minimum(a as float, b as float)
	if a<b then exitfunction a
endfunction b
function drawGraphPaper(PPU as integer, midx as integer, midy as integer, subUnits as integer)
   if subunits<=0 then subunits=1
   if PPU<=0 then ppu=40
   subDist as float
   subDist=PPU*1.0/subUnits
   counter1 as float =
   counter1=-PPU
   counter2 as float = 0
   w=screen width()
   h=screen height()
   ink rgb(0,0,40)
   line 0,h/2,w,h/2
   ink rgb(0,40,0)
   line w/2,0,w/2,h
   while (counter1+midx<w) or (midx-counter1>0)
      inc counter1, PPU
      counter2=0
      while (counter2<PPU)
         inc counter2, subDist
         ink rgb(150,150,150),0
         line midx+counter1+counter2,midy-2,midx+counter1+counter2,midy+2
         line midx-counter1+counter2,midy-2,midx-counter1+counter2,midy+2
      endwhile
      ink rgb(255,255,255),0
      line midx+counter1,midy-5,midx+counter1,midy+5
      line midx-counter1,midy-5,midx-counter1,midy+5

endwhile
   counter1=-PPU
   counter2=0
   while (counter1+midy<w) or (midy-counter1>0)
      inc counter1, PPU

counter2=0
      while (counter2<PPU)
         inc counter2, subDist
         ink rgb(150,150,150),0
         line midx-2,midy+counter1+counter2,midx+2,midy+counter1+counter2
         line midx-2,midy-counter1-counter2,midx+2,midy-counter1-counter2
      endwhile
      ink rgb(255,255,255),0
      line midx-4,midy+counter1,midx+4,midy+counter1
      line midx-4,midy-counter1,midx+4,midy-counter1
   endwhile
endfunction

function ln(x as float)
   n=16
   ret#=3.141593/(2*agm(1,4/(x*2^n)))-n*.6931472
endfunction ret#

`arithmetic-geometric mean
function agm(x as float, y as float)
   a#=(x+y)/2
   g#=sqrt(x*y)
   for n=1 to 10
      at#=(a#+g#)/2
      gt#=sqrt(a#*g#)
      a#=at#
      g#=gt#
   next n
   ret#=(a#+g#)/2
endfunction ret#

+ Code Snippet

global pixels_per_unit as float
pixels_per_unit=100
sync on
r1 as float
r1=0
r2 as float
r2=0
m1 as float
m1=0
m2 as float
m2=0
G as float
G=0
set display mode 1280,800,32

do
   cls
   from#=-4
   ato#=8
  ` F1(d)=G*m1*m2/d^2
	`F2(d)=M1*M2(1/(r1+r2)^2-1/d^2)
	inc r1,(keystate(16)-keystate(30))/100.0
	inc r2,(keystate(17)-keystate(31))/100.0
	inc m1,(keystate(18)-keystate(32))/100.0
	inc m2,(keystate(19)-keystate(33))/100.0
	inc G,(keystate(20)-keystate(34))/100.0
	text 0,0,"r1="+str$(r1)
	text 0,16,"r2="+str$(r2)
	text 0,32,"m1="+str$(m1)
	text 0,48,"m2="+str$(m2)
	text 0,64,"G="+str$(G)
	text 0,80,"density1="+str$(m1/(3.14159265*r1*r1))
	text 0,96,"density2="+str$(m2/(3.14159265*r2*r2))
   for x#=from# to ato# step 0.005
      dot x#*pixels_per_unit+screen width()/2,-G*m1*m2/x#^2*pixels_per_unit+screen height()/2
   next x#
   for x#=from# to ato# step 0.005
      dot x#*pixels_per_unit+screen width()/2,-(G*m1*m2/x#^2)*(1-(r1+r2)^2/x#^2)*pixels_per_unit+screen height()/2
   next x#
   drawGraphPaper(pixels_per_unit,screen width()/2,screen height()/2,4)
   sync
loop

end
function minimum(a as float, b as float)
	if a<b then exitfunction a
endfunction b
function drawGraphPaper(PPU as integer, midx as integer, midy as integer, subUnits as integer)
   if subunits<=0 then subunits=1
   if PPU<=0 then ppu=40
   subDist as float
   subDist=PPU*1.0/subUnits
   counter1 as float =
   counter1=-PPU
   counter2 as float = 0
   w=screen width()
   h=screen height()
   ink rgb(0,0,40)
   line 0,h/2,w,h/2
   ink rgb(0,40,0)
   line w/2,0,w/2,h
   while (counter1+midx<w) or (midx-counter1>0)
      inc counter1, PPU
      counter2=0
      while (counter2<PPU)
         inc counter2, subDist
         ink rgb(150,150,150),0
         line midx+counter1+counter2,midy-2,midx+counter1+counter2,midy+2
         line midx-counter1+counter2,midy-2,midx-counter1+counter2,midy+2
      endwhile
      ink rgb(255,255,255),0
      line midx+counter1,midy-5,midx+counter1,midy+5
      line midx-counter1,midy-5,midx-counter1,midy+5

   endwhile
   counter1=-PPU
   counter2=0
   while (counter1+midy<w) or (midy-counter1>0)
      inc counter1, PPU

      counter2=0
      while (counter2<PPU)
         inc counter2, subDist
         ink rgb(150,150,150),0
         line midx-2,midy+counter1+counter2,midx+2,midy+counter1+counter2
         line midx-2,midy-counter1-counter2,midx+2,midy-counter1-counter2
      endwhile
      ink rgb(255,255,255),0
      line midx-4,midy+counter1,midx+4,midy+counter1
      line midx-4,midy-counter1,midx+4,midy-counter1
   endwhile
endfunction


function ln(x as float)
   n=16
   ret#=3.141593/(2*agm(1,4/(x*2^n)))-n*.6931472
endfunction ret#

`arithmetic-geometric mean
function agm(x as float, y as float)
   a#=(x+y)/2
   g#=sqrt(x*y)
   for n=1 to 10
      at#=(a#+g#)/2
      gt#=sqrt(a#*g#)
      a#=at#
      g#=gt#
   next n
   ret#=(a#+g#)/2
endfunction ret#

[edit]
well neither of those worked at all, and everything goes haywire when I try to implement that collision conservation of momentum stuff.

Is't life, I ask, is't even prudence, to bore thyself and bore thy students?

Back to top

Profile PM Email

Sorry your browser is not supported!

Code Snippets / [GDK] O(n^2) n-body gravitational simulation

Attachments