Celestial simulation

Hey all,

first, let me apologize in advance for the length of this post, and thank anyone willing to read and reply to it.

For quite some time now I've been wanting to make a program that simulates the movement of various celestial objects, mostly stars, as a result of gravity. I made such a program quite some time ago, using plain C and Allegro; stars had mass, density, radius, position and velocity, and were represented as little circles on the screen. When they got too close together, I would "turn off" gravity. Simple-as-a-pimple.

One day I though I would make a more advanced version, so I started to read up on star collisions. I quickly learned that simulating a collision between two stars is something one can not do on an average home computer with an average knowledge of C++. So I came up with a workaround; the star would be made of particles. That meant that the computer would think of a star as a star when nothing was in close proximity, but when two stars got close enough to each other, the stars would be treated as lots of particles. I think a crude illustration is in order:

Illustration #1

or something like that anyway. If this ever works out, I would like to introduce pressure into the whole program, specifically to treat the stars as a bunch of particles only if the gravitational stress would be enough to "rip the star apart".

Anyway, I decided that a stars mass would be generated, followed by it's density. From that, the radius would arise. And the number of particles the star was made up would be determined by the number of pixels it would occupy when drawn to the screen. In other words, one pixel=one particle. Here I encountered the first big problem. When a circle is drawn to the screen, it's made up of a bunch of squares. It's not a perfect circle. What I needed to find out was how to calculate the number of pixels set when drawing a circle with a radius of n pixels to the screen. And I had no idea how to do that. So I came here, asked my question, and was answered with the Midpoint Circle Algorithm. Unfortunately, I couldn't really figure out the math behind it, but I though I had grasped the theory, so I started working on my own, drawing up circles on squared paper, and eventually came up with the formula 2*r^2 + 6*r - 19, where r is the radius. It seemed to work with all the circles I had drawn up that had a radius bigger than three. So, yay, I had overcome the big obstacle. I chose SDL as the library I would use, and started learning to use it. Now, here comes the reason why I'm writing this article anyway. This whole thing depends on the fact that when the stars are far away from each other, I tell SDL to render a circle, and when the stars are close together, I tell SDL to render lots of pixels. But a few days ago I learned that, as opposed to Allegro, SDL does not have functions that handle drawing primitive shapes to the screen. Which means I was sc***ed. So I had a look at Google, and came up with something called SDL_Draw, which is an "add on" to SDL containing functions to draw primitive shapes to the screen. Problem was, I couldn't get it to work (I'm using Code::Blocks, and I'm terrible at making libraries run). As if it wasn't enough, I opened Microsoft Paint, made a couple of circles, counted the squares in them and compared them to what my formula gave me, just out of curiosity. And guess what; my formula was off. Which means that either Microsoft Paint uses some other algorithm that I'm not aware of, or I understood the algorithm wrong and drew up my circles wrong. So, here come the questions:

Should I stick with SDL, or use another library? If so, which one? If not, can someone help me get that SDL_Draw think working?
Does anyone know how to calculate the number of pixels that will be set when a circle of a radius r is drawn to the screen using SDL?

There are lots of other question I'd like to ask, but I'll post those in another topic, since they are irrelevant if I can't get this to work.

Cheers,

Gabe

To draw a tesselated sphere:

Angles
x = cos(alpha) * sin(beta) * radius
y = sin(alpha) * cos(beta) * radius
z = cos(beta) * radius

Alpha ranges 0 to 360.
Beta ranges 0 to 180 (or -90 to +90 ...I forget)

Angle increment for alpha: number_of_stacks / 360
Angle increment for beta: number_of_slices / 360

Gravitational force between two objects:

F = (m1 * m2 * G) / (d * d)

Where G is 6.67 * 10 ^ - 11 n/m^2
m1 and m2 are the masses of the objects expressed in Newtons
d is the distance between the objects expressed in meters

You will also need some type of integrator to make this all come together. If you are going to alter a body's position according to the gravitational formula then you will not need to use a hierarchical model for the system otherwise you will.

If you are planning on making each star be composed of particles equidistant around a center...then don't. This is going to kill any type of framerate. You should model your stars as meshes and then if you want to deform them just deform the mesh and/or split it when they collide. There are mesh splitting and deformation algos all over the net. This is not a simple task but it will probably be a better solution that using particles to represent entire objects.

Quote:

---

Originally Posted by Bubba

To draw a tesselated sphere:

Angles
x = cos(alpha) * sin(beta) * radius
y = sin(alpha) * cos(beta) * radius
z = cos(beta) * radius

Alpha ranges 0 to 360.
Beta ranges 0 to 180 (or -90 to +90 ...I forget)

Angle increment for alpha: number_of_stacks / 360
Angle increment for beta: number_of_slices / 360

---

Erm...I'm talking 2D here...and have no idea what that means anyway :D

Quote:

---

Originally Posted by Bubba

Gravitational force between two objects:

F = (m1 * m2 * G) / (d * d)

Where G is 6.67 * 10 ^ - 11 n/m^2
m1 and m2 are the masses of the objects expressed in Newtons
d is the distance between the objects expressed in meters

---

I know all this, I said I've already made the program. I just want to make it better.

Quote:

---

Originally Posted by Bubba

You will also need some type of integrator to make this all come together. If you are going to alter a body's position according to the gravitational formula then you will not need to use a hierarchical model for the system otherwise you will.

---

I can't do calculus, so to be completely honest I have again no idea what you're talking about. The way I did it was calculate F, get the acceleration using f=m/a, get velocity using v=a [delta]t, where [delta]t is a fixed constant. After that I calculate the displacement using s=v [delta]t.

Quote:

---

Originally Posted by Bubba

If you are planning on making each star be composed of particles equidistant around a center...then don't. This is going to kill any type of framerate. You should model your stars as meshes and then if you want to deform them just deform the mesh and/or split it when they collide. There are mesh splitting and deformation algos all over the net. This is not a simple task but it will probably be a better solution that using particles to represent entire objects.

---

Again, as I said, this is 2D..and why will it kill framerate? I know that computationally when they collide it will be a lot more difficult, but when they aren't colliding it shoudn't be a problem.

Please understand, I'm still in high school, so calculus isn't going to tell me much. And anyway, not that I want to sound rude, but I was asking about something else...

Cheers,

Gabe

Quote:

---

I can't do calculus, so to be completely honest I have again no idea what you're talking about. The way I did it was calculate F, get the acceleration using f=m/a, get velocity using v=a [delta]t, where [delta]t is a fixed constant. After that I calculate the displacement using s=v [delta]t.

---

This is correct except in a simulation you have to integrate the forces over time. High school math falls short in that they assume you are working with constants. Unfortunately in simulations you don't have constant time and you have variable acceleration over time which is normally your frame delta and/or your accumulated time due to the time delta. Please understand...what you are trying to do requires more than high school math.
Also, there was no mention of 2D or 3D.

Quote:

---

Anyway, I decided that a stars mass would be generated, followed by it's density. From that, the radius would arise. And the number of particles the star was made up would be determined by the number of pixels it would occupy when drawn to the screen. In other words, one pixel=one particle. Here I encountered the first big problem. When a circle is drawn to the screen, it's made up of a bunch of squares. It's not a perfect circle. What I needed to find out was how to calculate the number of pixels set when drawing a circle with a radius of n pixels to the screen. And I had no idea how to do that.

---

Radius arises out of density? Ok.
One pixel = one particle = a ton of particles for just 1 object. Not probable.

Quote:

---

When a circle is drawn to the screen, it's made up of a bunch of squares. It's not a perfect circle.

---

Because you are using particles instead of primitives. If you follow the equation of a circle which I assume you know since you seem to know everything now anyways then your circle will be a circle. Again I would show you how to create a 'sphere or circle' of particles using the equation but you probably already know how to do it. r^2 = x^2 + y^2. That is all you need.

Quote:

---

What I needed to find out was how to calculate the number of pixels set when drawing a circle with a radius of n pixels to the screen. And I had no idea how to do that.

---

You just said that radius is determined by density. This is absurd but hey it is your simulation and not mine. If radius is determined by density and you know density then you know the radius. I'm lost here.

If you say your max particle density of any one star is 10000 pixels then the number of pixels for any star in your system with density expressed as a normalized float is:

num_particles = density * max_particles

When density is 1.0f you get max particles. When density is 0.5f you get half of max_particles and so on. You just need to decide on a theoretical maximum of particles. If you are going to render these particles like this:

Code:

---

for (unsigned int i = 0;i < m_NumParticles; ++i)
{
  //Draw particle
}

---

This would be a very bad idea. A better idea would be to transform a portion of the particles and then draw. After that you are always drawing what you just transformed. This way you don't get graphics pipeline stalls. If you just transform and then render the entire buffer of particles then the video card is locked out at several points thus causing stalls.

And, BTW, density has nothing to do with the size of an object. You can take a teaspoon of material from a white dwarf and it will have a mass of a million newtons and you can take a teaspoon of material from our own sun and it will have a very small mass.

A huge boulder of granite hass less mass than a small boulder of solid steel.

Quote:

---

And anyway, not that I want to sound rude, but I was asking about something else...

---

Not to be rude but I don't think you know exactly what you are asking.

If you want to render your 'stars' as empty circles then you can use the circle drawing algo. The angle increment for a circle is:
angle_incr = PI / number_of_facets

Where the number of facets is the desired approximation of a circle that you want. If you want to draw perfect circles you will have to use another algo but I do not recommend it. Circles with 32 facets are usually a good approximation.

Okay; first, I'd appreciate if you'd leave your remarks out of this. If you think I'm incompetent, fine, but please, keep it to yourself. I asked a polite question, you politely answered. However, you didn't answer to the questions I asked, which I politely pointed out (if for any reason it sounded aggressive or offensive, I apologize) in my answer.
Second, 90% of what you say there is based on the complete misunderstanding of what I wrote.

Now I know you can't get an exact simulation using only what I know, but I'm not doing it for NASA; I'm just doing it for fun. I don't need a perfect pinpoint position of every single object at any given time. All I want is to just have fun. As for the 2D/3D, I though it would be obvious from my constant references to SDL (not to be confused with OpenGL), which is, the last time I checked, a 2D graphics library. Anyway, yes, I'm talking about 2D.
You get the radius from the mass and density using this equation.
As for the whole particle thing, read about the midpoint circle algorithm. It will explain better what I mean when I say "not a perfect circle". And I would only be rendering it particle by particle if a collision was in process. Otherwise it would just be

Code:

---

makecircle(x, y, r)

---

Now, to re-post my questions:

Quote:

---

Should I stick with SDL, or use another library? If so, which one? If not, can someone help me get that SDL_Draw think working?
Does anyone know how to calculate the number of pixels that will be set when a circle of a radius r is drawn to the screen using SDL?

---

Hope this clears things up.

Cheers,

Gabe

If you know the radius of the circle in pixels then you can calculate the amount of pixels in a square of the same size, and use simple trigonometry to shave off the unecessary pixels. This would not be expensive because you would only need to calculate this once.

Also, I see absolutely wrong with what bubba is saying, and he is not wrong. He may not be answering your questions, but he is telling what you are doing wrong. It is worthless to ignore constructive criticism such as he has provided.

SDL is a slow library when it comes transformations, distortions, and other CPU intensive graphics. HGE is a DirectX based engine that works pretty well. If you use HGE you will want to recompile the engine, and change a few settings in the code. Like set it up to use Hardware Vertex processing, and a few other things.

How would I "shave off" the pixels? A quick example/explanation will suffice. And I'm not ignoring the constructive part of what bubba said, I'm ignoring the part that was based on misunderstanding my post. But thanks to both of you ;-)

Cheers,


Gabe

Quote:

---

Originally Posted by G4B3

How would I "shave off" the pixels?
Gabe

---

Illustration

Code:

---

void PutPixel(SDL_Surface *Surf_Dest, int X, int Y, Uint32 Pixel)
{
        if(Surf_Dest == NULL || X < 0 || X >= Surf_Dest->w || Y < 0 || Y >= Surf_Dest->h)
                return;
       
        int bpp = Surf_Dest -> format -> BytesPerPixel;
        Uint8 *p = (Uint8*) Surf_Dest -> pixels + Y * Surf_Dest -> pitch + X * bpp;
       
        switch(bpp)
        {
                case 1:
                        *p = Pixel;
                        break;
                case 2:
                        *(Uint16*) p = Pixel;
                        break;
                case 3:
                        if(SDL_BYTEORDER == SDL_BIG_ENDIAN)
                        {
                                p[0] = (Pixel >> 16) & 0xFF;
                                p[1] = (Pixel >> 8) & 0xFF;
                                p[2] = (Pixel & 0xFF);
                        }
                        else
                        {
                                p[0] = (Pixel & 0xFF);
                                p[1] = (Pixel >> 8) & 0xFF;
                                p[2] = (Pixel >> 16) & 0xFF;
                        }
                        break;
                case 4:
                        *(Uint32*) p = Pixel;
                        break;
        }
}

SDL_Surface *DrawCircle(int Radius, Uint32 Color)
{
        const SDL_VideoInfo *Info = SDL_GetVideoInfo();
        if(Info == NULL)
                return NULL;
       
        SDL_Surface *Surf_Circle = SDL_CreateRGBSurface(SDL_HWSURFACE, Radius << 1, Radius << 1, Info -> vfmt -> BitsPerPixel, Info -> vfmt -> Rmask, Info -> vfmt -> Gmask, Info -> vfmt -> Bmask, Info -> vfmt -> Amask);
       
        if(Surf_Circle == NULL)
                return NULL;
       
        Uint32 Transparent = SDL_MapRGBA(Surf_Circle -> format, 0, 0, 0, 0);
        if(SDL_SetColorKey(Surf_Circle, SDL_SRCCOLORKEY, Transparent) < 0)
                return NULL;
       
        if(SDL_MUSTLOCK(Surf_Circle))
                SDL_LockSurface(Surf_Circle);
       
        for(int Y = -Radius;Y < Radius;Y++)
        {
                int half_row_width = int(sqrt(double((Radius * Radius) - (Y * Y))));
                for(int X = -half_row_width;X < half_row_width;X++)
                        PutPixel(Surf_Circle, Radius + X, Radius + Y, Color);
        }
       
        if(SDL_MUSTLOCK(Surf_Circle))
                SDL_UnlockSurface(Surf_Circle);
       
        return Surf_Circle;
}

---

With the code above, just count how many times PutPixel is called by DrawCircle, and you will know how many pixels has been drawn by SDL, it's inflexible and very slow though, but it's just like you want, hope this helps :D

wow, thanks guys :D Okay, this gets me somewhere.

Cheers,

gabe

Quote:

---

As for the whole particle thing, read about the midpoint circle algorithm. It will explain better what I mean when I say "not a perfect circle". And I would only be rendering it particle by particle if a collision was in process.

---

I understand the midpoint algorithm. I used to use it way back in the DOS days prior to using Direct3D. What I am recommending is to forget about the circles and approximate them with a series of line segments composed of points equidistant from a center...IE: an approximated circle. As to what library you should use is completely up to you and comes down to simple API calls that can be done in any library. However if you don't understand the concept of what you are trying to do you could use any API and still be completely off base.

If I were to continue to allow you to use your own flawed algorithms and data structures to represent objects then I would be helping you fail. Your particle idea is at best confusing to me. Are you representing the stars as a group of particles equidistant from a central point and then transforming the entire system to world space or are you just wanting to use circles? Particles threw me way off here. In 2D or 3D you do not want to use perfect circles. They are slow, cumbersome, not easily deformable, and so on. Half the battle of programming graphics is figuring out which way NOT to do it rather than exactly how to do it.


Hierarchichal Solar system representation
There are multiple ways to represent a solar system. The easiest is a hierarchy.
Your star system could easily be represented as a hierarchy even in 2D. If you assume the center of the star system resides at 0,0 in local space then all planets reside at offsets from 0,0.

Code:

---

struct SystemNode
{
  Matrix4x4 toRoot;
  Matrix4x4 localOrientation;
  Vector3 pos;
  SystemNode *pParent;
  SystemNode *pChild;
};

---

At load you would walk up the hierarchy from a leaf node or one with no child. You would concatenate all the matrices up the hierarchy and store that in toRoot. This means that the toRoot matrix will get you from the current node to the root transform. After this you multiply the resulting matrix by the model matrix of the planet.

Since you are doing a solar system your hierarchy will usually just be one node deep. But if you have any moons orbiting planets you would have two or more nodes deep. So the ending hierarchy might look like this for our solar system:

Root - Sun
- Mercury - SystemNode->pChild = 0; SystemNode->pParent = Sun
- Venus - SystemNode->pChild = 0;SystemNode->pParent = Sun
- Earth - SystemNode->pChild = Moon;SystemNode->pParent = Sun
- Moon - SystemNode->pChild = 0;SystemNode->pParent = Earth

The toRoot transform for Moon:

Code:

---

Moon->toRoot = Earth->localOrientation * Sun->localOrientation

---

The transform needed to correctly render Moon:

Code:

---

Matrix4x4 objectMatrix = Moon->localOrientation * Moon->toRoot * SolarSystem->modelMatrix;

---

Without the toRoot functionality you would do this every frame:

Code:

---

Matrix4x4 matObject;

while (pNode->pParent)
{
  matObject *= pNode->localOrientation;
  pNode = pNode->pParent;
}

matObject *= SolarSystem->modelMatrix();

---

Where the model matrix of the solar system represents the current rotation, translation, and
scaling of the entire solar system.

True gravity representation
Another way is to represent all of the objects as just that...unrelated objects that are grouped together simply because of gravitational attraction. Keep in mind this is very expensive. Since you want to effect the velocity vector of an object based on the formula for gravitational attraction between two objects you will have to compute every object's influence on one object.

So if you want to affect the moon's course due to gravity:

Code:

---


Vector3 vecVelocity;

for (int i = 0;i < numObjects; ++i)
{
  Vector3 vecGravity = computeAttraction(this,Objects[i]);

  //Grab speed from vector
  float speed = vecVelocity.length();

  //Normalize the attraction vector
  vecGravity.Normalize();

  //Alter our velocity vector based on gravitational vector
  vecVelocity += vecGravity;
}

Vector3 computeAttraction(const SolarBody &Body1,const SolarBody &Body2)
{
  //Compute gravity vector
  Vector3 vecAttraction = Body2.Pos - Body1.Pos;
 
  //Get distance squared
  float distance = toSolarBody2.length();
  distance *= distance;

  //Normalize vector
  toSolarBody2.normalize();
 
  //extend gravitation vector towards body based on formula
  return toSolarBody2 * (Body1.Mass * Body2.Mass * gravity_constant) / (distance);
}

---

Once you finally compute this then you can place it into your existing F = ma integrator with F being the force you just calculated.

Code:

---

Acceleration = Force / mass;
newVelocity = Velocity + Acceleration * delta;
speed = newVelocity.length();
newPos = Pos + newVelocity * delta;

Pos = newPos;
Velocity = newVelocity;

---

Keep in mind this is ignoring any angular momentum and mass moments of inertia. It also does not account for drag but since this is space you shouldn't have any.

Bob would be more suited to the task. He could definitely chime in. This is very basic but that is what you said you were after.
If you just wanted to know which graphics API to use...you took a long time in your post getting there. Also what you posted concerning F = ma is your integrator so you probably already have this done and working.