Solar System Simulation Part 2 Asteroids Code

Posted on

Solar System Simulation Series: Part 1 Part 2 Part 3

As was promised in the previous post, an explanation of how my program’s logic worked was upcoming. Without further ado:

First, I needed code to actually generate the asteroids. I also wanted each asteroid to have a randomly generated color and position on the 1280×720 resolution picturebox on the form.  This was achieved by declaring public lists of a class called GravAsteroids (with a namespace “Gravity”) and of Colors. Then, in the form load event, a for loop creates new “GravAsteroids” and assigns both a randomized color and position to each instance of the class:

public List<Gravity.GravAsteroids> asteroids = new List<Gravity.GravAsteroids>();
public List<Color> asteroidcolor = new List<Color>();
double numasteroids = trackBar1.Value;

for (int x = 0; x <= numasteroids - 1; x++)
{
Random s = new Random(Guid.NewGuid().GetHashCode());
Random t = new Random(Guid.NewGuid().GetHashCode());
Random randomGen = new Random(Guid.NewGuid().GetHashCode());
KnownColor[] names = (KnownColor[])Enum.GetValues(typeof(KnownColor));
KnownColor randomColorName = names[randomGen.Next(names.Length)];
Color randomColor = Color.FromKnownColor(randomColorName);
asteroids.Add(new Gravity.GravAsteroids(this, s.Next(1,720), t.Next(1,1280)));
}


Then, I needed the program to find the force vector (in my case comprised of an angle and magnitude) acting on any single asteroid – this is the purpose of the GravAsteroids class and Gravity namespace. The first part was accomplished with the atan2 function under the built in Math.h library to find the angle between each and every pair of two asteroids/planets/sun. Atan2 is defined as the following:

$\text{atan2}(y,x) = \begin{cases} \arctan(\frac y x) &\text{if } x > 0, \\ \arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y > 0, \\ \arctan(\frac y x) – \pi &\text{if } x < 0 \text{ and } y < 0, \\ +\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0,\\ -\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\ \text{undefined} &\text{if } x = 0 \text{ and } y = 0. \end{cases}$

,where x and y signify the coordinate values of a point/vector, and results in an output of the angle in radians between the positive portion of the x-axis to the point/vector. Because I was using doubles and real distances in kilometers before drawing to the screen, the last possibility is extremely unlikely, and when it does happen, my code deletes and then regenerates that particular instance of GravAsteroids.

In the code, I found the difference between the coordinates/position vector of each asteroid and every other asteroid and planets/sun to yield the vector difference between each pair of positions and inserted this into the atan2 function. This essentially translates one of the points/the difference vector to the origin (0,0) for calculation purposes.
Note slight syntactical differences depending on if calculating for asteroid or planet/sun:

//if finding angle between two asteroids, "other" referring to other asteroid.
double Angle = Math.Atan2((other.Y - Y), (other.X - X));
//if finding angle between asteroid and established astronomical bodies - sun and 8 planets
double Angle = Math.Atan2((planetcoor[1] - Y), (planetcoor[0] - X));


Next, I needed the magnitude of the force vector. Remember the equation for the gravitational attraction between two masses $ma = F = \frac{Gm_1*m_2}{r^2}$ from my first solar system post?
The magnitude of the gravitational force vector is an adaptation of that equation in code. The *1000 multiplier is to account for the trackbar controlling the asteroids’ masses having a real range of 1-100, I multiplied all instances of the asteroidmass by 1000 to more accurately simulate a medium-low mass asteroid. The last part of the equation converts the number from force in kiloNewtons to a more program friendly (unnamed) force unit consisting of pixels as written in the following equation: $\frac{kg*km}{s^2} = \frac{kg*pixels}{timesfull*s^2}$ timesfull represents the number of number of kilometers is represented by each side of a square pixel (approximately 13 million km when zoomed out to Neptune’s orbit). It changes if the user zooms into the inner solar system.
Again, note the slight syntactical differences depending on if calculating for asteroid or planet/sun:

//for two asteroids
double Magnitude = (gravconstant * Form.asteroidmass * 1000 * other.Form.asteroidmass * 1000 / Math.Sqrt(((Math.Pow((other.X - X), 2)) + (Math.Pow((other.Y - Y), 2))))) / Form.timesfull;
//for asteroid to planet/sun attraction
double Magnitude = (gravconstant * Form.asteroidmass * 1000 * Form.massarray[which] / Math.Sqrt(((Math.Pow((planetcoor[0] - X), 2)) + (Math.Pow((planetcoor[1] - Y), 2)))))/Form.timesfull;


To put the angle and magnitude into a vector recognizable by C Sharp:
Complex is under System.Numerics, and Vector is under System.Windows.

Complex c = Complex.FromPolarCoordinates(Magnitude, Angle);
Vector r = new Vector(c.Real, c.Imaginary);


To convert from force in $\frac{kg*pixels}{timesfull*s^2}$ to velocity in $\frac{pixels}{iteration}$ – pixels per iteration, where iterations are the number of processing cycles that my program has accomplished – a sort of “timer” for my program.

For this, I needed to convert from the iterations that currently was my program’s “time” to real SI units of years/days/hours/second etc. I noticed that Earth completes an orbit on my program every 250 iterations, so with a little math, there are approximately 250/365/24/3600 ≈ 0.000007927 iterations per second or 1/0.000007927 ≈ 126144 seconds per iteration.

So therefore, I need to multiply by seconds per iteration and divide by the mass of the two objects in question combined (timesfull cancelled out the km): $\frac{kg*pixels}{s^2} * \frac{s}{iter*kg} = \frac{pixels}{s*iteration}$. I can’t cancel out the last seconds in the denominator, so I just make it a user adjustable value to speed up or slow down the simulation. The following method “Update” is called in every do loop iteration, the details of which you can find in my previous post. The color c refers to the randomized color generated in the above code. Grav and planetinter refer to methods that return the force vectors.

Note that my code also includes a feature that counts the number of impacts that each major solar system body (sun + 8 planets) have sustained since the program started running, by checking the diameter of each of the bodies and if the asteroid has come close enough to impact. It then regenerates the asteroid at a random location with zero velocity as to prevent an undefined solution to the aforementioned atan2 function. This “respawn” of the asteroid also happens when an asteroid goes out of bounds of the 1280×720 picturebox.

Finally, after calculating velocity, it uses the drawPoint method in the main form (also detailed in my previous post) to draw the current position of the asteroid on the bitmap that is displayed on the picturebox of the main form.

The method “Update” is the following:

public void Update(Color c)
{
Vector Force = new Vector(0, 0);
foreach (GravAsteroids asteroid in Form.asteroids)
{
if (this != asteroid)
Force += Grav(asteroid);
}
foreach (double[] planetsun in Form.coords)
{
int counter = 0;
if (X != planetsun[0] && Y != planetsun[1])
{
Force += planetinter(planetsun, counter);
}
else
{
Random s = new Random(Guid.NewGuid().GetHashCode());
int hit = s.Next(1,Convert.ToInt32(Form.timesfull));
if(hit<=Form.radius[counter]*2) { Form.labelarray[counter].Invoke((System.Windows.Forms.MethodInvoker)delegate { Form.labelarray[counter].Text = (Convert.ToInt32(Form.labelarray[counter].Text) + 1).ToString(); }); Random u = new Random(Guid.NewGuid().GetHashCode()); Random t = new Random(Guid.NewGuid().GetHashCode()); X = u.Next(1, 1279); Y = t.Next(1, 719); Velocity.X = 0; Velocity.Y = 0; } else Force += planetinter(planetsun, counter); } counter++; } double TimeElapsedSinceLastUpdate = (Form.iterations*126144); Vector acceleration = Force / ((Form.asteroidmass * 1000 * Form.asteroids.Count)*Form.massarray.Sum); Velocity += acceleration * TimeElapsedSinceLastUpdate; X = (X + Velocity.X * VelocityScale); Y = (Y + Velocity.Y * VelocityScale); if (X > 1279 || X <= 1 || Y> 1279 || Y<= 1)
{
Random s = new Random(Guid.NewGuid().GetHashCode());
Random t = new Random(Guid.NewGuid().GetHashCode());
X = s.Next(1, 1279);
Y = t.Next(1, 719);
Velocity.X = 0;
Velocity.Y = 0;
}
double[] xy = new double [] { X, Y };
Pen pen = new Pen(c, 1);
Form.drawPoint(xy, pen, 2);
}


Here is another screenshot of results after 50000 iterations, with 3 impacts into the Sun! (and fortunately no impacts on Earth, or any of the planets for that matter, but it’s just a matter of time!)

And here’s the video up to 2500 iterations:

To do list with more stuff crossed off!

• Asteroids that are implied from the options panel, affected by the gravity of all the planets and the sun, with separate force vectors from each planet and sun to the tune of $F=\frac{Gm_{asteroid}m_{planetorsun}}{r^2}$ for each major body in the solar system. The next position of the asteroid can be implied by the sum of these vectors.
• Increase in orbital velocity as when orbital radius decreases (planets move faster during their closest approach and slower when farther away), instead of constant velocity. Can probably be done with a sinusoidal function of some sort.
• Make the Sun one of the foci of the elliptical orbits a lá Kepler’s Laws of Planetary Motion, instead of the center of the orbit. Can probably be done by shifting the orbit over by the difference between the planet’s perihelion (closest approach) and aphelion (furthest distance).
• The Sun isn’t stationary relative the the planets either! It is affected by the gravity of the planets in our solar system and moves over time, though such movement is relatively tiny (a few solar diameters/decade). Fun fact, this stellar movement is the “wobble” that astronomers use to detect some of the 1-2 thousand or so exoplanets (planets outside the solar system) they’ve found so far.
• Would require a lot of work, for very little noticeable change – these perturbations won’t be seen on the 1280×720 resolution picturebox that I’ve used. But, for the sake of accuracy, this is down the pipeline, though at the bottom.
• New Goal: Gravitational anomaly in the bottom left – asteroids tend to be attracted to an area of empty space in the bottom left – a planetary level mass apparently exists there. This needs to be fixed.
• And a new goal: number of Earth days passed since the simulation started, this will be a little more complicated since now with the asteroid simulation, the number of iterations per second decreases as a function of the number of asteroids because of processing power limits.
• New Goal, add number of impacts of asteroids on the sun and planets

Expect an update that includes the second, third and fifth bullet point in Part 3!