C++

Forum

 
gravity simulator

closed account (zwA4jE8b)
Hey guys, I am trying to write a gravity simulator in c++ and sdl.

Can anyone, who knows physics, help me with the main code.
I am doing this mostly from memory and the help of my physics 131 book. My main problem, I believe, is in the merge function. If you run it, you will see that sometimes two colliding bodies will have a greater velocity than they really should. Any input will be appreciated

For collision .cpp and .h
SDL_collide
http://sdl-collide.sourceforge.net/

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
 
/*
 * 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 += aX * DELTA_TIME;
			vY += aY * DELTA_TIME;

			x += vX * DELTA_TIME;
			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;
}
KRAkatau (197)
It looks like your new velocity computation is overly complicated. If you think about the momentum you will see that the momentum of the resulting body should be equal to the sum of momentums of the original bodies. So msum vsum = m1 v1 + m2 v2. So to compute the velocity of resulting object, compute the momentums of original objects, sum them up (as vectors) and divide by the mass of resulting object. This should do the job.
Last edited on
Major Tom (55)
There are two ways to simulate gravity. You can use a functional gravity simulator for projectile motion, that will generate a function for the trajectory, or you can go the real way an use a numerical integrator.
Projectile Motion:
http://electron9.phys.utk.edu/phys135d/modules/m3/projectile%20motion.htm

Numerical Integration of Forces:
http://gafferongames.com/game-physics/integration-basics/

I would highly recommend the numerical integration method, but it requires knowledge of basic calculus.
Glenn Fiedler's articles on integration and physics are fantastic; if your looking to build an engine, read those.
Last edited on
closed account (zwA4jE8b)
thank you two. Im going to go work on it. Ill let post how it turns out.
great article btw
rapidcoder (1010)
Nice article, until to this point:


Rather than introduce you to the vast array of different integrators that exist, I’m going to cut to the chase and go straight to the best. This integrator is called the Runge Kutta order 4 integrator aka RK4. This is the standard integrator used for numerical integration these days and is sufficiently accurate for just about anything required in game physics, given an appropriate timestep.


This is utter nonsense. There is no thing like the "best" integrator. RK4 is good only for very small timesteps (read it: it is slow) and for larger timesteps it suffers from huge instabilities. In games it is much better to use a stable integrator which is slightly inaccurate (even the inaccurate Euler one is good for many problems), than an accurate integrator that suffers from instability problems. A good second order integrator (like BDF23) with dynamic timestep beats pants off RK4 IMHO.
Major Tom (55)
I agree. But for a simple gravity program, a Euler integrator would be the easiest to implement, and produce decent enough results. A backwards differentiation formula for a gravity problem is slightly overkill xD
Last edited on
closed account (zwA4jE8b)
Heya, I decided to start with a simpler program.

This one allows the user to click anywhere on the screen and a planet is spawned. It is then sucked into the 'blackhole' at the center of the screen. You can click on the main planet to select it and move it around. Its kind of cool becuase some planets will start to orbit.

I am however getting a 'floating point exception' error in 1) when more than 139 is specified as a command line argument and 2) when clicking planets, the count gets to 209. Can anyone tell me what is happening here?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
 
/*
 * 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 = 1; 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, 1))
			{
				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;
}
Topic archived. No new replies allowed.