How Does Euler Integration Affect Collision Handling in a C++ Gravity Simulator?

[creativemfs]

Hi everyone, I am new here, nice to meet you.

For fun I am trying to make a gravity simulator in c++ and using SDL.

I believe there is a problem in my merge function when I calculate the new velocity and trajectory of two colliding bodies.

Also, does my main code look accurate. The part where I calculate the gravitational force, determine acceleration, and sum the acceleration from one body to all the other bodies.

This is my first time attempting a project like this and my professors are gone for the summer so I am kind of in the dark.

I am using euler integration for now, it is all I know and can implement.

Thank you,
Michaela

if you want to run the code
the collision code is here
http://sdl-collide.sourceforge.net/

/*
 * main.cpp
 *
 *  Created on: Jul 28, 2012
 *      Author: michaela
 *      Fnet = (G * m1 * m2) / r^2
 */

#include "SDL_collide.h"
#include <SDL/SDL.h>
#include <SDL/SDL_gfxPrimitives.h>
#include <list>
#include <math.h>
#include <iostream>

const int SCREEN_HEIGHT = 600;
const int SCREEN_WIDTH = 800;
const int SCREEN_BPP = 32;

const double G = 6.67e-11;
const double DELTA_TIME = 0.1;

//Hardcoded number of bodies, used if no argument is specified
int bodies = 100;

struct celestial_body
{
		int count;
		double mass;
		double radius;
		double x, y;
		double vX, vY;
		double aX, aY;
		SDL_Color color;

		celestial_body(int newCount)
		{
			count = newCount;

			//Random mass
			mass = rand() % 100000000;

			radius = 0.5;

			//Random position on screen
			x = fmod((double)rand(), (double)SCREEN_WIDTH);
			y = fmod((double)rand(), (double)SCREEN_HEIGHT);

			vX = vY = 0;

			//Random color and avoid very dark colors
			color.r = (int)mass % (rand() % 10000) % 255;
			color.g = (int)mass % (rand() % 10000) % 255;
			color.b = (int)mass % (rand() % 10000) % 255;

			if(color.r <= 20 && color.g <= 20 && color.b <= 20)
			{
				color.r += 25;
				color.g += 25;
				color.b += 25;
			}
		}
		void draw(SDL_Surface* screen)
		{
			filledCircleRGBA(screen, x, y, radius, color.r, color.g, color.b, 255);
		}
		void update()
		{
			vX = vX + (aX * DELTA_TIME);
			vY = vY + (aY * DELTA_TIME);

			x = x + (vX * DELTA_TIME);
			y = y + (vY * DELTA_TIME);
		}
		void print()
		{
			std::cout << "Mass: " << mass << "\nX: " << x << "\nY: " << y <<
						 "\nVx : " << vX << "   Vy: " << vY <<
					     "\nColor: " << (int)color.r << " " << (int)color.g << " " << (int)color.b << "\n\n";
		}
};

void merge(celestial_body &a, celestial_body &b)
{
	//Radius and color merge is arbitrary (for now)
	a.radius += b.radius;
	a.color.r = (a.color.r + b.color.r) / 2;
	a.color.g = (a.color.g + b.color.g) / 2;
	a.color.b = (a.color.b + b.color.b) / 2;

	//Total velocity for a and b
	double vA = sqrt((a.vX * a.vX) + (a.vY * a.vY));
	double thetaA = atan2(a.vY, a.vX);
	if(thetaA < 0)
		vA *= -1.0;

	double vB = sqrt((b.vX * b.vX) + (b.vY * b.vY));
	double thetaB = atan2(b.vY, b.vX);
	if(thetaB < 0)
		vB *= -1.0;

	//New velocity of collided bodies
	double velocity = ((a.mass * vA) + (b.mass * vB)) / (a.mass + b.mass);

	//Angle of trajectory
	double theta = atan2(((a.mass * a.vY) + (b.mass * b.vY)), ((a.mass * a.vX) + (b.mass * b.vX)));

	std::cout  << "vA: " << vA << "   vB: " << vB << "\n" << velocity << "  :  " << theta << std::endl;

	//Vector components of velocity
	a.vX = velocity * cos(theta);
	a.vY = velocity * sin(theta);

	//Masses are added together (mass final)
	a.mass += b.mass;

	a.print();
}

void cleanup(std::list<celestial_body> &cb)
{
	cb.clear();
	SDL_Quit();
}

int main(int argc, char** argv)
{
	//Process arguments
	if(argc == 2)
	{
		int tempBodies = atoi(argv[1]);
		if(tempBodies > 0)
			bodies = tempBodies;
	}

	SDL_Surface* screen;
	SDL_Event event;

	std::list<celestial_body> cb;
	std::list<celestial_body>::iterator it;
	std::list<celestial_body>::iterator ij;

	double distance, diffX, diffY, tempAx, tempAy;
	double Fnet;	//Net Force on body
	double theta;	//Angle between two points in 2-D space
	double accel;	//Net acceleration of body

	bool quit = false;

	if(SDL_Init(SDL_INIT_EVERYTHING) == -1)
		quit = true;

	SDL_WM_SetCaption("Gravity Simulator", NULL);
	SDL_ShowCursor(1);

	screen = SDL_SetVideoMode(SCREEN_WIDTH, SCREEN_HEIGHT, SCREEN_BPP, SDL_HWSURFACE | SDL_DOUBLEBUF);
	if(screen == NULL)
		quit = true;

	for(int i = 0; i < bodies; ++i)
	{
		celestial_body newBody(i);
		cb.push_back(newBody);
		newBody.print();
	}

	std::cout << "Number of bodies: " << bodies << "\n\n";

	while(!quit)
	{
		while(SDL_PollEvent(&event))
		{
			if(event.type == SDL_QUIT)
				quit = true;
			if(event.type == SDL_KEYDOWN)
				if(event.key.keysym.sym == SDLK_ESCAPE)
					quit = true;
		}

//		SDL_FillRect(screen, &screen->clip_rect, 0x000000);

		for(it = cb.begin(); it != cb.end(); ++it)
			(*it).draw(screen);

		SDL_Flip(screen);

		for(it = cb.begin(); it != cb.end(); ++it)
		{
			tempAx = tempAy = 0;
			for(ij = cb.begin(); ij != cb.end(); ++ij)
			{
				//Simply makes sure a body is not calculated against itself
				if((*it).count != (*ij).count)
				{
					diffX = (*it).x - (*ij).x;
					diffY = (*it).y - (*ij).y;

					if(!SDL_CollideBoundingCircle((*it).x, (*it).y, (*it).radius, (*ij).x, (*ij).y, (*ij).radius, 0))
					{
						//Determine the distance between two bodies
						distance = sqrt( (diffX * diffX) + (diffY * diffY) );

						//Determine the net gravitational force between the two
						Fnet = ((G * (*it).mass * (*ij).mass)) / (distance * distance);

						//Determine the angle from a to b in 2-Dimensional space
						theta = atan2(diffY, diffX);

						//Determine the total acceleration
						accel = Fnet / (*it).mass;

						//Find acceleration from vector components
						//Add them together for each affecting body
						tempAx += -(accel * cos(theta));
						tempAy += -(accel * sin(theta));
					}
					else
					{
						if((*it).mass > (*ij).mass)
						{
							merge((*it), (*ij));
							ij = cb.erase(ij);
						}
						else
						{
							merge((*ij), (*it));
							it = cb.erase(it);
						}
						std::cout << "Number of bodies: " << cb.size() << "\n";
					}
				}
			}
			(*it).aX = tempAx;
			(*it).aY = tempAy;
		}

		for(it = cb.begin(); it != cb.end(); ++it)
			(*it).update();
	}

	cleanup(cb);

	return 0;
}

 

[chiro]

Hey creativemfs and welcome to the forums.

I've noticed you are calculating velocity vectors but you are not using forces per se: in other words, instead of introducing forces that are updated at every time interval and introduced upon either a collision event or by some non-contact force (like gravity) you are modifying the velocity vector directly.

Usually what happens in a physics engine is that the velocity vector doesn't change unless a force is applied (just like Newtons Laws of motion). So in a contact situation you can calculate the force on each body and then use a numerical integrator (like Euler's method) to use the force over some small dt to get the new velocity and then use this new velocity to update the position of said object.

So in short: you have three things for each object (vectors): position, velocity, acceleration. You only modify acceleration (unless you want to say do debugging, initialization or to hack your variables) and any introduced acceleration gets applied immediately within the time delta. The integrator will take the acceleration components and apply them within the time-delta and then set the acceleration back to zero. If a new event is called that applies more acceleration (like another collision or by having a non-contact force like gravity) then the next time step will incorporate that new force.

This is the basic idea of how a Newtonian physics engine works: you typically look at forces in one form or another and then you apply these forces to generate velocity changes and finally changes in displacement.

The question then remains: what kind of collision do you want? It seems like you were just averaging the two velocities together with the corresponding masses to be used as weights.

If you want to retain this, then you should calculate the impulse and then consider that relative to the time frame (i.e. the time delta in your simulation) to get the force. (So if your impulse is I then calculate I/dt where dt is your delta to get the force).

Now you take this force and then apply it to your object: I noticed you've merged objects together so you can just delete the other one and apply the force to the one in question.

In terms of masses, you don't supply these to the numerical integrator: the only thing that the integrator should use is the forces: you can use the mass to create forces, but ultimately it's the job to just look at forces and how these translate into positions at each time-step.

 

[Mark44]

Michaela,
Thanks for putting your code inside HTML code tags! It's rare that new members think to do this.

 

[creativemfs]

thank you for your reply, I did make this modified version. I decided starting with a simpler program will be best.

This is more like a black hole that sucks everything in.

So chiro, basically I want two bodies that collide to stick together (mass adds together) and travel at the new velocity in x and y. So I don't know if I need to even calculate a 'theta' because an individual x and y velocity will be traveling at some angle.

Can you give me an equation, or some code examples, for calculating the velocity of the collided bodies? How do I code impulse?

I have only taken one physics class in college, but I remember a little about impulse, I will look it up in my book and see if I can figure out where to go.

oh and your welcome mark. I spend a lot of time on cplusplus.com so code tags are nothing new to me :)

/*
 * main.cpp
 *
 *  Created on: Jul 28, 2012
 *      Author: michaela
 *      Fnet = (G * m1 * m2) / r^2
 *
 *      Optional parameters ./gravity_simulator [--random] [num bodies]
 *
 *      As of now, only the main body affects the other bodies, they do not affect each other
 *      Click on the main body and move the mouse to move it, click again to release
 */

#include "SDL_collide.h"
#include <SDL/SDL.h>
#include <SDL/SDL_gfxPrimitives.h>
#include <list>
#include <math.h>
#include <iostream>

const int SCREEN_HEIGHT = 600;
const int SCREEN_WIDTH = 800;
const int SCREEN_BPP = 32;

const double G = 6.67e-1;		//Scaled down by e-10
const double DELTA_TIME = 0.1;	//Arbitrary
const double RADIUS = 3.0;

//Hardcoded number of bodies, used if no argument is specified
int bodies = 100;

struct celestial_body
{
		int count;
		double mass;
		double radius;
		double x, y;
		double vX, vY;
		double aX, aY;
		SDL_Color color;

		celestial_body(int newCount = 0, double newRadius = 0.5, double newMass = 0, int newX = 0, int newY = 0)
		{
			count = newCount;

			if(newMass == 0)
				mass = (rand() % 1000) + 1;
			else
				mass = newMass;

			radius = newRadius;

			if(newX == 0)
				x = fmod((double)rand(), (double)SCREEN_WIDTH);
			else
				x = newX;
			if(newY == 0)
				y = fmod((double)rand(), (double)SCREEN_HEIGHT);
			else
				y = newY;

			vX = vY = 0;
			aX = aY = 0;

			//Random color and avoid very dark colors
			color.r = (int)mass % (rand() % 10000) % 255;
			color.g = (int)mass % (rand() % 10000) % 255;
			color.b = (int)mass % (rand() % 10000) % 255;

			if(color.r <= 20 && color.g <= 20 && color.b <= 20)
			{
				color.r += 25;
				color.g += 25;
				color.b += 25;
			}
		}
		void draw(SDL_Surface* screen)
		{
			filledCircleRGBA(screen, x, y, radius, color.r, color.g, color.b, 255);
		}
		void update()
		{
			vX += aX * DELTA_TIME;
			vY += aY * DELTA_TIME;

			x += vX * DELTA_TIME;
			y += vY * DELTA_TIME;
		}
		void print()
		{
			std::cout << "Count: " << count << "\nMass: " << mass << "\nX: " << x << "\nY: " << y <<
					     "\nColor: " << (int)color.r << " " << (int)color.g << " " << (int)color.b << "\n\n";
		}
};

void spawn_body(std::list<celestial_body>& cb, int newX = 0, int newY = 0)
{
	static int count = 1;
	celestial_body new_body(count++, RADIUS, 0, newX, newY);
	cb.push_back(new_body);
	new_body.print();
}

void cleanup(std::list<celestial_body>& cb)
{
	cb.clear();
	SDL_Quit();
}

int main(int argc, char** argv)
{
	SDL_Surface* screen;
	SDL_Event event;

	std::list<celestial_body> cb;
	std::list<celestial_body>::iterator it;
	std::list<celestial_body>::iterator ij;

	double distance, diffX, diffY;
	double Fnet;	//Net Force on body
	double theta;	//Angle between two points in 2-D space
	double accel;	//Net acceleration of body

	bool quit = false;
	bool random = false;
	bool selected = false;

	//Process arguments
	if(argc == 3)
	{
		if(strcmp(argv[1], "--random") == 0)
			random = true;
		int tempBodies = atoi(argv[2]);
		if(tempBodies > 0)
			bodies = tempBodies;
	}

	if(SDL_Init(SDL_INIT_VIDEO) == -1)
		quit = true;

	SDL_WM_SetCaption("Gravity Simulator", NULL);
	SDL_ShowCursor(1);

	screen = SDL_SetVideoMode(SCREEN_WIDTH, SCREEN_HEIGHT, SCREEN_BPP, SDL_HWSURFACE | SDL_DOUBLEBUF);
	if(screen == NULL)
		quit = true;

	//Initialize main static body
	celestial_body main(0, 8.0, 1000, SCREEN_WIDTH / 2, SCREEN_HEIGHT / 2);
	main.print();

	if(random)
		for(int i = 0; i < bodies; ++i)
			spawn_body(cb);

	while(!quit)
	{
		while(SDL_PollEvent(&event))
		{
			if(event.type == SDL_QUIT)
				quit = true;
			if(event.type == SDL_KEYDOWN)
				if(event.key.keysym.sym == SDLK_ESCAPE)
					quit = true;
			if(event.type == SDL_MOUSEBUTTONDOWN)
			{
				if(SDL_CollideBoundingCircle(main.x, main.y, main.radius, event.button.x, event.button.y, 1, 0))
					selected = !selected;
				else
					spawn_body(cb, event.button.x, event.button.y);
			}
		}

		if(selected)
		{
			main.x = event.motion.x;
			main.y = event.motion.y;
		}

		SDL_FillRect(screen, &screen->clip_rect, 0x000000);

		main.draw(screen);
		for(it = cb.begin(); it != cb.end(); ++it)
			(*it).draw(screen);

		SDL_Flip(screen);

		for(it = cb.begin(); it != cb.end(); ++it)
		{
			diffX = main.x - (*it).x;
			diffY = main.y - (*it).y;

			//Determine the distance between two bodies
			distance = sqrt((diffX * diffX) + (diffY * diffY));

			//Determine the net gravitational force between the two
			Fnet = ((G * (*it).mass * main.mass)) / (distance * distance);

			//Determine the angle from a to b in 2-Dimensional space
			theta = atan2(diffY, diffX);

			//Determine the total acceleration
			accel = Fnet / (*it).mass;

			//Find acceleration from vector components
			(*it).aX = accel * cos(theta);
			(*it).aY = accel * sin(theta);

			//Check for collision between new_body and main
			//main gains the mass of colliding body
			if(SDL_CollideBoundingCircle(main.x, main.y, main.radius, (*it).x, (*it).y, (*it).radius, 0))
			{
				main.mass += (*it).mass;
				it = cb.erase(it);
			}
		}

		for(it = cb.begin(); it != cb.end(); ++it)
			(*it).update();

		SDL_Delay(10);
	}

	cleanup(cb);

	return 0;
}

 

[creativemfs]

I see that (F*dT) / M = dV but how do I get this to apply to two colliding bodies?

Also was my method of calculating the acceleration due to the gravitational pull of the other bodies and summing it correct?