Samstag, 18. Mai 2013

Planet Gravitation Map Simulator

Screenshot CUDA / OpenGL Demo of Gravitation Map Simulator
This is just another fun demo for CUDA and OpenGL. It is a slightly modified version of the so called "random oscillating magnetic pendulum" (ROMP). A pixel represents a start position of one particle. This object becomes attracted by all the planets around due to their gravity. After a while the object is going to hit a planet. Now that corresponding pixel gets the color of that planet. You see an evolving map of "gravity" structures in realtime (more or less *cough*).


Download: (tested on Linux and Windows (VS2010), needs CUDA, OpenGL, GLUT, GLEW)

Download Source Code v1.1
Download Source Code v1.0



Positions and masses of the planets can be changed by the user as well as the scale of the map.

The simulation uses Euler or Runge-Kutta integration method for solving the differential equation. It is also possible to use doubles instead of floats for better precision.

The idea for this demo comes from a friend, he also implemented the algorithm in processing. I ported this to CUDA, so the stuff is entirely computed and rendered on the GPU. You need a Nvidia CUDA capable GPU for that.
Screenshot 2.1, different coloring
Screenshot 2.2

Computation

Given a particle at position $ \vec{y}(t)=\left(x(t), y(t)\right)$ at a time t. This particle has a (positive) mass m. In our space there are n planets with their fixed positions $ \vec p_1,\ldots,\vec p_n$ and (positive) masses $ m_1,\ldots,m_n$.
The following equation represents the gravitational force of planet i acting on our particle. This force comes from Newton's law of universal gravitation:

$\displaystyle \vec F_i = G\frac{m_i\cdot m}{r_i^2}\vec e_{r_i}$
G is just a gravitation constant and can be neglected for our purposes. r is the distance between the particle and the planet, $ \vec e_r$ is the force direction vector:
$\displaystyle \vec e_r = \frac{\vec p_i - \vec y}{\Vert \vec p_i - \vec y\Vert}$
Hence, there is the following force at a time t:
$\displaystyle \sum_{i=1}^n\vec F_i(t) = \sum_{i=1}^n m\cdot \frac{m_i\cdot\left(\vec p_i - \vec y(t)\right)}{\Vert\vec p_i-\vec y(t)\Vert^3}$
Newton's second law gives the equality:

$\displaystyle \vec F=m\cdot\vec a\Rightarrow \sum_{i=1}^n\vec F_i(t) = m\cdot \frac{\mathrm d^2 \vec y(t)}{\mathrm dt^2}$
The mass m of the particle can be eliminated in the equation. This ordinary differential equation system of order 2 need to be integrated two times to obtain the position function:

\begin{displaymath}\begin{split}\vec v(t) =& \frac{\mathrm dy(t)}{\mathrm dt}=\i...
...c y(t) =& \int\limits_0^t \vec v(\tau)\mathrm d\tau \end{split}\end{displaymath}    

This integration can be solved numerically, e.g. by Euler's method:

\begin{displaymath}\begin{split}\vec v_{k+1} =& \vec v_k + h\cdot\sum_{i=1}^n\fr...
...y_k\Vert^3}\\ \vec y_{k+1} =& \vec y_k + h\cdot v_k \end{split}\end{displaymath}    

with $ y_k{=}y(t_k),\, v_k{=}(t_k),$ and $ y(t_0){=}y_0,\, v(t_0){=}v_0$ as the initial values. h is the size of every step and v represents the velocity of the particle.


Here you see some videos, which are showing the simulation more or less in realtime.


More:
Article by Ingo Berg, 2006, on magnetic pendulum with implementation (using Beeman's algorithm):
http://www.codeproject.com/Articles/16166/The-magnetic-pendulum-fractal

Mathematica implementation:
http://nylander.wordpress.com/2007/10/27/magnetic-pendulum-strange-attractor/

Some more stuff, videos:
http://magnetmfa.wikispaces.com/pendula?responseToken=05ae0f7708a5c4c989455051dc970f570
http://www.youtube.com/watch?v=duy8s8C7-Uc
http://www.youtube.com/watch?v=QXf95_EKS6E