Solar System Simulation Series: Part 1 Part 2 Part 3

Like almost everything in this universe, nothing is perfect, including the planetary orbits of our solar system. Though some are more circular than others, all of them are ellipses/ovals, instead of perfect circles. My program parametrically plots ellipses with the intent of approximating the planetary orbits using their semi-major axis and approximating the semi-minor axis with the perihelion (closest approach to the sun).

For those of you that have read my previous post, I decided to abandon the derivative and vector based approach since this new method is much less computationally demanding.

An ellipse (or sometimes circle) can be defined parametrically by the following:

[latex]x = a+bcos(t)[/latex]

and

[latex]y = c+dsin(t)[/latex]

where a and c are the translation factors (i.e. if a is 30, the ellipse’s center has an x-coordinate of 30) and b and d are the semi major and semi minor axes of the ellipse, with the smaller one being the semi-minor axes. If b=d it becomes a circle – though since planetary orbits are ellipses, we won’t have any of these.

As such I took these equations and codified them in C Sharp, inserting the values of the planet’s semi-major axis for b and the approximated semi-minor axis for d. Then I translated the ellipse being drawn to the center of the Picturebox, where I put the sun.

Example with only the first four planets: Mercury, Venus, Earth, and Mars

To draw the orbits in real time, I created a thread specifically (hence the multiple “Invoke” commands below) for processing the elliptical orbits in the Form Load event, which then calls a function named simulator: It runs a do loop that is paused by a checkbox, and also calls two other methods: orbit, to calculate the orbits, and drawpoint, which actually draws the orbits:

public void simulation() { do { mercurypos = orbit(mercurymass); venuspos = orbit(venusmass); earthpos = orbit(earthmass); marspos = orbit(marsmass); jupiterpos = orbit(jupitermass); saturnpos = orbit(saturnmass); uranuspos = orbit(uranusmass); neptunepos = orbit(neptunemass); coords[0] = solpos; coords[1] = mercurypos; coords[2] = venuspos; coords[3] = earthpos; coords[4] = marspos; coords[5] = jupiterpos; coords[6] = saturnpos; coords[7] = uranuspos; coords[8] = neptunepos; for (int x = 1; x <= 8; x++) { drawPoint(coords[x], prevcoords[x], penarray[x]); } pictureBox1.Invoke((MethodInvoker)delegate { pictureBox1.Image = bmp; }); } while (checkBox1.Checked); }

The method orbit basically calculates the next position of each of the planets per iteration. Originally, the planets’ angular rate of change were equivalent, which is unrealistic. So, I put in a way for the values for orbital velocity to be calculated with respect to Mercury’s – for example, for Earth, I multiplied the default orbital velocity by the ratio of the real world average orbital velocity of Mercury over Earth (29.78 / 47.87 km/s) as well as the ratio of Mercury’s orbital circumference to Earth’s. This yielded reasonably accurate orbital velocities for all of the planets based on real data (Jupiter takes 12 Earth years to orbit, Neptune takes 165 etc). The variable iterations keeps track of “time” and is used in drawPoint. divider simply scales the values to the 1280×720 resolution picturebox.

public double[] orbit(double planetmass) { double divider; if (radioButton3.Checked) divider = timesfull; else divider = timesjup; double track4 = 0; trackBar4.Invoke((MethodInvoker)delegate { track4 = trackBar4.Value; }); if(radioButton3.Checked) iterations = iterations + (0.006 * track4 / 30 *5/12); else iterations = iterations + (0.006*track4/30); label16.Invoke((MethodInvoker)delegate { label16.Text = Convert.ToInt32(iterations).ToString(); }); double mercurycircum = 2 * Math.PI * 57.9; double venuscircum = 2 * Math.PI * 108.2; double earthcircum = 2 * Math.PI * 149.6; double marscircum = 2 * Math.PI * 227.9; double jupitercircum = 2 * Math.PI * 778.3; double saturncircum = 2 * Math.PI * 1427; double uranuscircum = 2 * Math.PI * 2871; double neptunecircum = 2 * Math.PI * 4497.1; if (planetmass == mercurymass) { double x = 640 + (46001200 *Math.Cos(iterations) / divider); double y = 360 + (57909050 *Math.Sin(iterations) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == venusmass) { double x = 640 + (107477000 * Math.Cos(iterations * (35.02 / 47.87) * (mercurycircum / venuscircum)) / divider); double y = 360 + (108208000 * Math.Sin(iterations * (35.02 / 47.87) * (mercurycircum / venuscircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == earthmass) { double x = 640 + (147095000 * Math.Cos(iterations * (29.78 / 47.87) * (mercurycircum / earthcircum)) / divider); double y = 360 + (149598023 * Math.Sin(iterations * (29.78 / 47.87) * (mercurycircum / earthcircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == marsmass) { double x = 640 + (206700000 * Math.Cos(iterations * (24.077 / 47.87) * (mercurycircum / marscircum)) / divider); double y = 360 + (227939200 * Math.Sin(iterations * (24.077 / 47.87) * (mercurycircum / marscircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == jupitermass) { double x = 640 + (740550000 * Math.Cos(iterations * (13.07 / 47.87) * (mercurycircum / jupitercircum)) / divider); double y = 360 + (778299000 * Math.Sin(iterations * (13.07 / 47.87) * (mercurycircum / jupitercircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == saturnmass) { double x = 640 + (1350000000 * Math.Cos(iterations * (9.69 / 47.87) * (mercurycircum / saturncircum)) / divider); double y = 360 + (1429390000 * Math.Sin(iterations * (9.69 / 47.87) * (mercurycircum / saturncircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == uranusmass) { double x = 640 + (2742000000 * Math.Cos(iterations * (6.81 / 47.87) * (mercurycircum / uranuscircum)) / divider); double y = 360 + (2875040000 * Math.Sin(iterations * (6.81 / 47.87) * (mercurycircum / uranuscircum)) / divider); double[] pos = { y, x }; return pos; } else if (planetmass == neptunemass) { double x = 640 + (4460000000 * Math.Cos(iterations * (5.43 / 47.87) * (mercurycircum / neptunecircum)) / divider); double y = 360 + (4504450000 * Math.Sin(iterations * (5.43 / 47.87) * (mercurycircum / neptunecircum)) / divider); double[] pos = { y, x }; return pos; } else { double[] sol = { 0, 0 }; return sol; } }

After telling the computer to crunch the numbers for the orbit, I told it to draw the numbers, scaled to a 1280×720 resolution, picturebox of the orbits. At first, after completing one orbit, you couldn’t tell where the planet was since it was simply redrawing the orbit in the same color, so I came up with two cosine functions to cycle through different colors. If you’ve read my previous article about detecting colors, you’ll have a sense of what RGB is. I basically told the computer to check to see if a color’s R, G, and B values were less than 99 or greater than 155, then I added a variable value defined by func and funcneg, which have ranges of 0 to 100 and -100 to 0 respectively. The previously mentioned variable iterations affects what func and funcneg’s output are. Finally, the method draws a 3×3 square surrounding the position of the planet with the given color determined by both func, funcneg, and the originally defined color. The Try Catch mechanism is here so that if a drawn orbit is outside the range of the picturebox, it still exists but is just not drawn (will be useful if asteroids are added later on).

public void drawPoint(double[] planetpos, double[] prevplanetpos, Pen pen) { Point currPoint = new Point(Convert.ToInt32(planetpos[1]), Convert.ToInt32(planetpos[0])); double func = 0; double funcneg = 0; trackBar4.Invoke((MethodInvoker)delegate { func = -50 * Math.Cos(iterations / 100 * Math.Log(750 - trackBar4.Value)) + 50; funcneg = 50 * Math.Cos(iterations / 100 * Math.Log(750 - trackBar4.Value)) - 50; }); Color c = pen.Color; int r = c.R; int g = c.G; int b = c.B; if (c.R < 155) r = c.R + Convert.ToInt32(func); else if(c.R > 99) r = c.R + Convert.ToInt32(funcneg); if (c.G < 155) g = c.G + Convert.ToInt32(func); else if (c.G > 99) g = c.G + Convert.ToInt32(funcneg); if (c.B < 155) b = c.B + Convert.ToInt32(func); else if (c.G > 99) b = c.B + Convert.ToInt32(funcneg); Color d = Color.FromArgb(r, g, b); try { for (int x = -1; x <= 1; x++) { for (int y = -1; y <= 1; y++) { bmp.SetPixel(currPoint.X + x, currPoint.Y + y, d); } } } catch { } }

This is the final result, with multiple options already operational – switching from a zoomed in view of the inner solar system (see above) or the full planetary view (below), changing the time multiplier between 1 day/second to 730 (2 years) per second, and the checkbox that pauses and resumes the simulation.

To do list – stuff not implemented yet:

- 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 [latex]F=\frac{Gm_{asteroid}m_{planetorsun}}{r^2}[/latex] for each major body in the solar system. The next position of the asteroid can be implied by the sum of these vectors. Edit:done click for part 2.
- 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.

If you have any ideas to improve the simulation, feel free to leave a comment below!

Edit – Form Load Code and Global Declared Variables:

public static double earthmass = 5.9722 * Math.Pow(10,24); public static double solmass = 333000 * earthmass; public static double mercurymass = 0.055 * earthmass; public static double venusmass = 0.815 * earthmass; public static double marsmass = 0.107 * earthmass; public static double jupitermass = 317.8 * earthmass; public static double saturnmass = 95.159 * earthmass; public static double uranusmass = 14.536 * earthmass; public static double neptunemass = 17.147 * earthmass; public static double gravconstant = 6.674 * Math.Pow(10,-11); public Pen[] penarray = new Pen[9]; public static Pen earthpen = new Pen(System.Drawing.Color.Aquamarine, 3); public static Pen mercurypen = new Pen(System.Drawing.Color.Brown, 3); public static Pen venuspen = new Pen(System.Drawing.Color.Coral, 3); public static Pen marspen = new Pen(System.Drawing.Color.Red, 3); public static Pen jupiterpen = new Pen(System.Drawing.Color.Orange, 3); public static Pen saturnpen = new Pen(System.Drawing.Color.PaleGoldenrod, 3); public static Pen uranuspen = new Pen(System.Drawing.Color.DarkSeaGreen, 3); public static Pen neptunepen = new Pen(System.Drawing.Color.DarkBlue, 3); public static Pen solpen = new Pen(System.Drawing.Color.Yellow, 3); public static double[] screenjupsize = { 1.4507377777778 * Math.Pow(10, 9), 816.04 * Math.Pow(10, 6) }; public static double[] screenfullsize = { 8071111111 , 4540000000}; double timesjup = 1000000; double timesfull = 13000000; public double[][] coords = new double[9][]; public static double[] earthpos = { 149598023, 0 }; public static double[] solpos = { 0, 0 }; public static double[] mercurypos = { 0 , 57909050 }; public static double[] venuspos = { 108208000, 0 }; public static double[] marspos = { 227939200, 0 }; public static double[] jupiterpos = { 778299000, 0 }; public static double[] saturnpos = { 1429390000, 0 }; public static double[] uranuspos = { 2875040000, 0 }; public static double[] neptunepos = { 4504450000, 0 }; public double[][] prevcoords = new double[9][]; public static double[] prevearthpos = { 149598023, 0 }; public static double[] prevsolpos = { 0, 0 }; public static double[] prevmercurypos = { 57909050, 0 }; public static double[] prevvenuspos = { 108208000, 0 }; public static double[] prevmarspos = { 227939200, 0 }; public static double[] prevjupiterpos = { 778299000, 0 }; public static double[] prevsaturnpos = { 1429390000, 0 }; public static double[] prevuranuspos = { 2875040000, 0 }; public static double[] prevneptunepos = { 4504450000, 0 }; public Thread drawer; public Thread asteroids; double xsin = Math.Sin(0); double ycos = Math.Cos(0); Graphics g; Bitmap bmp; double iterations = 0; private void Form10_Load(object sender, EventArgs e) { bmp = new Bitmap(pictureBox1.Width, pictureBox1.Height); g = Graphics.FromHwnd(pictureBox1.Handle); coords[0] = solpos; coords[1] = mercurypos; coords[2] = venuspos; coords[3] = earthpos; coords[4] = marspos; coords[5] = jupiterpos; coords[6] = saturnpos; coords[7] = uranuspos; coords[8] = neptunepos; prevcoords[0] = prevsolpos; prevcoords[1] = prevmercurypos; prevcoords[2] = prevvenuspos; prevcoords[3] = prevearthpos; prevcoords[4] = prevmarspos; prevcoords[5] = prevjupiterpos; prevcoords[6] = prevsaturnpos; prevcoords[7] = prevuranuspos; prevcoords[8] = prevneptunepos; penarray[0] = solpen; penarray[1] = mercurypen; penarray[2] = venuspen; penarray[3] = earthpen; penarray[4] = marspen; penarray[5] = jupiterpen; penarray[6] = saturnpen; penarray[7] = uranuspen; penarray[8] = neptunepen; bmp = new Bitmap(pictureBox1.Width, pictureBox1.Height); bmp.SetPixel(640,360,Color.Yellow); bmp.SetPixel(641, 361, Color.Yellow); bmp.SetPixel(641, 360, Color.Yellow); bmp.SetPixel(641, 359, Color.Yellow); bmp.SetPixel(640, 361, Color.Yellow); bmp.SetPixel(640, 359, Color.Yellow); bmp.SetPixel(639, 359, Color.Yellow); bmp.SetPixel(639, 360, Color.Yellow); bmp.SetPixel(639, 361, Color.Yellow); pictureBox1.Image = bmp; drawer = new Thread(simulation); drawer.Priority = ThreadPriority.AboveNormal; drawer.IsBackground = true; drawer.Start(); asteroids = new Thread(asteroid); asteroids.Priority = ThreadPriority.Normal; asteroids.IsBackground = true; asteroids.Start(); }