# Particle Spring Simulation Tutorial

A tutorial for the course Digital World 2015 demonstrating a bare simple particle spring simulation for Processing/Python. To run the code below, download Processing and install the Python editing mode. The tutorial contains a series of suggestions/challenges for improving this basic template.

``````#--
#-- Particle Spring Dynamics Simulation
#--

#-- Window Dimensions
#--
dx = 800
dy = 600

#-- Definition of Complex Types ---------------------------------------------
#--
#-- As not sure if you are familiar with such concepts as a class or a tuple,
#-- I use below vanilla lists to express semantically complex data forms such
#-- as a vector, a particle and a spring.
#--
#-- The convention used is: Any complex type of information will be stored in
#-- a list each element thereof representing a property in orderly manner.
#--
#-- The name of the property will be globally defined as a constant containing
#-- the index of the property value in the list.
#--
#-- Example: A point/vector is a complex entity which contains a number of
#-- coordinates. For 2D vectors we need two entries the first one being X
#-- and second Y coordinate. Below is the definition of X = 0 and Y = 1
#-- Therefore I can get/set the X coordinate of point P using P[X]
#--

#-- Vector Definition -------------------------------------------------------
#--
#-- A vector represents a position, displacement or force in two dimensions.
#-- Property definitions (X,Y coordinates) and basic algebra for vectors are
#-- expressed below as global functions.
#--
#-- CHALLENGE: Modify the definion to support 3D points/vectors (easy)
#--
X = 0
Y = 1

#-- Create a New Vector
#--
def vnew( x, y ):
return [x, y]

#--
return vnew( u[X] + v[X], u[Y] + v[Y] )

#-- Subtract Two Vectors
#--
def vsub( u, v ):
return vnew( u[X] - v[X], u[Y] - v[Y] )

#-- Vector Scale
#--
def vmul( u, s ):
return vnew( u[X] * s, u[Y] * s )

#-- Vector Dot Product
#--
def vdot( u, v ):
return u[X] * v[X] + u[Y] * v[Y]

#-- Vector Length
#--
def vlen( u ):
return sqrt( vdot( u, u ) )

#-- Particle Definition -----------------------------------------------------
#--
#-- A particle represent a dimensionless mass with internal dynamics state of
#-- position and velocity. The force acting upon the particle can be also
#-- though of as its acceleration. In other words the particle captures a set
#-- of information about a vector valued function and its first and second
#-- derivatives.
#--
MASS     = 0
POSITION = 1
VELOCITY = 2
FORCE    = 3

#-- Create Particle
#--
def pnew( m, p, v, f ):
return [m, p, v, f]

#-- Spring Definition -------------------------------------------------------
#-- A spring represents an idealized elastic link between two particles that
#-- obays Hooke's Law of springyness: F = k * d, where F is the force magnitude
#-- k is the spring constant, and d is the displacement.
#--
#-- Alternatively you can think of this as: F = k * ( RestLength - CurrentLength )
#-- where RestLength is the length at which the spring excerts zero force to both
#-- particles (aka ends of the spring locations) and CurrentLength is the length
#-- which is either bigger (spring is stretched) or smaller (spring is compressed)
#-- than the rest length (if they are equal then F = 0).
#--
#-- In terms of properties we store the source and target locations (the end-points)
#-- of the spring such that we can connect them with other springs or fix them etc.
#-- We obviously need the spring constant and rest length.
#--
#-- NOTE: The source and target properties are expressed as indices (integers)
#--       instead of particle objects (list of properties). The conventions is
#--       made such that we can store multiple particles in one list (each having
#--       a known index) and define springs that associate/point-into the particle
#--       collection by index (instead of duplicating particles per spring).
#--
#-- CHALLENGE #1: Can we change the representation from the high-school notion
#--               of the spring, to the proper engineering idea of elasticity? (easy)
#--               HINT: google stress/strain and replace k with the elastic modulus e
#--
#-- CHALLENGE #2: The spring, apart from its direct reaction to loading expressed
#--               by Hooke's Law, it also has a force component, not used above, which
#--               captures its internal damping. Otherwise a spring would oscillate
#--               indefinately if think about it. Can express this also (medium)
#--               HINT: google linear spring damping
#--
SOURCE    = 0  #-- Source Particle INDEX
TARGET    = 1  #-- Target Particle INDEX
STIFFNESS = 2  #-- Spring Stiffness Constant
LENGTH    = 3  #-- Rest Length

def snew( s, t, k, l ):
return [s, t, k, l]

#-- Particle and Spring Lists -----------------------------------------------
#--
particles = []
springs   = []

#-- Create Initial Model ----------------------------------------------------
#--
#-- This function is required by processing in order to prepare the program.
#-- It will be the first thing to be executed before any drawing takes place.
#-- In here we will construct a chain of particles and connect them with
#-- springs. They will be stored in the global variables "particles" and
#-- "springs" for later use in the simulation/animation cycle aka draw( )
#--
#-- NOTE: There is a conceptual leap here... it is only a convention or some
#--       sort of approximation involved in the assumption that a sequentially
#--       connected set of particles represents a "physical" chain. In the
#--       fashion we may "define" a rectangular grid of connected particles
#--       as a "fabric" sheet. More on that later, meanwhile think of other
#--       "things" that can be represented by such configurations of particles
#--       as springs.
#--
def setup( ):

#-- Window Size
#--
size( dx, dy )
frameRate( 24 )

#-- Import Globals
#--
global particles
global springs

#-- Number of Particles
#--
count = 12

#-- Construct Particles ---------------------------------
#-- Create a sequence of particles along the x-direction
#-- NOTE: We need to provide mass, possition and velocity
#--       The last two items are important as they define
#--       the dynamics state of the particle. You can set
#--       velocity to non-zero if you like.
#--
for index in range( 0, count ):
mass = 20.0

x = index * 50 + 100
y = 100

position = vnew( x, y )
velocity = vnew( 0, 0 )
force    = vnew( 0, 0 )

particles.append( pnew( mass, position, velocity, force ) )

#-- Construct Springs -----------------------------------
#-- Connect previously constructed particles in order
#--
for index in range( 1, count ):
source_index = index - 1
target_index = index - 0

stiffness   = 15.0
rest_length = 50.0

springs.append( snew( source_index, target_index, stiffness, rest_length ) )

#-- By own convention, particles with zero mass are assumed fixed in place
#-- There are other ways to define such a constraints but for simplicity
#-- we define here as such. The first and last particles of the chain are fixed.
#--
particles[        0][MASS] = 0.0
particles[count - 1][MASS] = 0.0

#-- Draw Window Contents -----------------------------------------------------
#-- This is also a standard function required by processing. It is executed
#-- every some odd milliseconds and it is supposed to draw the screen (aka FPS)
#--
#-- In here we perform two tasks: 1. We run the simulation of particles and
#-- springs. 2. We draw them in the window as we are supposed.
#--
#-- Logic of Simulation: the steps taken for simulation are...
#-- 1. Force Calculation/Accumulation     F = m * a  ->
#--                                       a = F / m
#--
#-- 2. Velocity update from acceleration  a = dv / dt ->
#--                                       dv = a * dt ->
#--                                       v_new - v_old = a * 1.0 ->  (!!! dv = v_new - v_old )
#--                                       v_new = v_old + a
#--
#-- 3. Position update from velocity      v  = dp / dt ->
#--                                       dp =  v * dt ->
#--                                       p_new - p_old = v * 1.0 ->  (!!! dp = p_new - p_old )
#--                                       p_new = p_old + v
#--
#-- 4. Enforce Constraints (such that particle dont fall off the screen)
#--
#-- NOTE: The above logic is called forward explicit Newton integration.
#--       The places with (!!!) are important because there are assumptions
#--       of approximating differential equations with their first linear
#--       term only (which will cause noticable inaccuracy, unless dt is small)
#--
#--       Newton integration is simple and fast but tends to spiral out of
#--       control due to numerical errors (we are just following the tangent
#--       line forward instead of the curved true trajectory, picture it).
#--
#--       The second assumption is that the time step dt = 1.0 which is a
#--       simplification that couples the animation FPS with simulation time.
#--       You may slow down the simulation by changing the FPS or actually
#--       multiplying velocity and acceleration by a small scalar eg. 0.0001
#--       Scaling the time step (dt) is advisable as it will improve accuracy.
#--
#-- CHALLENGE: Add a proper time step component. HINT: basic vector math (easy)
#--
#-- CHALLENGE: There other more sophisticated integration methods which are
#--            also explicit but more accurate eg. Mid-Point (easy), Verlet (easy)
#--            RK4 (tricky), as well as implicit methods using linear algebra (hard)
#--
#-- CHALLENGE: So we know from theory that the shape of handing chain is
#--            expressed-by / has-a closed form equation known as the catenoid.
#--            Draw it over the simulated shape (same parameters) and compare the
#--            results. What is the significance of the error?
#--
#-- CHALLENGE: So we have been using abstract units for space (pixels), mass(units)
#--            time (fps). What if you apply real units eg. meters, kilograms, seconds
#--            such that you can compare the real shape of a physical chain with the
#--            simulation. How does the real shape (take photo and overlay), theory
#--            and simulated form compare?
#--
def draw( ):

#-- Import Globals
#--
global particles
global springs

#-- Reset Force Accumulator ---------------------------------
#--
for particle in particles:
particle[FORCE] = vnew( 0.0, 0.0 )

#-- Accumulate Gravity Forces -------------------------------
#-- Assuming earth gravity constant but note y is positive
#-- because y-screen direction is pointing downwards already.
#--
GRAVITY = 9.81

for particle in particles:
gravity = vnew( 0, GRAVITY * particle[MASS] )
particle[FORCE] = vadd( particle[FORCE], gravity )

#-- Accumulate Damping Forces -------------------------------
#-- Fictional linear drag opposite to particle's velocity
#-- for stability. Air drag is quadratic to velocity and
#-- proportional to fluid density as well as surface area.
#--
#-- CHALLENGE: Change this unrealistic drag force to a proper
#--            air drag force (easy). HINT: google drag force
#--            when is doubt about areas & constants use 1.0
#--
DAMPING = 2.00
for particle in particles:
damping = vmul( particle[VELOCITY], -DAMPING )
particle[FORCE] = vadd( particle[FORCE], damping )

#-- Accumulate Spring Forces -------------------------------------------
#--
for spring in springs:

#-- Extract Particle Indices
#--
source_index = spring[SOURCE]
target_index = spring[TARGET]

#-- Extract Particles (themselves)
#--
source_particle = particles[source_index]
target_particle = particles[target_index]

#-- Extract Positions
#--
source_position = source_particle[POSITION]
target_position = target_particle[POSITION]

#-- Hook's Law (over-simplified)
#-- The proper formulation includes a damping
#-- component proportional to velocity along the
#-- spring direction. NOTE: I am dividing by the
#-- distance between the spring's end-points aka
#-- the current length, because the force magnitude
#-- is assuming a normalized directional vector
#--
stiffness   = spring[STIFFNESS]
rest_length = spring[LENGTH]
direction   = vsub( target_position, source_position )
distance    = vlen( direction )
magnitude   = stiffness * ( distance - rest_length ) / distance

source_particle[FORCE] = vadd( source_particle[FORCE], vmul( direction,  magnitude ) )
target_particle[FORCE] = vadd( target_particle[FORCE], vmul( direction, -magnitude ) )

#-- Newton Integration -----------------------------------------
#-- Discrete integration by assuming time step of one unit.
#-- For better accuracy and stability either scale acceleration
#-- and velocity vectors by a time step eg. dt = 0.001 or use
#-- better integration approximation.
#--
for particle in particles:

#-- Zero mass particles are ignored
#--
mass = particle[MASS]

if( mass > 0.0 ):
#-- F = ma <=> a = F / m
#--
acceleration = vmul( particle[FORCE], 1.0 / mass )

#-- acceleration = dv/dt ~> v_new = v_old + acceleration * dt
#--
particle[VELOCITY] = vadd( particle[VELOCITY], acceleration )

#-- velocity = dp/dt ~> p_new = p_old + velocity * dt
#--
particle[POSITION] = vadd( particle[POSITION], particle[VELOCITY] );

#-- Constraint Enforcement -----------------------------------
#--
#-- Force particles to remain within window boundaries by
#-- reversing their position to previous and flipping velocity
#--
#-- This is kind of a perfect elestic collision responce since
#-- velocity is preserved 100%
#--
#-- The reason for reversing position to previous is because
#-- the particle may have escaped significantly outside. So
#-- better safe than sorry.
#--
#-- Constraints are better expressed as "force constraints"
#-- rather than "penalties". If you think about it, the code
#-- below is not 100% realistic because every time we move
#-- particles back and flip their velocity direction we are
#-- fiddling with the energy/momentum of the system in ways
#-- that may violate constitutional laws of coservation.
#--
#-- implement a constraint such as particle-on-circle (hard)
#--
for particle in particles:

#-- Particle Properties
#--
position = particle[POSITION]
velocity = particle[VELOCITY]
radius   = particle[MASS] * 0.5  #-- Drawing artifact, see below

#-- Bound Checking
#--
if( ( position[X] < radius ) or ( position[X] > ( dx - radius ) ) ):
position[X] = position[X] - velocity[X]
velocity[X] =-velocity[X]

if( ( position[Y] < radius ) or ( position[Y] > ( dy - radius ) ) ):
position[Y] = position[Y] - velocity[Y]
velocity[Y] =-velocity[Y]

#-- Prepare Frame -------------------------------------------
#-- Clears the background, otherwise it will look psychedelic
#-- try commenting out the background line below and enjoy
#--
fill( 0 )
background( 255 )

#-- Draw Springs --------------------------------------------
#--
for spring in springs:

source_index = spring[SOURCE]
target_index = spring[TARGET]

source_particle = particles[source_index]
target_particle = particles[target_index]

source_position = source_particle[POSITION]
target_position = target_particle[POSITION]

stroke( 0, 0, 255 )
line( source_position[X], source_position[Y],
target_position[X], target_position[Y] );

#-- Draw Particles ---------------------------------------------
#--
noStroke( )
index = 0
for particle in particles:

position = particle[POSITION]
mass     = particle[MASS]

#-- Fixed particles (with zero mass) have green color
#-- RGB=(0,255,0), unconstrained particles have red
#-- color (0,255,0) and radius = half their mass
#--
if( mass > 0.0 ):
fill( 255, 0, 0 )
else:
fill( 0, 255, 0 )