powered by NetLogo

view/download model file: SolarSystem.nlogo

This model mimics the solar system. Setup creates the sun, five planets, and a "comet". Unlike in the real solar system,
these orbiting bodies are in the same plane, and setup distributes the planets randomly. However, the planets are at the
correct relative distances from the sun and have correct relative masses. The comet can be placed up to 100 AU (an AU is
the distance from earth to sun, or Astronomical Unit) from the sun and assigned an "ellipticity" which adjusts its initial
velocity to result in a orbit that can range from a circle to an extreme ellipse.

As the model is run, the ratio of the current distance from the sun to the initial distance can be plotted for earth and
the comet. If the gravitational effects of all planets, in addition to the sun, are computed, earth's orbit can be seen to
vary slightly in ellipticity. If only the sun's gravity is considered, earth's orbit is still elliptical, but does not vary.
The orbit of the comet can vary greatly, depending upon how close and at what speed it approaches other planets or the sun.
Depending upon the initial positions of the planets and comet, and the comet's ellipticity, the comet could orbit indefinitely,
immediately be ejected from the solar system, or be knocked around into other orbits before being ejected.

The third plot, "Comet - Closest Planet Distance" plots the distance from the comet to Jupiter and/or Saturn when the distance
between them is within the range set on the sliders "Comet-Jupiter" and "Comet-Saturn" in AU. Whenever the comet comes within
this range of either planet, if a plot for neither planet was already in progress, the previous plot is cleared and a new one
started. This can give a good view of how the comet is moving with respect to these planets, and how close it actually gets,
when within the specified ranges.

To Run the Model

Just press Setup to use the default speed, distance of comet to sun, and comet orbit ellipticity. Other settings can be
changed during the run. Setup resets Zoom to 15, so setting it to another value initially isn't useful.

Press Go

Notes:

If you want a different speed, initial comet position, and/or ellipticity, make these changes before pressing Setup. Changing
these later has no effect.

The speed adjustment actually changes the size of time intervals at which the new planet positions are updated. This affects
the accuracy of the orbits. At a speed setting of 10 there are about 1000 updates of the position of the earth per
revolution of the sun, but only about 250 for Mercury. At a speed setting of 1 there are ten times as many updates per
revolution. The effects of speed settings on the orbit calculations can be seen by setting the Sun_Only switch to "on" and
observing the orbital plots for both earth and comet at different speed settings. The ellipticity of the earth's orbit
decreases as the number of updates per year increases. At higher ellipticity settings for the comet, which bring it closer
to the sun, its orbit should be repetitive if the Sun_Only switch is on, but at higher speed settings it may vary and even
become unstable due to insufficient updates where the comet's orbital path is making sharp turns. Note that when the Sun_Only
switch is "off" so that gravitational effects of the planets as well as the sun on the comet are computed, the comet may
come close enough to any planet, or the sun, to cause sufficiently abrupt changes in orbital direction that the time interval
represented by higher speed settings may be too large to compute its orbit correctly. Also, no provision has been made to
note the "collision" of the comet with a planet or the sun, so the resulting acceleration will almost certainly fling the
comet out of the solar system.

Setting the comet's initial distance from the sun to be near that of Saturn or Jupiter makes it easy to see the effects of these
large planets on the comet's orbit. The farther out the comet starts, the longer it is likely to take to see an obvious
interaction with a planet.

The Earth - Sun Distance plot is scaled according to whether the Sun_Only switch is on and the current Speed setting, so that
the changes are easily visible. The scaling is based on trial and error observation.

While the Model is Running

To increase the speed of execution without affecting computation accuracy, turn the display off by pressing Display On/Off.
This is a toggle which alternately stops the display from showing the changing planet positions, or resumes showing their
positions. When the display is off, you can follow the progress by looking at the "earth years" monitor which tracks the
number of orbits earth has made, and by looking at the plots of the earth's and comet's orbits.

You can change the distance between the comet and Jupiter and the comet and Saturn, with sliders "Comet-Jupiter" and
"Comet-Saturn", which will trigger the following:

If the display is off and the comet approaches within the slider set distances of either Jupiter or Saturn,
the display is turned on and execution is stopped for two seconds to aid in starting to visually track the comet's
orbit from that point.

The "Comet - Closest Planet Distance" plot will show the distance between comet and Jupiter and/or comet and Saturn,
either, or both as long as the distance is less than the slider setting. If the distance to either planet comes
within range and a plot is not already in progress, the previous plot is cleared and a new one started. This is
independent of whether the display was on or off.

To follow a comet's orbit when it goes off-screen, lower the Zoom value. Otherwise, when any orbiting body passes off-
screen, it is hidden until it comes back on-screen. Also, if you increase the Zoom setting so that you view fewer planets,
say only through Mars, execution speed increases.

The Plotting On/Off switch is obvious.

The Labels On/Off switch, can be used to display a label for each object in the solar system. The model then runs much slower.

The Sun_Only switch can be used to observe how the planets' orbits behave when all gravitational interactions are computed, or
when only the sun's gravitational pull is computed. For example, the comet's orbit will not be affected by close encounters
with any of the planets when the Sun_Only switch is "on".

Calculations

Computing the acceleration of object #1 towards object #2 due to gravity:

From F = ma = -GMm/r^2 and solving for the x and y components of acceleration a we get accx = -GMx/r^3 and accy = -GMy/r^3
where accx and accy are the accelerations in the x and y directions, x and y are the distances in x and y directions between
the objects #1 and #2, M is the mass of object #2 (the mass m of object #1 cancels out), G is the gravitational constant, and
r = sqrt(x^2 + y^2). G is expressed in terms of Astronomical Units, AU and earth masses, Me, giving units of AU^3/(Me sec^2).

Computing the inital conditions

After the planets are randomly placed at their proper orbiting distances from the sun, the initial velocities are computed
from the equation v = -sqrt(G*M/r) where v = tangential velocity. The velocities in the x and y directions are then
x-velocity = v * sin (angle) and y-velocity = v * cos (angle) where angle is reported from (towards-nowrap sun - 180).

The initial accelerations are computed as described above, with each planet having the effects of the sun's and all the other
planets' gravity on itself computed.

At each time interval a new x and y position for each planet and the comet is computed based on the x and y velocities for
that time period. Since the velocities at the beginning of the time interval increase by (acceleration * time) the
average velocity for each time period is used. For the first time period this is computed by adding to the initial velocity
the acceleration * time/2.

Updating planets' and comet's positions

This is just a matter of repeating the computations above of new accelerations added to current velocities in x and y
directions.

PROCEDURES

breed [ planets ]
 
turtles-own [mass xvelocity yvelocity accx accy newx newy radius avg-radius radius-ratio oldy orbitcount mylabel]
globals [sun earth mars mercury venus jupiter saturn time scale G display_on NearPlanet dist2jupiter dist2saturn jupiterplot saturnplot]
 
to setup
 display
 ca
 clear-output
 set display_on true
 set NearPlanet false
 set Zoom 13
 set scale 1
 set scale Zoom * max-pxcor / 150
 set time Speed / 100
 set G 1.205 * 10 ^ (0 - 10) ;Gravitation constant in units of newton*m^2/kg^2 converted to units of Astronomical Units, 
 ;AU and earth masses, Me, giving units of AU^3/(Me sec^2)
 crt 1
 create-planets 7  ;the seventh planet is processed as turtle 7, and is an extra orbiting body, i.e. "comet",
 ; that might be affected by other planets, if it gets close enough
 ask planets [set color green]
 set sun turtle 0
 ask sun [
 set color yellow
 set size 5
 set mass 329390 
 set mylabel "Sun "
 ]
 ask turtles [set shape "circle"]
 set mercury turtle 1
 ask mercury [
 set mass .0549 
 set color brown
  set avg-radius 0.388 
 set mylabel "Mer "
 ]
 set venus turtle 2
 ask venus [
 set mass .8073 
 set color pink
 set avg-radius .722 
 set mylabel "V "
 ]
 set earth turtle 3
 ask earth [
  set mass 1 
 set avg-radius 1.00 
 set mylabel "E "
 set color blue
 ]
 set mars turtle 4
 ask mars [
 set mass .1065 
 set color red
 set avg-radius 1.53 
 set mylabel "Mar "
 ]
  set jupiter turtle 5
 ask jupiter [
 set mass 314.5 
 set avg-radius 5.2 
 set shape "jupiter"
 set color orange
 set mylabel "J "
 ]
 set saturn turtle 6
 ask saturn [
 set mass 94.07 
 set shape "saturn"
 set avg-radius 9.54 
 set color green
 set mylabel "Sat "
 set heading 45
 ]
 
 ask turtle 7 [
 set mass .00001 
 set avg-radius Comet_Start 
 set color white 
 set mylabel "Com "
 set newx avg-radius 
 set newy 0 ;sets comet at y distance of zero, and at distance from sun specified by the slider
 ]
 
 ;planets are randomly positioned, except for the comet, but at the correct distance from the sun, in Astromonical Units. Masses are in earth-masses.
 ;calculations are done with newx and newy as actual coordinates, then translated into NetLogo screen units scaled to show more or fewer planets orbiting, as desired
 ask planets [set size 3 + log mass 10 ]
 ask planets with [who < 7][
  ifelse random 2 = 0 
 [set newx random-float avg-radius]
 [set newx 0 - random-float avg-radius] 
 ifelse random 2 = 0 
 [set newy sqrt (avg-radius ^ 2 - newx ^ 2)]
  [set newy 0 - sqrt (avg-radius ^ 2 - newx ^ 2)]
 ]
 ask planets [
 ifelse abs (newx * scale) > max-pxcor or abs (newy * scale) > max-pxcor 
 [set hidden? true setxy newx * scale newy * scale]
  [set hidden? false setxy newx * scale newy * scale]
 ] 
 
 ask planets [set radius sqrt(newx ^ 2 + newy ^ 2) ]
 ifelse ViewLabels = false 
 [ask turtles [set label ""]]
 [ask planets [set label mylabel]]
 
 ;Compute the initial velocity in x and y directions to produce an approximate circular orbit of radius calculated above regardless of where start point is.
 ;This comes from a = v^2/r for tangential velocity and a = -G*M/r^2. Solve for v at x = r. Then v = -sqrt(G*M/r). Then compute v in x and y directions.
 ;Using sin ((towards-nowrap sun) - 180) and cos ((towards-nowrap sun) - 180) to compute y and x components of velocity, respectively, gives the correct 
 ;sign for the resulting value regardless of which quadrant the planets are located.
 
 ask planets with [who < 7][
 set yvelocity sqrt (([mass] of sun * G) / radius) * sin ((towards-nowrap sun) - 180)
 set xvelocity 0 - sqrt (([mass] of sun * G) / radius) * cos ((towards-nowrap sun) - 180)
 ]
 
 ask turtle 7 [
 set xvelocity 0
 ;;The higher the Comet_Ellipticity setting the lower the initial velocity relative to that needed for a circular orbit, so the orbit is more elliptical
 ;;Since the comet can be initially off-screen the towards-nowrap primitive cannot be used. Since it is always at -90 degrees with xvelocity = 0,
 ;;the value sin ((towards-nowrap sun) - 180) is always sin (-270).
 set yvelocity sqrt (([mass] of sun * G) / radius) * sin (-270) * ((100 - Comet_Ellipticity) / 100)
 ]
 
 ;compute acceleration in x and y directions due to sun's gravity only
 ask planets [
 set accx 0 - ((newx * [mass] of sun * G) / radius ^ 3) 
 set accy 0 - ((newy * [mass] of sun * G) / radius ^ 3)
 ]
 
 ;compute acceleration in x and y directions due to gravities of other planets - each planet computes influence of other planets excluding itself
 ;and combine accelerations due to sun and other planets
 if Sun_only = false ;;switch allows turning off effects of other planets and use gravity of sun only 
 [ask planets [
 let my-accx accx
 let my-accy accy
 ask planets with [who != [who] of myself] [
 set my-accx my-accx - (((mass * G ) * ([newx] of myself - newx))/(sqrt ( (newx - [newx] of myself ) ^ 2 + (newy - [newy] of myself) ^ 2)) ^ 3)
 set my-accy my-accy - (((mass * G ) * ([newy] of myself - newy))/(sqrt ( (newx - [newx] of myself ) ^ 2 + (newy - [newy] of myself) ^ 2)) ^ 3)
 ]
 set accx my-accx
 set accy my-accy
 ]
 ]
 ;compute the average velocity that will apply for the time period used to compute the next position, estimated by computing velocity at one half the time period
 ;and applying it for the entire time period. The new velocity is the initial velocity plus the acceleration applied over t/2 seconds.
 ask planets [
 set xvelocity xvelocity + accx * time / 2 
 set yvelocity yvelocity + accy * time / 2
 ]
 
 ; setup-plots
end
 
 
to go
 orbit-sun
 draw-orbits
 ifelse comet? [ ask turtle 7 [set size 2
 if Orbits = "All" [pendown] ] ]
  [ask turtle 7 [set size 0.01 ]]
end
 
to draw-orbits
 ifelse Orbits = "All" 
 [ask planets [pendown] ask turtle 7 [penup]]
 [ ifelse Orbits = "Mars" 
 [ask planets [penup] ask mars [pendown] ]
  [ifelse Orbits = "Venus" 
 [ask planets [penup] ask venus [pendown] ]
 [ifelse Orbits = "Jupiter"
 [ask planets [penup] ask Jupiter [pendown]] 
  [ifelse Orbits = "Mercury" 
 [ask planets [penup] ask mercury [pendown ]]
 [ifelse Orbits = "Saturn" 
 [ ask planets [penup] ask saturn [pendown] ] 
 [ ask planets [penup] ] ] ] ] ] ]
end
 
to orbit-sun 
 set time Speed / 100
 set scale Zoom * max-pxcor / 150
 ask planets [ ;;compute new positions for all orbiting bodies
 set newx newx + xvelocity * time 
 set newy newy + yvelocity * time
 ;;hide any turtle whose new position is off-screen
 ifelse abs (newx * scale) > max-pxcor or abs (newy * scale) > max-pxcor 
 [set hidden? true]
 [set hidden? false]
 set radius sqrt(newx ^ 2 + newy ^ 2)
 ;;compute new accelerations for all orbiting bodies based on sun's gravity. Equations are the same as in setup.
 set accx 0 - ((newx * [mass] of sun * G) / radius ^ 3) 
 set accy 0 - ((newy * [mass] of sun * G) / radius ^ 3)
 
 if Sun_only = false ;;unless only the sun's gravity is to be considered, each orbiting body asks all other orbiting bodies for their contribution to its acceleration
 [ 
 let my-accx accx
 let my-accy accy
  ask planets with [who != [who] of myself] [
 set my-accx my-accx - (((mass * G) * ([newx] of myself - newx))/(sqrt ( (newx - [newx] of myself ) ^ 2 + (newy - [newy] of myself) ^ 2)) ^ 3)
 set my-accy my-accy - (((mass * G) * ([newy] of myself - newy))/(sqrt ( (newx - [newx] of myself ) ^ 2 + (newy - [newy] of myself) ^ 2)) ^ 3)
 ]
 set accx my-accx
 set accy my-accy 
 ]
 ]
 ifelse Center ="Sun" [ ask planets [ setxy newx * scale newy * scale ]
 ask sun [setxy 0 0 ] ]
 [ let earth-x [newx] of earth
  let earth-y [newy] of earth
 ask planets [setxy (newx - earth-x) * scale (newy - earth-y) * scale]
 ask sun [setxy (- earth-x) * scale (- earth-y) * scale ] ]
 
 ask planets [ ;;compute new velocities
 set xvelocity xvelocity + accx * time 
 set yvelocity yvelocity + accy * time
 ]
 
 ifelse ViewLabels = false 
 [ask turtles [set label ""]]
 [ask planets [set label mylabel]]
 
 ;;If comet comes within .2 AU of Jupiter or Saturn and the display is off, set variable NearPlanet to true, which
 ;;causes the display to be turned on with the comet and planets in their updated positions and then held there for
 ;;two seconds so you can see what has happened and begin following the comet. 
 ;;Also, whether the display is on or not, while the comet is within the specified distance of either Jupiter or
 ;;Saturn, the distance between either or both and the comet is plotted.
 
 ask turtle 7 [
 ;;compute distance between comet and jupiter and saturn
 set dist2jupiter sqrt((newx - [newx] of jupiter) ^ 2 + (newy - [newy] of jupiter) ^ 2)
 set dist2saturn sqrt((newx - [newx] of saturn) ^ 2 + (newy - [newy] of saturn) ^ 2) 
 
 ifelse (dist2jupiter) < Comet-Jupiter ;if comet is within range of jupiter -
 [if display_on = false ;if display was off, set flag to turn it on below
 [set NearPlanet true]
 if jupiterplot = 0 and saturnplot = 0 ;if no plot for either Jupiter or Saturn is in progress
 [set-current-plot "Comet - Closest Planet Distance" ;clear the current plot
 clear-plot
 set jupiterplot 1 ;set flag indicating a plot is in progress
 ]
 set-current-plot "Comet - Closest Planet Distance"
 set-current-plot-pen "Jupiter" ;plot distance from comet to Jupiter
 plot dist2jupiter
 ]
 [set jupiterplot 0 ;if comet is out of range, set flag indicating no plot in progress
 ]
 ifelse (dist2saturn) < Comet-Saturn ;procedure below is same for Saturn as Jupiter
 [if display_on = false
 [set NearPlanet true]
 if saturnplot = 0 and jupiterplot = 0
 [set-current-plot "Comet - Closest Planet Distance"
 clear-plot
 set saturnplot 1
 ]
 set-current-plot "Comet - Closest Planet Distance"
 set-current-plot-pen "Saturn"
 plot dist2saturn
 ]
 [set saturnplot 0
 ]
 ]
 if NearPlanet = true ;check flag to see if display is to be turned on
 [set NearPlanet false ;after a two second delay
 ToggleDisplay
 wait 2]
 
 ask planets [ ;;keep track of number of orbits each planet makes
 if oldy < 0 and newy > 0 
 [set orbitcount orbitcount + 1]
 set oldy newy
 ]
 
 ;ask planets [set oldy newy]
 
 ;;keep track of ratio of current radius to initial average radius for both earth and comet so they can be plotted
 ask earth [set radius-ratio ((radius / (avg-radius)))]
 ask turtle 7 [set radius-ratio ((radius / (avg-radius)))]
 
 ; if Plotting = true [do-plot]
end 
 
to do-plot
 set-current-plot "Earth - Sun Distance"
 ask earth [set-current-plot-pen "radius-ratio earth"]
 ask earth [plot radius-ratio]
 set-current-plot "Comet - Sun Distance"
 ask turtle 7 [set-current-plot-pen "radius-ratio turtle 7"]
 ask turtle 7 [plot radius-ratio]
end
 
 
to setup-plots
 set-current-plot "Comet - Sun Distance"
 set-plot-y-range .8 1.2
 set-plot-x-range 0 1000
 set-current-plot "Earth - Sun Distance"
 ifelse Sun_Only = true ;;scale the earth orbit plot based on anticipated ranges. Determined by previous observations
 [set-plot-y-range precision (1 - (time)/ 10000) 4 precision (1 + (2 * time)/ 10000) 4]
 [set-plot-y-range precision (1 - (2 * sqrt(time))/ 1000) 4 precision (1 + (2 * sqrt(time))/ 1000) 4]
 set-plot-x-range 0 1000
 set-current-plot "Comet - Closest Planet Distance"
 set-plot-y-range 0 1
 set-plot-x-range 0 1000
end
 
to ToggleDisplay 
 ifelse display_on = true
 [no-display set display_on false]
 [display set display_on true]
end