This page was automatically generated by NetLogo 5.0.
The applet requires Java 5 or higher. Java must be enabled in your browser settings. Mac users must have Mac OS X 10.4 or higher. Windows and Linux users may obtain the latest Java from Sun's Java site.
powered by NetLogo
view/download model file: solar_system.nlogo
This is a model of a simple solar system with one sun and one planet, for the book entitled “Physicomimetics: Physics-Based Swarm Intelligence.”
Two particles use F = ma and Newton’s gravitational force law to simulate a planet orbiting a sun. The sun has a much higher mass than the planet.
Click SETUP AGENTS to initialize the sun and planet, and click MOVE AGENTS to have them move. The sun is the yellow dot and the planet is the white dot.
The ANGULAR_MOTION slider allows you to impart a spin to the system at initialization, causing the planet to orbit the sun. The amount of ANGULAR_MOTION has a very large impact on the shape of the orbit. The spin is established very carefully to make sure that there is no linear momentum at initialization. Changing this slider while the simulation is running will have no effect.
The MASS_OF_SUN slider allows you to control the mass of the sun (the planet has a mass of one). This slider is continually monitored, so moving the slider while the simulation is running will affect the simulation.
All other sliders will affect the simulation when it is running.
This simulation serves to teach you a simple model of a solar system, as well as to continue with more advanced topics, such the Conservation of Linear and Angular Momenta and the Conservation of Energy. This is covered in detail in Chapter 2 of the book.
There is no random component to this code, so it will run the same way each time, if you don’t change the settings.
There are additional monitors in this simulation. One computes the current total energy of the system. Another displays the velocity of the planet. The third shows the magnitude of the force on the planet.
This simulation introduces two new sliders, FORCE_MAXIMUM and VELOCITY_MAXIMUM, to reflect that when working with robots, they can only achieve some maximum velocity. When using a 1/r^2 force law, the force can potentially go to infinity, creating velocities that no real robot can obtain. But what if the robot can’t go that fast? How should we deal with that? We investigate two ways: (1) limiting the velocity, and (2) limiting the force.
Note how lowering the TIME_STEP causes the total energy monitor to fluctuate less, indicating that the Conservation of Energy is holding better. In the graph, the total energy is shown in brown, the kinetic energy is in green, and the potential energy is in blue. The total energy stays relatively constant while there is a constant tradeoff between potential and kinetic energy. This is best seen when FORCE_MAXIMUM is set to one and ANGULAR_MOTION is set to 10. Watch the simulation - when is kinetic energy high? When is potential energy high?
The red dot in the simulation shows the center of mass of the system. If the Conservation of Linear Momentum holds in both the x- and y-dimensions the red dot will not move. The graph shows that the linear momenta stay extremely close to zero, which is what we should expect if the system is programmed properly. This simulation also includes a monitor for the Angular Momentum and you will see that it does not change over time. If, however, the planet crosses the boundary of the world (re-entering from the other side) the conservation laws can be broken, because the standard physics assumption of an Euclidean geometry no longer holds (as explained in Chapter 2).
See how the GRAVITATIONAL_CONSTANT changes behavior.
Try different values of ANGULAR_MOMENTUM (even zero).
Change the FORCE_MAXIMUM and VELOCITY_MAXIMUM. What happens?
What if you change the MASS_OF_SUN?
It is somewhat hard to set the sliders so that you get orbits that stay within the graphics pane. If you are running the simulation on your computer directly (by using “solar_system.nlogo”) you can increase the pane size by clicking on “settings” and then increasing max-pxcor and max-pycor.
Introduce a second planet. How difficult is it to create a stable solar system?
Note, in order to change any NetLogo simulation, you must have the source code (i.e., “solar_system.nlogo”) downloaded to your computer, as well as NetLogo itself. You can not change the code when you are running the simulation with your browser.
Since we are using a patch size of one, we wanted the particles to be more visible. This is done with “set size 5” in the code. However, they are still considered to be point particles (with no size) in the simulation.
This simulation introduces the “pen down” (pd) command to draw the orbit of the planet. You could make a similar change to the “spring2D” simulation also, if you want. The new CLEAR button allows you to erase this orbit.
Note how the “do-plots” procedure draws the Energy graph and the Momenta graph.
Now we can see why we didn’t use the built-in NetLogo commands to model springs. By modeling springs from first principles, it was trivial to change the previous two-dimensional spring model to use Newton’s gravitational force law instead of Hooke’s spring law. Working from first principles allows us to have more flexibility in what we are modeling.
This is our third simulation, which builds on the two-dimensional spring model. It will be generalized more and more throughout the book.
For a similar analysis of potential energy, see:
Spears, W. M., Spears, D. F., Hamann, J., and Heil, R. (2004) Distributed, physics-based control of swarms of vehicles. Autonomous Robots, 17 (2-3).
If you mention this model in an academic publication, we ask that you include these citations for the model itself and for the NetLogo software:
- Spears, William M. and Spears, Diana F. (eds.) Physicomimetics: Physics-Based Swarm Intelligence, Springer-Verlag, (2011).
- Wilensky, U. (1999). NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
Copyright 2011 William M. Spears. All rights reserved.
Permission to use, modify or redistribute this model is hereby granted, provided that both of the following requirements are followed:
a) this copyright notice is included, and
b) this model will not be redistributed for profit without permission from William M. Spears. Contact William M. Spears for appropriate licenses for redistribution for profit.
; William M. Spears September 2011
; Solar System Tutorial Code
; For research and educational use only
globals [total_lmx total_lmy total_angular_mom total_ke total_pe total_energy
FMAX VMAX center_of_mass_x center_of_mass_y rprime G DeltaT S]
turtles-own [hood deltax deltay r F Fx Fy v vx vy dvx dvy mass ke
lmx lmy theta lever_arm_x lever_arm_y lever_arm_r angular_mom]
to setup
clear-all ; Clear everything
crt 2 ; Create two turtles (particles)
update-info
ask turtles [setup-turtles] ; Set up the two particles
; Compute initial separation
set S sqrt (((([xcor] of turtle 1) - ([xcor] of turtle 0)) *
(([xcor] of turtle 1) - ([xcor] of turtle 0))) +
((([ycor] of turtle 1) - ([ycor] of turtle 0)) *
(([ycor] of turtle 1) - ([ycor] of turtle 0))))
; Computes center of mass and displays location
set center_of_mass_x (([mass] of turtle 0) * ([xcor] of turtle 0) +
([mass] of turtle 1) * ([xcor] of turtle 1)) /
(([mass] of turtle 0) + ([mass] of turtle 1))
set center_of_mass_y (([mass] of turtle 0) * ([ycor] of turtle 0) +
([mass] of turtle 1) * ([ycor] of turtle 1)) /
(([mass] of turtle 0) + ([mass] of turtle 1))
ask patch (round center_of_mass_x) (round center_of_mass_y)
[ask patches in-radius 4 [set pcolor red]]
reset-ticks
end
to run-and-monitor
if (count turtles < 1) [user-message "Please click HALT and then SETUP AGENTS first" stop]
update-info
ask turtles [ap]
ask turtles [move]
; Computes center of mass and displays location
set center_of_mass_x (([mass] of turtle 0) * ([xcor] of turtle 0) +
([mass] of turtle 1) * ([xcor] of turtle 1)) /
(([mass] of turtle 0) + ([mass] of turtle 1))
set center_of_mass_y (([mass] of turtle 0) * ([ycor] of turtle 0) +
([mass] of turtle 1) * ([ycor] of turtle 1)) /
(([mass] of turtle 0) + ([mass] of turtle 1))
ask patch (round center_of_mass_x) (round center_of_mass_y)
[ask patches in-radius 4 [set pcolor red]]
set total_lmx sum [lmx] of turtles ; Total linear momentum, x-component
set total_lmy sum [lmy] of turtles ; Total linear momentum, y-component
set total_angular_mom sum [angular_mom] of turtles
set total_ke sum [ke] of turtles ; Total kinetic energy of both particles
; The following code is used to calculate the potential energy of the system
; The paper mentioned in the References of the Information tab provides a similar
; derivation of potential energy
set rprime sqrt (G * ([mass] of turtle 0) * ([mass] of turtle 1) / FMAX)
ifelse (S >= rprime)
[set total_pe ((G * ([mass] of turtle 0) * ([mass] of turtle 1) * ((1 / rprime) - (1 / S))) +
(FMAX * rprime))]
[set total_pe (FMAX * S)]
set total_energy (total_ke + total_pe) ; Total energy of two particle system
tick
do-plots
end
to setup-turtles
home set shape "circle" set size 5 set vx 0 ; Start with no x motion on either particle
ifelse (who = 0) ; To allow for unequal masses of the two particles
[set color yellow set mass Mass_of_Sun set vy ((- Angular_Motion) / Mass_of_Sun)]
[set color white set xcor 15 set mass 1 set vy Angular_Motion pd]
set theta 0
end
to ap ; Run artificial physics
if (who = 0) [set mass Mass_of_Sun] ; Monitor this slider always
set hood [who] of other turtles ; Get the IDs of your neighbors
foreach hood [
set deltax (([xcor] of turtle ?) - xcor)
set deltay (([ycor] of turtle ?) - ycor)
set r sqrt (deltax * deltax + deltay * deltay)
set S r
; Newton's gravitational force law
set F (G * mass * ([mass] of turtle ?) / (r * r))
if (F > FMAX) [set F FMAX]
set Fx (F * (deltax / r)) ; The x-component of force
set Fy (F * (deltay / r)) ; The y-component of force
]
set dvx DeltaT * (Fx / mass)
set dvy DeltaT * (Fy / mass)
set vx (vx + dvx) ; The x-component of velocity
set vy (vy + dvy) ; The y-component of velocity
set v sqrt (vx * vx + vy * vy)
if (v > VMAX) [set vx (VMAX * vx / v) set vy (VMAX * vy / v) set v VMAX]
set deltax DeltaT * vx
set deltay DeltaT * vy
if ((deltax != 0) or (deltay != 0))
[set heading (atan deltax deltay)] ; Because heading = 0 means turtle faces straight up!
end
to move ; Move the turtle
fd (sqrt (deltax * deltax + deltay * deltay))
set lmx (mass * vx) ; Linear momentum of the turtle
set lmy (mass * vy)
set ke (v * v * mass / 2) ; Kinetic energy of the turtle
set lever_arm_x (xcor - center_of_mass_x)
set lever_arm_y (ycor - center_of_mass_y)
set lever_arm_r sqrt (lever_arm_x * lever_arm_x + lever_arm_y * lever_arm_y)
if (((vx != 0) or (vy != 0)) and ((lever_arm_x != 0) or (lever_arm_y != 0)))
[set theta (atan (mass * vy) (mass * vx)) - (atan lever_arm_y lever_arm_x)]
set angular_mom (lever_arm_r * mass * v * (sin theta)) ; Angular momentum of the turtle
end
to update-info ; Update information from the sliders
set G Gravitational_Constant
set FMAX Force_Maximum
set VMAX Velocity_Maximum
set DeltaT Time_Step
end
to do-plots
set-current-plot "Energy" ; Select the Energy plot
set-current-plot-pen "Total" ; Select the Total Energy Pen
plot total_energy ; Plot the total_energy
set-current-plot-pen "Potential"
plot total_pe ; Plot the potential energy
set-current-plot-pen "Kinetic"
plot total_ke ; Plot the kinetic energy
set-current-plot "Linear and Angular Momenta"
set-current-plot-pen "Lmx"
plot total_lmx ; Plot the linear momentum, x-component
set-current-plot-pen "Lmy"
plot total_lmy ; Plot the linear momentum, y-component
set-current-plot-pen "Angular"
plot total_angular_mom ; Plot the angular momentum
end