6
\$\begingroup\$

After a relative hiatus of around two years, I've decided to start delving back into the programming world; to start, I've appropriated one of my older projects, CLIGL, and have made a solar system simulator that runs exclusively within a console interface.

There are several classes of importance, whose purposes and functions should be highlighted before the presentation of the code itself:

  • Body: this class is responsible for storing the rendering and physics properties of a gravitational body, performing numerical integration of the position and velocity properties (via IntegratePosition and IntegrateVelocity), and rendering the body itself.
  • BodySystem: this class is responsible for maintaining a list of gravitational bodies and simulating them properly; it will also render each body by calling their Render methods. Additionally, it will also render a set of physics properties for each simulated body in the top left corner of the console window.
  • EntryPoint: fairly self-explanatory; this class contains the initialization for and the main loop of the simulation, as well as a few constants. It also contains two helper methods for generating a decorative field of stars in the background.

EntryPoint.cs

using System;
using System.Diagnostics;
using System.Collections.Generic;
using CLIGL;

namespace GravitySimulator
{
    /// <summary>
    /// Class containing the program entry point, several constants/statics, and several
    /// helper functions which have no other place.
    /// </summary>
    public class EntryPoint
    {
        public const int WINDOW_WIDTH = 159;
        public const int WINDOW_HEIGHT = 60;
        public static (bool, char)[] STAR_FIELD = new (bool, char)[WINDOW_WIDTH * WINDOW_HEIGHT];

        /// <summary>
        /// Randomly populate the star field with boolean values; false indicates no star, whereas
        /// true indicates that a star is present. A character value will also be randomly generated.
        /// </summary>
        /// <param name="chanceOfStar">The chance that a star will be generated, from 0-100.</param>
        /// <param name="starTypes">The different types of stars.</param>
        public static void PopulateStarField(float chanceOfStar, char[] starTypes)
        {
            Random randomGenerator = new Random();
            for(int i = 0; i < WINDOW_WIDTH * WINDOW_HEIGHT; i++)
            {
                STAR_FIELD[i] = (
                    randomGenerator.Next(0, 100) <= chanceOfStar, 
                    starTypes[randomGenerator.Next(0, starTypes.Length)]
                );
            }
        }

        /// <summary>
        /// Render the previously-generated random star field.
        /// </summary>
        /// <param name="buffer">The buffer to which the star field is to be rendered.</param>
        public static void RenderStarField(ref RenderingBuffer buffer)
        {
            for(int i = 0; i < WINDOW_WIDTH * WINDOW_HEIGHT; i++)
            {
                if(STAR_FIELD[i].Item1)
                {
                    buffer.SetPixel(i, new RenderingPixel(STAR_FIELD[i].Item2, ConsoleColor.DarkGray, ConsoleColor.Black));
                }
            }
        }

        /// <summary>
        /// Application entry point.
        /// </summary>
        /// <param name="args">Command line arguments.</param>
        public static void Main(string[] args)
        {
            RenderingWindow window = new RenderingWindow("Gravity Simulator", WINDOW_WIDTH, WINDOW_HEIGHT);
            RenderingBuffer buffer = new RenderingBuffer(WINDOW_WIDTH, WINDOW_HEIGHT);

            Stopwatch timeAccumulator = new Stopwatch();
            timeAccumulator.Start();
            float oldTime = (float)timeAccumulator.Elapsed.TotalSeconds;
            float newTime;
            float dt;

            PopulateStarField(1, new char[] { '.', ',', '`', '\'' });
            BodySystem system = new BodySystem(
                new List<Body>() {
                    new Body("Star", '@', ConsoleColor.Yellow, ConsoleColor.Red, 10000, 0, 0, 0, 0),
                    new Body("Planet 1", '@', ConsoleColor.Gray, ConsoleColor.DarkGray, 10, 20, 0, 0, 35),
                    new Body("Planet 2", '@', ConsoleColor.DarkGray, ConsoleColor.Gray, 10, -30, 0, 0, -25),
                    new Body("Planet 3", '@', ConsoleColor.White, ConsoleColor.Gray, 10, -50, 0, 0, -15)
                }
            );

            while(true)
            {
                newTime = (float)timeAccumulator.Elapsed.TotalSeconds;
                dt = newTime - oldTime;

                buffer.ClearPixelBuffer(RenderingPixel.EmptyPixel);
                RenderStarField(ref buffer);

                buffer.SetString(0, 0, $"DT :: {dt}", ConsoleColor.White, ConsoleColor.Black);
                system.RenderBodies(ref buffer);
                window.Render(buffer);

                system.SimulateBodies(dt);
                oldTime = newTime;
            }
        }
    }
}

Body.cs

using System;
using CLIGL;

namespace GravitySimulator
{
    /// <summary>
    /// Contains position, velocity, and rendering information pertaining to a gravitational body,
    /// simulated by the System class.
    /// </summary>
    public class Body
    {
        public const float G = 0.1f;

        public string Name { get; set; }
        public char Character { get; set; }
        public ConsoleColor ForegroundColor { get; set; }
        public ConsoleColor BackgroundColor { get; set; }

        public float Mass { get; set; }
        public float X { get; set; }
        public float Y { get; set; }
        public float Vx { get; set; }
        public float Vy { get; set; }
        public float Ax { get; set; }
        public float Ay { get; set; }

        /// <summary>
        /// Constructor for the Body class.
        /// </summary>
        /// <param name="name">The name of the body.</param>
        /// <param name="character">The rendering character for the body.</param>
        /// <param name="foregroundColor">The foreground color of the body.</param>
        /// <param name="backgroundColor">The background color of the body.</param>
        /// <param name="mass">The mass of the body.</param>
        /// <param name="x">The initial X position of the body.</param>
        /// <param name="y">The initial Y position of the body.</param>
        /// <param name="vx">The initial X velocity of the body.</param>
        /// <param name="vy">The initial Y velocity of the body.</param>
        public Body(
            string name, 
            char character, 
            ConsoleColor foregroundColor, 
            ConsoleColor backgroundColor, 
            float mass, 
            float x, 
            float y, 
            float vx, 
            float vy
        )
        {
            this.Name = name;
            this.Character = character;
            this.ForegroundColor = foregroundColor;
            this.BackgroundColor = backgroundColor;
            this.Mass = mass;
            this.X = x;
            this.Y = y;
            this.Vx = vx;
            this.Vy = vy;
            this.Ax = 0.0f;
            this.Ay = 0.0f;
        }

        /// <summary>
        /// Integrate the position of the Body; this is done by adding the current velocities,
        /// Vx and Vy, to the current positions, X and Y.
        /// </summary>
        /// <param name="dt">The amount of time since the last frame, delta time.</param>
        public void IntegratePosition(float dt)
        {
            this.X += this.Vx * dt;
            this.Y += this.Vy * dt;
        }

        /// <summary>
        /// Integrate the velocity of the Body; this is done by calculating the attractive force
        /// between the current body and a provided body, then deriving the accelerations, and then
        /// adding the derived accelerations to the current velocities.
        /// </summary>
        /// <param name="otherBody">The body to which attraction is calculated.</param>
        /// <param name="dt">The amount of time since the last frame, delta time.</param>
        public void IntegrateVelocity(Body otherBody, float dt)
        {
            float distance = this.DistanceBetween(otherBody);
            float directionX = (otherBody.X - this.X) / distance;
            float directionY = (otherBody.Y - this.Y) / distance;

            float attractiveForce = (G * this.Mass * otherBody.Mass) / distance;
            float acceleration = attractiveForce / this.Mass;

            float aX = acceleration * directionX;
            float aY = acceleration * directionY;
            this.Ax = aX;
            this.Ay = aY;

            this.Vx += this.Ax * dt;
            this.Vy += this.Ay * dt;
        }

        /// <summary>
        /// Find the distance between the current body and another body.
        /// </summary>
        /// <param name="otherBody">The other body.</param>
        /// <returns>A distance value.</returns>
        public float DistanceBetween(Body otherBody)
        {
            return (float)Math.Sqrt(
                Math.Pow(otherBody.X - this.X, 2) +
                Math.Pow(otherBody.Y - this.Y, 2)
            );
        }

        /// <summary>
        /// Render the body to a provided buffer.
        /// </summary>
        /// <remarks>This will render the body relative to the center of the buffer, not the top left corner.</remarks>
        /// <param name="buffer">The buffer to which the body is rendered.</param>
        /// <param name="renderName">Whether or not to render the name of the body.</param>
        public void Render(ref RenderingBuffer buffer, bool renderName = false)
        {
            int roundedX = (int)Math.Round(buffer.BufferWidth / 2.0f + this.X);
            int roundedY = (int)Math.Round(buffer.BufferHeight / 2.0f + this.Y);
            buffer.SetPixel(roundedX, roundedY, new RenderingPixel(this.Character, this.ForegroundColor, this.BackgroundColor));

            if(renderName)
            {
                buffer.SetString(roundedX + 1, roundedY, this.Name, this.BackgroundColor, ConsoleColor.Black);
            }
        }
    }
}

BodySystem.cs

using System;
using System.Collections.Generic;
using CLIGL;

namespace GravitySimulator
{
    /// <summary>
    /// This class is responsible for simulating a system of bodies.
    /// </summary>
    public class BodySystem
    {
        public List<Body> Bodies { get; set; }

        /// <summary>
        /// Constructor for the BodySystem class.
        /// </summary>
        /// <param name="bodies">A list of bodies to simulate.</param>
        public BodySystem(List<Body> bodies)
        {
            this.Bodies = bodies;
        }

        /// <summary>
        /// Simulate the bodies.
        /// </summary>
        /// <param name="dt">The amount of time since the last frame, delta time.</param>
        public void SimulateBodies(float dt)
        {
            foreach(Body body in this.Bodies)
            {
                body.IntegratePosition(dt);
            }

            foreach(Body otherBody in this.Bodies)
            {
                foreach(Body body in this.Bodies)
                {
                    if(body.Name == otherBody.Name)
                    {
                        continue;
                    }

                    body.IntegrateVelocity(otherBody, dt);
                }
            }
        }

        /// <summary>
        /// Render all the simulated bodies.
        /// </summary>
        /// <param name="buffer">The buffer to which to render.</param>
        public void RenderBodies(ref RenderingBuffer buffer)
        {
            foreach(Body body in this.Bodies)
            {
                body.Render(ref buffer, true);
            }

            int offset = 1;
            foreach(Body body in this.Bodies)
            {
                buffer.SetString(
                    0, 
                    offset, 
                    $"{body.Name} :: " +
                    $"P(x,y) = {body.X.ToString("F1")},{body.Y.ToString("F1")} :: " +
                    $"V(x,y) = {body.Vx.ToString("F1")},{body.Vy.ToString("F1")} :: " +
                    $"A(x,y) = {body.Ax.ToString("F2")},{body.Ay.ToString("F2")}", 
                    ConsoleColor.White, 
                    ConsoleColor.Black
                );
                offset++;
            }
        }
    }
}

I am not as concerned with the correctness of my code prima facie (although suggestions for improvements in this regard would be appreciated nonetheless); rather, I am more concerned with the correctness of my implementation of the physics of gravitation, and any issues that may be present in it.

A few side notes:

  • This project links a previous project of mine, CLIGL; if you wish to test out this simulator, you will need to download and compile CLIGL to a .DLL and link it accordingly.
  • I highly recommend that, should you wish to test this project for yourself, that you set the font of your console window to a raster font with a size of 12x16.
  • For those who may not wish to go to the effort of downloading and compiling CLIGL and this project, I have uploaded two videos demonstrating the simulator: one without labels and one with labels.
  • Finally, I am well aware that the value of the gravitational constant \$G\$ is \$6.67408\times10^{-11}\$ and not \$0.5\$, as I have implemented it here; using the correct value of \$G\$ results in a simulation that runs far too slow to appreciate.
\$\endgroup\$
2
\$\begingroup\$

Gravitational Constant

In your description you mention the gravitational constant is 0.5 for performance reasons, in the code it is actually 0.1.

This Pointer

There are times when the this pointer is necessary, but it is not necessary in the code under review (EntryPoint.cs, Body.cs and BodySystem.cs).

Variable Names

To some extent this is personal preference, but the 2 letter variable names are generally not descriptive enough. I can see x and y as positions, but I don't necessarily understand dt (deltaTime), Vx, Vy (?velocity?), Ax, Ay (?acceleration?), G (?gravitational constant?). Variable names should be descriptive so that comments aren't necessary.

Comments

The param comments for name, mass, foregroundColor and backgroundColor aren't necessary, the other param comments would not be necessary if better variable names were used. Keep in mind the more comments there are the more the text in the source code needs to change when maintenance is performed. By using self documenting code and only commenting important algorithm information, edits to comments are kept to a minimum.

Private Versus Public

There is at least one function in Body (float DistanceBetween(Body otherBody)) that could be private rather than public.

\$\endgroup\$
  • 1
    \$\begingroup\$ With regards to the terse variable names, my rationale was simply to follow the existing conventions of physics notation (I was conflicted on the issue of their brevity though, for what it's worth). Regarding comments, I 100% agree; I've always had somewhat of an obsessive tendency to over-comment. I very much appreciate the answer! \$\endgroup\$ – Ethan Bierlein Nov 4 '19 at 20:30
  • \$\begingroup\$ In the context of physics it's ok to have short variable names, as long as these variables are understandable without any additional context. Ideally there should be a comment near those variable names that explains the names. Linking to the Wikipedia article about gravity should be enough. \$\endgroup\$ – Roland Illig Feb 13 at 23:17
6
\$\begingroup\$

Finally, I am well aware that the value of the gravitational constant 𝐺 is 6.67408×10−11 and not 0.5, as I have implemented it here; using the correct value of 𝐺 results in a simulation that runs far too slow to appreciate.

It is strange to me that you chose to change the gravitational constant rather than simply choosing a different time step; can you say more about why this seemed like the right thing to do? (Bonus points for justifications via references to ST:TNG.)

I am more concerned with the correctness of my implementation of the physics of gravitation, and any issues that may be present in it.

You are right to be concerned. Your approach is where everyone starts, but it is not accurate.

First off, the obvious way in which it is inaccurate is that you're doing Newtonian dynamics but that does not reflect reality. The orbit of Mercury was noticeably "wrong", for example, until Einstein determined that identical clocks run at a different rate when orbiting closer to the Sun than the Earth. But let's assume that you do not wish to model factors due to general relativity and stick to the world as it was known to Newton.

So let's just consider your method here, which as I mentioned is the naive way of numerically simulating n bodies:

  • Start with known positions and velocities for all bodies
  • Compute forces from positions
  • Compute accelerations from forces
  • Positions update based on velocities
  • Velocities update based on accelerations
  • Step forward one time unit
  • Repeat

FYI this is a variation of Euler's Method.

This is a good first approximation, but it has many problems starting with it is not conservative. That is, it is possible that the system has more energy or less energy after each step than the previous step, which is not physically plausible. But the real problem is that if you ever get into a situation where energy is being consistently added over several subsequent steps, you can end up with a runaway scenario where a body gets violently flung away from the system at high speed.

Learning how to adapt your method to a more physically accurate simulation will take you very deep indeed into the study of numerical methods, so I recommend getting a good book on the subject. But here are some thoughts to get you going:

  • You have the same time step for every particle, and I assume, the same timestep for each step. This is potentially very bad. The faster a particle is moving (or accelerating), the smaller its timestep should be to help avoid unphysicalities showing up in the simulation. Slow the timestep down when particles are moving fast, or even better, have every particle keep track of its own timestep which gets shorter the faster it goes. Coming up with a good data structure to do the computations in the right order when every particle has a changing time step is a great exercise.

  • You can be more clever about computing the change in position and velocity than simply getting the acceleration at a point in time. Suppose for example you have the position and velocity of a particle at time t. You compute the acceleration based on the state at t and from that you compute the increment to the state that gets you a new state at t + dt. That's what you're doing now. But here's the trick. Now suppose you compute the equivalent increment by doing two steps of dt/2. That increment will be different, and both will be wrong, but if you are clever about it then you can average the increments and be much closer to the accurate solution. This is the Runga-Kutta method; I wrote an RK4-5 solver to solve differential equations when I was in university and it is quite straightforward code to write.

That will get you started I hope. If you are interested in going deeper down this rabbit hole, this is a good overview of where to go next:

http://www.scholarpedia.org/article/N-body_simulations


One other thing I just thought of. I assume that in your graphics system you have the ability to position the camera at an arbitrary position and orientation. It would be an interesting test of your system to model three bodies: the Sun, Jupiter, and a Trojan near but not at L4, then put the camera focus on Jupiter and watch the Trojan "orbit" the stable L4 point.

Alternatively, put the camera on the sun and rotate the camera at the same rate as Jupiter moves; that would look like this: http://sajri.astronomy.cz/asteroidgroups/hitrfix.gif

\$\endgroup\$
  • \$\begingroup\$ Op I hope you who Eric is... I would look at what he says. \$\endgroup\$ – Gilad Feb 13 at 22:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.