Gravitas Simulator

[ejc.cryptography]

From gravitas.sourceforge.net:

Gravitas is a fully scriptable, real-time 3D gravity simulator. Scripts written in the integrated development environment (IDE) are interpreted by the simulator and allow users to create any imaginable planetary setup.

Once inside the simulated system, the user can move around in 3D space to observe the gravitational interactions of the celestial bodies. With the aid of the head-up display (HUD), users can track bodies and view live information about the simulation as it progresses.



Gravitas makes use of the powerful yet simple-to-use Lua scripting language, which offers the perfect introduction to programming for beginners, and extensive features for more experienced users. The rendering of simulations is handled by OpenGL, the industry standard interactive 3D graphics environment.

Gravitas is an Open Source application which means that you can download the full source code for the simulator. It is programmed and compiled using FreeBASIC.

• Gravitas website
• Download
• Documentation

The website also features a page which gives a LUA script to simulate the Solar System; it is updated from the NASA JPL Horizons system every 4 hours - Solar System Script.

[h4tt3n]

Wow! This looks very impressive, ejc. More than any other fb developer you're capable of giving your software a profesional look and feel.

As for the functionality of the code, I took a brief look at the gravitycalc sub, where it seems all gravitational interaction and integration is done, and have a few ideas for (huge) improvement.

(1) The overall structure of iterating through the planets can be optimized quite a lot. Instead of doing

Code: Select all

For i As Integer = 1 to Ubound(planet)
	for j As Integer = 1 to Ubound(planet)
		if j <> i
			''	calc gravity force
		end if
	next j
next i

I suggest that you do like this (notice the small but important difference)

Code: Select all

For i As Integer = 1 to Ubound(planet)-1
	for j As Integer = i+1 to Ubound(planet)
		''	calc gravity force
	next j
next i

In your version gravity is calculated for each planet pair twice, making the code much slower than needed (and perhaps straight wrong ?). On top of that the if-then statement in the inner loop takes up some cpu time aswell. You'll experience a noteworthy increase in performance by doing like I suggested. This is also called the "handshake-problem".

(2) I'd suggest that you get rid of all instances of the exponentiate operator ( ^ ), since I'm afraid they're rather cpu-heavy and un-optimized (please correct me if I'm wrong here). As an example, instead of writing

Code: Select all

D = (planet(i).x-planet(j).x)^2+(planet(i).y-planet(j).y)^2+(planet(i).z-planet(j).z)^2

just write

Code: Select all

D = (planet(i).x-planet(j).x*planet(i).x-planet(j).x)+(planet(i).y-planet(j).y*planet(i).y-planet(j).y)+(planet(i).z-planet(j).z*planet(i).z-planet(j).z)

And likewise in similar cases. This should improve performance a tad.

(3) You're using the simple, first order symplectic Euler integration method. This is ok for games, but imho it is not accurate enough for an otherwise impressive simulator like yours.
You should at least let the user pick a symplectic second order or fourth order integration method for higher precision and stability. Since this is a somewhat advanced topic, I'll just add a simple example on how to apply up to sixth and eighth order methods for you to read at your pleasure. Take it or leave it :-)

Good luck with the further development. I won't mind discussing this further with you.

(I'll add the examples as soon as I've dug them out )

Edit: Here they are...


vec2f.bi (vector library)

Code: Select all

''*******************************************************************************
''		
''		Freebasic 2d floating point vector library
''   	version 0.4b, may 2009, Michael "h4tt3n" Nissen, jernmager@yahoo.dk
''		
''   	function syntax:
''		
''   	(return type) (function name) (argument type (, ...))
''		
''	 	function list:
''		
''   	vector absolute          (vector)          - absolute value
''   	vector normal            (vector)          - normal vector
''   	vector normalised        (vector)          - normalised vector
''   	scalar magnitude         (vector)          - magnitude
''   	scalar magnitudesquared  (vector)          - magnitude squared
''   	scalar distance          (vector, vector)  - distance between vectors
''   	scalar distancesquared   (vector, vector)  - distance between vectors squared
''   	scalar dot               (vector, vector)  - dot product
''   	scalar dotnormal         (vector, vector)  - normal dot product
''   	vector project           (vector, vector)	 - project vector a on vector b
''		
''	 	all functions exist both in a member and non-member version.
''		
''*******************************************************************************

type float as single

''  2d float vector structure
type vec2f
	
  ''  variables
  as float x, y
	
  ''  constructor declarations
	declare constructor ()
  declare constructor (byval X as float, byval Y as float)
	
  ''  compound arithmetic member operator declarations
  declare operator += (byref rhs as vec2f)
  declare operator -= (byref rhs as vec2f)
  declare operator *= (byref rhs as vec2f)
  declare operator *= (byref rhs as float)
  declare operator /= (byref rhs as float)
  declare operator let (byref rhs as vec2f)
	
	''  member function declarations
  declare function absolute() as vec2f
  declare function normal() as vec2f 
  declare function normalised() as vec2f
  declare function magnitude() as float
  declare function magnitudesquared() as float
  declare function distance(byref rhs as vec2f) as float
  declare function distancesquared(byref rhs as vec2f) as float
  declare function dot(byref rhs as vec2f) as float
  declare function dotnormal(byref rhs as vec2f) as float
  declare function project(byref rhs as vec2f) as vec2f
  
end type

''  unary arithmetic non-member operator declarations
declare operator - (byref rhs as vec2f) as vec2f

''  binary arithmetic non-member operator declarations
declare operator + (byval lhs as vec2f, byref rhs as vec2f) as vec2f
declare operator - (byval lhs as vec2f, byref rhs as vec2f) as vec2f
declare operator * (byval lhs as vec2f, byref rhs as vec2f) as vec2f
declare operator * (byval lhs as float, byref rhs as vec2f) as vec2f
declare operator * (byval lhs as vec2f, byval rhs as float) as vec2f
declare operator / (byval lhs as vec2f, byval rhs as float) as vec2f

''  non-member function declarations
declare function absolute (byref lhs as vec2f) as vec2f
declare function normal (byref lhs as vec2f) as vec2f
declare function normalised (byref lhs as vec2f) as vec2f
declare function magnitude (byref lhs as vec2f) as float
declare function magnitudesquared (byref lhs as vec2f) as float
declare function distance (byval lhs as vec2f, byref rhs as vec2f) as float
declare function distancesquared (byval lhs as vec2f, byref rhs as vec2f) as float
declare function dot (byval lhs as vec2f, byref rhs as vec2f) as float
declare function dotnormal (byval lhs as vec2f, byref rhs as vec2f) as float
declare function project (byval lhs as vec2f, byref rhs as vec2f) as vec2f

''  constructors
constructor vec2f(): this.x = 0.0: this.y = 0.0: end constructor
constructor vec2f(byval x as float, byval y as float): this.x = x: this.y = y: end constructor

''  compound arithmetic member operators
operator vec2f.+= (byref rhs as vec2f): this.x += rhs.x: this.y += rhs.y: end operator
operator vec2f.-= (byref rhs as vec2f): this.x -= rhs.x: this.y -= rhs.y: end operator
operator vec2f.*= (byref rhs as vec2f): this.x *= rhs.x: this.y *= rhs.y: end operator
operator vec2f.*= (byref rhs as float): this.x *= rhs: this.y *= rhs: end operator
operator vec2f./= (byref rhs as float): this.x /= rhs: this.y /= rhs: end operator
operator vec2f.let (byref rhs as vec2f): this.x = rhs.x: this.y = rhs.y: end operator

''  member functions
function vec2f.absolute() as vec2f: return vec2f(abs(x), abs(y)): end function
function vec2f.normal() as vec2f: return vec2f(y, -x): end function 
function vec2f.normalised() as vec2f: return this/magnitude(): end function
function vec2f.magnitude() as float: return sqr(magnitudesquared()): end function
function vec2f.magnitudesquared() as float: return (x*x+y*y): end function
function vec2f.distance(byref rhs as vec2f) as float: return sqr(distancesquared(rhs)): end function
function vec2f.distancesquared(byref rhs as vec2f) as float: return (x-rhs.x)*(x-rhs.x)+(y-rhs.y)*(y-rhs.y): end function
function vec2f.dot(byref rhs as vec2f) as float: return (x*rhs.x+y*rhs.y): end function
function vec2f.dotnormal(byref rhs as vec2f) as float: return dot(rhs.normal()): end function
function vec2f.project(byref rhs as vec2f) as vec2f: return (dot(rhs)/magnitudesquared())*rhs: end function

''  unary arithmetic non-member operators
operator - (byref rhs as vec2f) as vec2f: return vec2f(-rhs.x, -rhs.y): end operator

''  binary arithmetic non-member operators
operator + (byval lhs as vec2f, byref rhs as vec2f) as vec2f: return vec2f(lhs.x+rhs.x, lhs.y+rhs.y): end operator
operator - (byval lhs as vec2f, byref rhs as vec2f) as vec2f: return vec2f(lhs.x-rhs.x, lhs.y-rhs.y): end operator
operator * (byval lhs as vec2f, byref rhs as vec2f) as vec2f: return vec2f(lhs.x*rhs.x, lhs.y*rhs.y): end operator
operator * (byval lhs as float, byref rhs as vec2f) as vec2f: return vec2f(lhs*rhs.x, lhs*rhs.y): end operator
operator * (byval lhs as vec2f, byval rhs as float) as vec2f: return vec2f(lhs.x*rhs, lhs.y*rhs): end operator
operator / (byval lhs as vec2f, byval rhs as float) as vec2f: return vec2f(lhs.x/rhs, lhs.y/rhs): end operator

''  non-member functions
function absolute (byref lhs as vec2f) as vec2f: return lhs.absolute(): end function
function normal (byref lhs as vec2f) as vec2f: return lhs.normal(): end function
function normalised (byref lhs as vec2f) as vec2f: return lhs.normalised(): end function
function magnitude (byref lhs as vec2f) as float: return lhs.magnitude(): end function
function magnitudesquared (byref lhs as vec2f) as float: return lhs.magnitudesquared(): end function
function distance (byval lhs as vec2f, byref rhs as vec2f) as float: return lhs.distance(rhs): end function
function distancesquared (byval lhs as vec2f, byref rhs as vec2f) as float: return lhs.distancesquared(rhs): end function
function dot (byval lhs as vec2f, byref rhs as vec2f) as float: return lhs.dot(rhs): end function
function dotnormal (byval lhs as vec2f, byref rhs as vec2f) as float: return lhs.dotnormal(rhs): end function
function project (byval lhs as vec2f, byref rhs as vec2f) as vec2f: return lhs.project(rhs): end function

main code

Code: Select all

''	symplectic integration algorithms
''	based on David Whysong's method:
''	http://www.aoc.nrao.edu/~dwhysong/prog/symp.cpp

#include once "vec2f.bi"

const Pi = 4*atn(1)				''	pi
const dt = 0.01						''	timestep, delta time
const G = 10							''	gravitational constant
const NumMasses = 10				''	number of planets
const ScrnWid		= 700			''	screen width
const ScrnHgt   = 700			''	screen height
const ScrnBrd   = 100			''	screen border

type MassType
	as single Mass, GM, Radius, RadiusSquared
	as vec2f Frc, Acc, Vel, Psn
end type

type SI4_Type
	
	public:
	
	declare sub Integrate()
	declare sub GetPosition(byval i as integer)
	declare sub GetVelocity(byval i as integer)
	
	private:
	
	const as double b 	= 2^(1/3)
	const as double a 	= 2-b
	const as double x0 = -b/a
	const as double x1 = 1/a
	as double d4(2) => { x1, x0, x1 }
	as double c4(3) => { x1/2, (x0+x1)/2, (x0+x1)/2, x1/2 }

end type

type SI6_Type
	
	public:
	
	declare sub Integrate()
	declare sub GetPosition(byval i as integer)
	declare sub GetVelocity(byval i as integer)
	
	private:
	
	const as double w1 = -0.117767998417887E1
	const as double w2 = 0.235573213359357E0
	const as double w3 = 0.784513610477560E0
	const as double w0 = (1-2*(w1+w2+w3))
	as double d6(6) => { w3, w2, w1, w0, w1, w2, w3 }
	as double c6(7) => { w3/2, (w3+w2)/2, (w2+w1)/2, (w1+w0)/2, (w1+w0)/2, (w2+w1)/2, (w3+w2)/2, w3/2 }

end type

type SI8_Type
	
	public:
	
	declare sub Integrate()
	declare sub GetPosition(byval i as integer)
	declare sub GetVelocity(byval i as integer)
	
	private:
	
	const as double W1 = -0.161582374150097E1
	const as double W2 = -0.244699182370524E1
	const as double W3 = -0.716989419708120E-2
	const as double W4 =  0.244002732616735E1
	const as double W5 =  0.157739928123617E0
	const as double W6 =  0.182020630970714E1
	const as double W7 =  0.104242620869991E1
	const as double W0 = (1-2*(W1+W2+W3+W4+W5+W6+W7))
	as double d8(14) => { W7, W6, W5, W4, W3, W2, W1, W0, W1, W2, W3, W4, W5, W6, W7 }
	as double c8(15) => { W7/2, (W7+W6)/2, (W6+W5)/2, (W5+W4)/2, (W4+W3)/2, (W3+W2)/2, _
											 (W2+W1)/2, (W1+W0)/2, (W1+W0)/2, (W2+W1)/2, (W3+W2)/2, (W4+W3)/2, _
											 (W5+W4)/2, (W6+W5)/2, (W7+W6)/2,  W7/2 }

end type

dim shared as MassType ptr Mass
dim shared as SI4_Type SI4
dim shared as SI6_Type SI6
dim shared as SI8_Type SI8

Mass = new MassType[NumMasses]

randomize timer

with Mass[0]
	.Mass = 1000
	.GM = G*.Mass
	.Radius = 12
	.RadiusSquared = .Radius*.Radius
	.Psn.X = ScrnWid\2
	.Psn.Y = ScrnHgt\2
end with

for i as integer = 1 to NumMasses-1
	with Mass[i]
		.Mass = 2 + rnd * 8
		.GM = G*.Mass
		.Radius = 4
		.RadiusSquared = .Radius*.Radius
		dim as single spawn_angle = rnd*2*pi
		dim as single spawn_distance = 50 + rnd * 250
		.Psn.X = Mass[0].Psn.X + spawn_distance * cos(spawn_angle)
		.Psn.Y = Mass[0].Psn.Y + spawn_distance * sin(spawn_angle)
		.Vel.X = Mass[0].Vel.X + sqr(Mass[0].GM/spawn_distance) * sin(spawn_angle)
		.Vel.Y = Mass[0].Vel.Y + sqr(Mass[0].GM/spawn_distance) * -cos(spawn_angle)
		
		Mass[0].vel.x -= .vel.x * (.GM/Mass[0].GM)
		Mass[0].vel.y -= .vel.y * (.GM/Mass[0].GM)
		
	end with
next

screeninfo ScrnWid, ScrnHgt
screenres ScrnWid, ScrnHgt, 16, 1

''	main loop
do
	
	SI4.Integrate
	'SI6.Integrate
	'SI8.Integrate
	
	screenlock
		cls
		for i as integer = 0 to NumMasses-1
			with Mass[i]
				circle (.Psn.X, .Psn.y), .Radius, rgb(64, 64, 255),,, 1, F
			end with
		next
	screenunlock
	
	sleep 1, 1
	
loop until multikey(1)

delete[] Mass

sub accelerate()
	
	for i as integer = 0 to NumMasses-1
		Mass[i].Acc = vec2f(0, 0)
	next
	
	for i as integer = 0 to NumMasses-2
		for j as integer = i+1 to NumMasses-1
			dim as vec2f DistVector  = Mass[i].Psn-Mass[j].Psn
			dim as single DistSquared = MagnitudeSquared(DistVector)
			if DistSquared > Mass[i].RadiusSquared+Mass[j].RadiusSquared then
				dim as single DistCubed = DistSquared*sqr(DistSquared)
				Mass[i].Acc -= Mass[j].GM*DistVector/DistCubed
				Mass[j].Acc += Mass[i].GM*DistVector/DistCubed
			end if
		next
	next

end sub

sub si4_type.integrate()

	for i as integer = 0 to 2
		GetPosition(i)
		accelerate()
		GetVelocity(i)
	next
	GetPosition(3)

end sub

sub si4_type.getposition(byval i as integer)

	dim as double Timestep = dt*c4(i)
	
	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Psn += .Vel*Timestep
		end with
	next

end sub

sub si4_type.getvelocity(byval i as integer)

	dim as double Timestep = dt*d4(i)

	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Vel += .Acc*Timestep
		end with
	next

end sub

sub si6_type.integrate()

	for i as integer = 0 to 6
		GetPosition(i)
		accelerate()
		GetVelocity(i)
	next
	GetPosition(7)

end sub

sub si6_type.getposition(byval i as integer)

	dim as double Timestep = dt*c6(i)
	
	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Psn += .Vel*Timestep
		end with
	next

end sub

sub si6_type.getvelocity(byval i as integer)

	dim as double Timestep = dt*d6(i)

	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Vel += .Acc*Timestep
		end with
	next

end sub

sub si8_type.integrate()

	for i as integer = 0 to 14
		GetPosition(i)
		accelerate()
		GetVelocity(i)
	next
	GetPosition(15)

end sub

sub si8_type.getposition(byval i as integer)

	dim as double Timestep = dt*c8(i)
	
	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Psn += .Vel*Timestep
		end with
	next

end sub

sub si8_type.getvelocity(byval i as integer)

	dim as double Timestep = dt*d8(i)

	for i as integer = 0 to NumMasses-1
		with Mass[i]
			.Vel += .Acc*Timestep
		end with
	next

end sub

Cheers,
Mike

[Frank Dodd]

I must say that I think you have done an excellent job with the presentation on this project. Your demo trailer was damn good fun to watch :D

Have you thought about pitching this at schools it would make a fun introduction in a science class I bet

[Lachie Dazdarian]

Yeah, the trailer is awesome. Where's that music from?

[ejc.cryptography]

Thanks for the comprehensive feedback, h4tt3n. I will be sure to incorporate some of the improvements you have suggested; I may contact you via email sometime soon if that's okay.

Yeah, the trailer is awesome. Where's that music from?

It is part of the score from the Blender Open Source movie 'Big Buck Bunny' (http://www.bigbuckbunny.org/index.php/c ... -download/). You can download the full score from the website for free.

[h4tt3n]

I may contact you via email sometime soon if that's okay.

Sure thing!

[Lachie Dazdarian]

I'm so happy for you two.