This week, we create our own universe.  Specifically, we create a universe with two physical laws:

F = ma, or Newton’s second law of motion

F = G(m1*m2)/r^2, or Newton’s law of universal gravitation

Using those laws, we’re going to simulate a ‘solar system’ of different bodies and plot the paths they make moving under gravitation.  To do this we’re going to use an iterative algorithm; calculating the forces that the planets exert on each other, moving them a little bit, re-calculating the forces based on their new positions, moving them again and so on.

First up, we’re going to create a class to represent the ‘planets’ in our solar system.  We’re going to do this separately using Visual Studio, then compile it as a DLL, a library of code that does nothing on its own but that can be used freely by other bits of code.  There are a couple of advantages to doing it this way – we can re-use this bit of code easily if we need to and we can minimise our time spent using grasshopper’s horrible built-in script editor.

To create the dll, open up Visual Studio/Visual Basic Express and create a new ‘class library’ project.  DLLs are used a lot – RhinoCommon itself is stored in a DLL, and since we want to use some of the RhinoCommon types (points, vectors etc.) our first task is to reference that dll to let Visual Studio know where it is and that we want to use it.  To do this, go to Project->Properties->References and add one to RhinoCommon.dll (should be somewhere in your grasshopper install directory).  Then add this code:

'Before we can use anything from RhinoCommon, we need to reference the .dll
'You can do this in Project->Properties->References

'This imports the Rhino.Geometry namespace so we don't have to type it out every time we use something from there
Imports Rhino.Geometry

'This is the class we will use to represent each 'planet'
'This could be adapted to any kind of agent, however
Public Class Planet

'Each planet has three 'physical' properties
Public position As Point3d
Public velocity As Vector3d
Public mass As Double

'So we can plot the trail of the planet we store all of it's previous points
Public path As List(Of Point3d)

'Constructor:
Public Sub New(ByVal initialPosition As Point3d, ByVal initialVelocity As Vector3d, _
ByVal initialMass As Double)

position = initialPosition
velocity = initialVelocity
mass = initialMass

path = New List(Of Point3d)
path.Add(New Point3d(position.X, position.Y, position.Z))

End Sub

''' <summary>
''' This sub gets called iteratively to update the position of the planet
''' </summary>
''' <remarks></remarks>
Public Sub Move()
position += velocity
path.Add(New Point3d(position.X, position.Y, position.Z))
End Sub

''' <summary>
''' This converts a force into a change in velocity
''' </summary>
''' <param name="force"></param>
''' <remarks></remarks>
Public Sub ApplyForce(ByVal force As Vector3d)
velocity += force / mass
End Sub

End Class

Then hit ‘Build’ to create our actual .dll. VS will put it in the folder for your project under Bin/Release.

Now for our grasshopper part.  Create a script component, right click on it, go to ‘Referenced assemblies’ and use the pop-up form to add a reference to our newly-created .dll.  This will allow us to use our ‘Planet’ class within the scripting component.

Set up the component’s inputs to take a list of points, P, a list of vectors, V, a list of doubles, M, and an integer, T.  These represent, respectively, our planet’s starting positions, initial velocities, masses and the length of time we want to simulate.  Then add this code:

Private Sub RunScript(ByVal P As List(Of Point3d), ByVal V As List(Of Vector3d), ByVal M As List(Of Double), ByVal T As Integer, ByRef Path As Object)

Dim planets As New List(Of Planet)

'First of all, we create our 'planets'.
'Here we recreate the behaviour Of the 'longest list' data matching option
For i As Integer = 0 To P.Count - 1
Dim position As Point3d = P(i)
'If values exist for the other parameters here then we use them
'otherwise, we use the last value in the list
Dim velocity As vector3d = V(math.Min(i, V.Count - 1))
Dim mass As Double = M(math.Min(i, M.Count - 1))
Dim planet As New Planet(position, velocity, mass)
planets.Add(planet)
Next

'All our simulation is going to take place in this loop
'Each step through this loop represents a discrete amount of time
'To simplify things, we've assumed this is equal to 1 time unit
For time As Integer = 0 To T

'First, we check every planet against every other planet and work out the forces each
'pair is exerting on each other
For j As Integer = 0 To planets.Count - 2

Dim planetA As Planet = planets(j)

'By checking each planet only against the planets later in the list, we ensure that each
'pair is only going to be checked once
For k As Integer = j + 1 To planets.Count - 1

Dim planetB As Planet = planets(k)
Dim F As Double
'Here we calculate the force using the formula:
'Force = (mass1 * mass2)/distance^2
'i.e. Newton's law of universal gravitation
F = (planetA.mass * planetB.mass) / planetA.position.DistanceTo(planetB.position) ^ 2
Dim direction As Vector3d = planetB.position - planetA.position
direction.Unitize
direction *= F

planetA.ApplyForce(direction)
planetB.ApplyForce(-direction)
Next

Next

'Now that we have worked out how each planet will be moving, we update their positions
For Each planet As Planet In Planets
planet.Move()
Next

Next

'NOTE: This part *should* work, but doesn't always due to a bug in RhinoCommon:
'=========================================================================================
'Dim outList As New List(Of Curve)
'For Each planet As Planet In Planets
'  Dim crv As Curve = Curve.CreateInterpolatedCurve(planet.path, 3)
'  outList.Add(crv)
'Next

'Path = outList
'=========================================================================================

'Since the above won't work reliably, we instead create a data tree, with the points that
'describe each planet's path stored in a separate branch.  This means we can feed them into
'an 'interpolate curve' component and get a separate curve for each planet.
Dim Tree As New DataTree(Of Point3d)

For z As Integer = 0 To planets.Count - 1
'For each planet we create a new path:
Dim TreePath As New GH_Path(z)
For Each pt As point3d In planets(z).path
'We then add each point to the tree, stipulating the path to use to store them
Tree.Add(pt, TreePath)
Next

Next

'And finally, we output our tree:
Path = Tree

End Sub 

Now plug the output path points into an ‘interpolate curve’ component, stick in some initial variables and you should have some interesting results, or at least interesting messes.  Note that the algorithm is quite sensitive to scale – if you don’t see anything it may just be that your planets are too far apart to have any effect on each other.

You can download the complete example here: http://www.vitruality.com/Teaching/PhysicsExample.zip

I used a script like this when we were in the initial stages of designing the Orbit as a way of defining the centreline of the lattice.  In the end however, this approach proved a bit to difficult to control.  If you play around with it a bit you will notice that it is extremely sensitive to initial conditions.  Move one of the starting positions even slightly and it will totally change the result you get.  Since for the Orbit we needed to have finer control over the form in order to get something that was structurally feasible (we were going for visually unstable, but not actually unstable), in the end we abandoned that idea and went for something a bit more manual.

But, while this script might not be particularly useful for this reason, this basic structure could be modified for any sort of physical simulation or agent-based generation (dynamic relaxation, boids, pheremone-based stuff, etc.).  The only things that you’d need to change would be the controlling equations (in this case, the gravitation formula)…