h4tt3n's spring physics simulation tutorial

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

h4tt3n's spring physics simulation tutorial

Post by h4tt3n »

1. INTRODUCTION

Hello everyone and welcome to this FreeBasic tutorial on spring simulation, or more exactly mass-spring-damper simulation. This is still somewhat of a rough draft covering the basics, and I intend to refine it and add some more stuff. Please do post a comment if you feel like it!

Cheers, Michael


2. THE BASICS OF SPRING PHYSICS

This is definitely not intended to be a physics lesson. In other words I'll keep the strictly scientific part at a minimum and focus on how to implement pysics in code rather than do a lengthy and boring lecture on mechanics.

Lets start by looking at how to represent a spring within a computer program. The simplest way to do this is by setting up two endpoint masses in a coordinate system and connect them to each other with a massless spring, like this:

Image



This is a very common way of doing it, and even in very advanced simulations it is most likely done in the same way. You can do this in one, two, or three dimensions depending on your needs, but the physics behind it is still the same.

The spring's job is to keep the two masses at a certain distance from each other at all times by pulling or pushing them. In the following we'll se how this is done from a physical viewpoint. In this tutorial we'll use Newton's three laws of motion when moving objects around:


1. Every object in a state of uniform motion (moving or at rest) remains in that state unless an external force is applied to it.

2. The relationship between an object's mass m, its acceleration a, and the applied force F is F = m * a.

3. For every force (action) there is an equal and opposite force (reaction).



Besides these three laws all you need to know is that


- Acceleration equals change of velocity per time unit: a = dv / dt

and

- Velocity equals change of position per time unit: v = dx / dt


When moving the endpoint masses around in virtual space our "path of calculation" is:

Force --> Acceleration --> Velocity --> Position

In relation to spring physics, this is how it works:

First we calculate the spring force, using the so called Hooke's law of elasticity. This equation approximates the force applied by a spring as a function of its deformation

F = -K * X

where F = Force, K = spring stiffnes constant, and X = the difference between the springs actual length and its neutral length or "rest length". A more reader friendly way of writing it might be:

Force = -stiffnes * (length - rest length)

The spring stiffnes constant is simply an arbitrary positive number which defines the springs hardness. A low value gives a soft, rubbery spring, and a high value results in rigid, steel-like behaviour.

The spring's rest length is usually a predefined constant. The spring force equation (Hooke's law) will try to keep the spring at that length by applying a force on the two endpoint masses, thus accelerating them towards or away from each other. This is how it works:

If the spring's length is identical to its rest length, then the expression

length - rest length = 0

which means no force is applied on the end masses. If the spring is squeezed shorter than its rest length then

length - rest length < 0

which means that F > 0 ie. an "outwards" force is applied to the two endpoint masses, pushing them away from each other. If on the other hand the spring is stretched longer than its rest length, then

length - rest length > 0

which results in a negative force value ie. an "inwards" force, which will try to pull the endpoint masses back towards each other. Here is an illustration:

Image

Once we've calculated the right force, we need to convert it into acceleration. From Newton's 2nd law we know that

Force = mass * acceleration

which means that

acceleration = Force / mass

Please note that if the two endpoints have different masses they will experience different acceleration, even though they are influenced by the same force. So, although we only need to calculate one force, we still need to calculate the acceleration of each mass separately:

acceleration(endpoint 1) = -Force / mass(endpoint 1)

acceleration(endpoint 2) = Force / mass(endpoint 2)

According to Newton's 3rd law each force is accompanied by an equal opposite force, hence the terms "Force" and "-Force" in the equations.

From the acceleration we need to find the change of velocity of the spring endpoint masses. Here we use the fact that acceleration equals change of velocity per time step

a = dv / dt

Which means that

dv = a * dt

Where a = acceleration, dv = delta velocity (change in velocity), and dt = delta time (time step).

Delta time or time step is simply a positive numerical value which decides how fast the simulation runs. Smaller values result in a slower but more accurate simulation, while bigger values means faster simulation speed at the cost of acccuracy.

Please note that we haven't actually calculated velocity, but change in velocity. The new velocity value is simply found by adding the delta velocity value to the current velocity

v(new) = v(old) + dv

Finally we need to find the new position of the masses based on velocity. Here we need to know that velocity equals change in position per time step

v = dx / dt

which means that

dx = v * dt

Where dx = delta position (or change in position), v = velocity, and dt = delta time.

Again, since we've only found the change in position, the actual new position is found by adding delta position to the current position value

x(new) = x(old) + dx

So, finally we've come all the way from spring deformation through force, acceleration, velocity, and to the resulting new positions of the spring endpoint masses!

And now I think it's time to look at an example program. Here is a simple one dimensional spring simulation which implements what we've gone through in the above.
Play around with it until you feel comfortable with the principles that makes it work. Also, try to find the limits of the simulation. For instance, try to find the highest acceptable time step and spring stiffnes values.

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's mass-spring-damper tutorial, October 2007
'
'   Example # 1
'
'   Simple one dimensional spring simulation. The spring is represented both 
'   graphically and numerically. Change the constants defining the spring and 
'   see what happens! (Press esc to quit)
'
'******************************************************************************'

''  set constants. experiment with these and see how the simulation reacts
const Pi = 4*Atn(1)             ''  Pi (better not change this ^^)
const Mass_1_mass = 100         ''  set mass of the spring's endpoints
const mass_2_mass = 20          ''  set mass of the spring's endpoints
const spring_distortion = 50    ''  how stretched the spring is at startup
const spring_rest_length = 200  ''  spring neutral length or rest length
const spring_stiffnes = 50      ''  spring stiffnes. Larger equals stiffer
const timestep = 0.01           ''  delta time. Smaller is slower but more accurate

''  dimension variables
Dim as Double time_passed, Spring_length, Force, _
  Mass_1_Acceleration, Mass_2_Acceleration, _
  Mass_1_Velocity, Mass_2_Velocity, _
  Mass_1_Position, Mass_2_Position, _
  Mass_1_Radius, Mass_2_Radius

''  set startup position of masses
Mass_1_Position = 300-(spring_rest_length+spring_distortion)/2
Mass_2_Position = Mass_1_Position + spring_rest_length + spring_distortion

''  set mass radius using the equation of a sphere's volume
Mass_1_Radius = ((Mass_1_mass/0.01)/((4/3)*pi))^(1/3)
Mass_2_Radius = ((Mass_2_mass/0.01)/((4/3)*pi))^(1/3)

''  set screen to 600 * 400 pixels * 16 bits per pixel
Screenres 600, 400, 16

''  main program loop
Do
  
  ''  first we find the length of the spring
  spring_length = Mass_2_Position - Mass_1_Position
  
  ''  then we calculate the force with Hooke's law of elasticity
  Force = -Spring_stiffnes * (spring_length - spring_rest_length)
  
  ''  now we can find the acceleration based on Newton's 2nd law, a = F/m
  ''  Since the two masses may be different we have to calculate their 
  ''  acceleration separately.
  ''  (According to Newton's 3rd law each force is accompanied by an equal 
  ''  opposite force, hence the terms "Force" and "-Force" in the equations.)
  Mass_1_Acceleration = -Force / mass_1_mass
  Mass_2_Acceleration = Force / mass_2_mass
  
  ''  then we find change in velocity and add it to the existing velocity
  mass_1_Velocity += mass_1_Acceleration * timestep
  mass_2_Velocity += mass_2_Acceleration * timestep
  
  ''  finally we find the change in position and add it to the existing position
  Mass_1_Position += mass_1_velocity * timestep
  Mass_2_Position += mass_2_velocity * timestep
  
  ''  keep track of time (just for fun. We don't need it for any calculations)
  time_passed += timestep
  
  ''  display the spring in numbers and graphics
  Screenlock
    
    ''  clear screen
    Cls
    
    ''  print various info
    Color RGB(64, 255, 64)
    Locate 2, 2,:  Print "green = constant"
    Locate 2, 55:  Print Using "  time step: #.### s"; timestep
    Color RGB(255, 255, 255)
    Locate 4, 2,:  Print "white = variable"
    Locate 4, 55:  Print Using "time passed: ###.# s"; time_passed
    
    ''  print spring data
    Color RGB(192, 192, 192)
    Locate 8, 26:  Print "------ SPRING DATA ------"
    Color RGB(64, 255, 64)
    Locate 10, 23: Print Using "     stiffnes: ######.## N/m"; spring_stiffnes
    Locate 12, 23: Print Using "  rest length: ######.## m"; Spring_rest_length
    Color RGB(255, 255, 255)
    Locate 14, 23: Print Using "spring length: ######.## m"; spring_length
    
    ''  print mass 1 data
    Color RGB(192, 192, 192)
    Locate 20, 7:  Print "------ MASS 1 DATA ------"
    Color RGB(64, 255, 64)
    Locate 22, 4:  Print Using "         mass: ######.## kg"; Mass_1_mass
    Color RGB(255, 255, 255)
    Locate 24, 4:  Print Using "        force: ######.## N"; -Force
    Locate 26, 4:  Print Using " acceleration: ######.## m/s^2"; Mass_1_Acceleration
    Locate 28, 4:  Print Using "     velocity: ######.## m/s"; Mass_1_Velocity
    Locate 30, 4:  Print Using "   ( position: ######.## )"; Mass_1_Position
    
    ''  print mass 2 data
    Color RGB(192, 192, 192)
    Locate 20, 45: Print "------ MASS 2 DATA ------"
    Color RGB(64, 255, 64)
    Locate 22, 42: Print Using "         mass: ######.## kg"; Mass_2_mass
    Color RGB(255, 255, 255)
    Locate 24, 42: Print Using "        force: ######.## N"; Force
    Locate 26, 42: Print Using " acceleration: ######.## m/s^2"; Mass_2_Acceleration
    Locate 28, 42: Print Using "     velocity: ######.## m/s"; Mass_2_Velocity
    Locate 30, 42: Print Using "   ( position: ######.## )"; Mass_2_Position
    
    ''  draw spring
    Line(Mass_1_Position, 320)-(Mass_2_Position, 320), RGB(255, 64, 64)
    Circle (Mass_1_Position, 320), Mass_1_Radius, RGB(192, 192, 192),,, 1, F
    Circle (Mass_2_Position, 320), Mass_2_Radius, RGB(192, 192, 192),,, 1, F
    
  Screenunlock
  
  ''  give the computer a break
  Sleep 1, 1
  
Loop Until Multikey(1)  ''  start all over with the two new position values

End
In the next part of the tutorial we'll go 2D and look further into the mysteries of friction and damping.

Cheers, Michael
Last edited by h4tt3n on Sep 12, 2010 8:58, edited 9 times in total.
notthecheatr
Posts: 1759
Joined: May 23, 2007 21:52
Location: Cut Bank, MT
Contact:

Post by notthecheatr »

Great... will this be in QBE as well?
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

Tbh the thought never crossed my mind until now ^^

But I suppose theres a good chance that I can polish this off and make it part 1 in a Freebasic mass-spring-damper tutorial series for QBE.

*edit*

Ok, until now I haven't realized how short QBE # 25 is on material. So, if you're there Pete I wouldn't mind compiling a FB spring tutorial Part 1 of XXX for your next issue. The subject is definitely too big for just one tutorial.

Cheers, Michael
notthecheatr
Posts: 1759
Joined: May 23, 2007 21:52
Location: Cut Bank, MT
Contact:

Post by notthecheatr »

I recommend you send it to qbexpress@gmail.com - it would especially help if you could format it properly in HTML using the QBE css. I don't know if they'll be able to get it into QBE 25 or not, but at least it'd be there for next issue. Pete seems to be unusually busy, though, so there might be time before he releases QBE 25.
roook_ph
Posts: 402
Joined: Apr 01, 2006 20:50
Location: philippines
Contact:

Post by roook_ph »

This is the best tutorial I have read so far spring could be easily modified to be like joints , bouncing ball, or anything moving or moved by force great job

you must teach more physics man!

sorry for late reply
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

Ok, here's finally part 2 of my spring tutorial. It will also appear in QBE # 27, which is due in a very short time. Enjoy.

1. INTRODUCTION

Hello and welcome to part 2 of my little spring tutorial!

In part 1 we went through all the underlying physics that we need to know about in order to create a realistic spring simulation, and we looked at a simple one dimensional code example.

In this next part we will go 2D and look into the mysteries of damping, which we need in order to make our springs behave in a more realistic way. In order to do all this we first need to know a little bit about vectors...

2. VECTORS

When you experimented with the code example from part 1 of this tutorial, you might have noticed how the force, acceleration and velocity values could be either positive or negative as the spring bounced back and forth. Surely there has to be something wrong here, because you can't drive at -20 mph! This is because these values arent "regular" numbers but in fact one dimensional vectors.

So, what's a vector? you might ask. Put very shortly, a vector is like a "regular" number, except it has both a value and a directon. A normal number, which we also call a scalar, only has a value. To stay with the driving analogy, the speed of a car is a scalar. If we want to say something about both the speed and direction of the car we can use a vector. Vectors are commonly depicted as arrows, like this:

Image

In 2D space a vector has a tail, usually at the coordianates (0, 0), and a head at the coordiates (x, y). The x and y values can be any negative or positive number, and they're called the horizontal component (x axis) and vertical component (y axis). In addition, any vector has a total length, also called the magnitude, which is defined by the distance between head and tail.

Thats a lot of things to have in mind, but all we need to do in order to use vectors in our programs is to store two numbers: the vector's length along the x and y axis, which define the position of the vector's head. Per definition the position of the tail is always (0, 0), so we don't need to store that. We don't always need to know the magnitude of the vector, so we don't calculate or store that either.

Now that we've been around the basics, let's take a look at how to use vectors in a program:


3. GOING FROM 1D TO 2D...

Once you understand the concept of vectors, going from 1D to 2D (or any dimension) is almost trivial. All that we do is use vectors instead of scalars when we calculate and store the vaules of force, acceleration, velocity and position.

Let's look at an example. This is a simple stripped-down simulation of masses and springs chained together, forming a rope-like structure:

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's FreeBasic mass-spring-damper tutorial, feb 2008
'
'   Example # 2.1
'
'   Simple two dimensional spring simulation. This program simulates an elastic 
'   rope made of serially connected springs and masses. (Press esc to quit)
'
'******************************************************************************'

''  set constants. experiment with these and see how the simulation reacts
Const Pi                =     4*Atn(1)          ''  pi (better not change ^^)
Const Timestep          =     0.01              ''  timestep
Const Num_Masses        =     5              -1 ''  number of masses in rope
Const Point_Mass        =     10                ''  mass of each point mass
Const Point_Density     =     0.01              ''  density of each point mass
Const Rope_Length       =     350               ''  rope length
Const Spring_Stiffnes   =     6000              ''  spring stiffnes
Const Grav_Acc          =     800               ''  gravitational acceleration
Const Rest_Length = Rope_Length/Num_Masses      ''  rest length of each spring

''  define types
''  two dimensional vector type
Type Vector_Type
  As Single X, Y
End Type

''  mass type including force-, acceleration-, velocity-, and position vectors
Type Mass_Type
  As Single Mass, Density, Radius
  As Vector_Type Frc, Acc, Vel, Pos
End Type 

''  dimension variables
Dim As Vector_Type Length, Normalised_Length, Vel
Dim As Mass_Type Mass(Num_Masses)
Dim As Single Force, Spring_Length
Dim As Integer i, i2, Scrn_Wid, Scrn_Hgt

''  set screen width, height, and bit depth
Scrn_Wid = 800
Scrn_Hgt = 600
ScreenRes Scrn_Wid, Scrn_hgt, 16

''  set startup conditions of masses
For i = 0 To Num_Masses
  
  With Mass(i)
    
    .Mass = Point_Mass
    .Density = Point_Density
    .Radius = ((.Mass/.Density)/((4/3)*pi))^(1/3)
    .Pos.x = (Scrn_Wid\2)+(i*Rest_Length)
    .pos.y = 100
    
  End With
  
Next

''  main program loop
Do
  
  ''  calculate the spring forces acting upon each point mass in the rope
  For i = 0 To Num_Masses-1
    
    i2 = i+1
    
    ''  spring length (Pythagoras' theorem)
    ''  (we need this to find the magnitude of the spring force)
    Length.X = Mass(i2).Pos.X-Mass(i).Pos.X
    Length.Y = Mass(i2).Pos.Y-Mass(i).Pos.Y
    Spring_Length = Sqr(Length.X*Length.X+Length.Y*Length.Y)
    
    ''  normalise the spring length vector
    ''  (same as finding sine and cosine to the angle between masses)
    ''  (we need this to find the direction of the spring force)
    Normalised_Length.X = Length.X/Spring_Length
    Normalised_Length.Y = Length.Y/Spring_Length
    
    ''  spring force (Hooke's law of elasticity)
    Force = -Spring_Stiffnes*(Spring_Length-Rest_Length)
    
    ''  split spring force into horizontal and vertical vector component
    ''  (each action has an equal opposite reaction, hence the - and + sign)
    Mass(i).Frc.X -= Force*Normalised_Length.X
    Mass(i).Frc.Y -= Force*Normalised_Length.Y
    
    Mass(i2).Frc.X += Force*Normalised_Length.X
    Mass(i2).Frc.Y += Force*Normalised_Length.Y
    
  Next
  
  ''  update acceleration, velocity, and position of point masses
  ''  (except for mass 0, which we'd like to keep fixed)
  For i = 1 To Num_Masses
    
    With Mass(i)
      
      ''  accelerate masses:
      ''  acceleration = force / mass
      .Acc.X = .Frc.X/.mass
      .Acc.Y = .Frc.Y/.mass
      
      ''  add gravity (downwards acceleration)
      .Acc.Y += Grav_Acc
      
      ''  update velocity:
      ''  delta velocity = acceleration * delta time
      ''  new velocity = old velocity + delta velocity
      .Vel.X += .Acc.X*Timestep
      .Vel.Y += .Acc.Y*Timestep
      
      ''  update position:
      ''  delta position = velocity * delta time
      ''  new position = old position + delta position
      .Pos.X += .Vel.X*Timestep
      .Pos.Y += .Vel.Y*Timestep
      
      ''  reset force vector
      .frc.x = 0
      .frc.y = 0
      
    End With
    
  Next
  
  ''  display the rope in graphics
  ScreenLock
    
    ''  clear screen
    Cls
    
    ''  draw springs
    For i = 0 To Num_Masses-1
      i2 = i+1
      Line (Mass(i).Pos.X, Mass(i).Pos.Y)-(Mass(i2).Pos.X, Mass(i2).Pos.Y), Rgb(64, 255, 64)
    Next i
    
    ''  draw masses
    For i = 0 To Num_Masses
      With Mass(i)
        Circle(.Pos.X, .Pos.Y), .Radius, Rgb(255, 64, 64),,, 1, f
      End With
    Next i
    
  ScreenUnlock
  
  Sleep 1, 1
  
Loop Until Multikey(1)

End
Essentially this code works exactly like the the example in part 1 of the tutorial, but the 2D environment presents some new obstacles, that we need to overcome:

Let's start out with an easy one. In order to calculate spring forces we need to know the distance between the two endpoint masses of the spring, which effectively is the same as the spring's length. In order to do that we use Pythagoras' theorem, which in short form states that for any right triangle the square on the hypotenuse is equal to the sum of the squares on the other two sides:

Image

This also means that the hypotenuse - or the spring's length - is equal to the square root of the sum of the two other sides squared. This is how it might look in code:

Lng.X = Mass(2).Pos.X-Mass(1).Pos.X
Lng.Y = Mass(2).Pos.Y-Mass(1).Pos.Y
Spring_Length = Sqr(Lng.X*Lng.X+Lng.Y*Lng.Y)

Now lets move on to a slightly harder obstacle. We only keep track of the forces acting upon each mass along the x and y axis, but the springs may be rotated in any non-vertical or non-horizontal angle. This means that the endpoint masses of the springs may be influeced by forces acting on them in any possible angle.

How do we split these angled forces up into vertical and horizontal force components, which the simulation can handle?

Once again we need to look at vectors and vector math, specifically on how to normalise a vector. Normalising a vector means changing its magnitude to 1 without changing its direction. When a vector has been normalised it is called a unit vector or "hat" vector, because it is written like this... Â, with a little hat on top. Normalising a vector is done by dividing all of its components by its magnitude. For 2d vector A it looks like this:

 = A/|A|

or if we look at each vector component separately

Â.x = A.x/|A|
Â.y = A.y/|A|

where x and y are the vector components, A is the vector, Â is the normalised vector, and |A| is the magnitude of the vector. In code form it might look something like this:

Vector_magnitude = SQR(Length.X^2+Length.Y^2)
Normalized_Vector.X = Length.X/Vector_magnitude
Normalized_Vector.Y = Length.Y/Vector_magnitude

A normalised vector is so to speak stripped of any information about magnitude, it only tells us direction. Now, since the force applied by a spring is always along the spring axis - either away from or towards the spring's centre - all we need to do is normalise the spring (as if it were a vector) in order to get the force direction along the x and y axis. Thats what happens in this part of the code example:

Length.X = Mass(i2).Pos.X-Mass(i).Pos.X
Length.Y = Mass(i2).Pos.Y-Mass(i).Pos.Y
Spring_Length = SQR(Length.X^2+Length.Y^2)

Normalised_Length.X = Length.X/Spring_Length
Normalised_Length.Y = Length.Y/Spring_Length

Then we multiply the spring force with the normalized x and y vector components to get the magnitude of the force vector along the x and y axis:

Mass(i).Frc.X -= Force*Normalised_Length.X
Mass(i).Frc.Y -= Force*Normalised_Length.Y

Mass(i2).Frc.X += Force*Normalised_Length.X
Mass(i2).Frc.Y += Force*Normalised_Length.Y

Remember Newton's 3rd law from part 1? It says that for each action - or force - there is an opposite and equal reaction. This is why it says -= for one mass and += for the other.

Ok, now we got that one solved too. Before moving on to the final part, here's another code example showing a somewhat different implementation of spring physics. In this case the example simulates a number of elastic balls, which jump around and interact. Spring forces are only applied when the balls touch, making them bounce away from each other.

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's FreeBasic mass-spring-damper tutorial, feb 2008
'
'   Example # 2.2
'
'   Simple two dimensional spring simulation. This program simulates a number
'   of elasic balls which jump around and interact. (Press esc to quit.)
'
'******************************************************************************'

''  set constants. experiment with these and see how the simulation reacts.
Const Pi                =     4*Atn(1)      ''  pi (better not change ^^)
Const Timestep          =     0.01          ''  timestep (delta time)
Const Num_Balls         =     12         -1 ''  number of elastic balls
Const Mass_Min          =     400           ''  smallest elastic ball mass
Const Mass_Max          =     8000         ''  biggest elastic ball mass
Const Density           =     0.01          ''  elastic ball density 
Const Spring_Stiffnes   =     5e5           ''  elastic stiffnes

''  define types
Type Vector_Type
  As Single X, Y
End Type

Type Ball_Type
  As Uinteger Col
  As Single Mass, Density, Radius
  As Vector_Type Frc, Acc, Vel, Pos
End Type 

''  dimension variables
Dim As Vector_Type Dst, Vel
Dim As Ball_Type Ball(Num_Balls)
Dim As Single Force, Ball_Distance, Cos_Angle, Sin_Angle, Dist_Min
Dim As Integer i, i2, Scrn_Wid, Scrn_Hgt

''  set screen width, height, and bit depth
Scrn_Wid = 800
Scrn_Hgt = 600
ScreenRes Scrn_Wid, Scrn_hgt, 16

''  use timer to generate random numbers
Randomize Timer

''  set startup conditions of elastic balls
For i = Lbound(Ball) To Ubound(Ball)
  
  With Ball(i)
    
    .Col = RGB(64+Rnd*192, 64+Rnd*192, 64+Rnd*192)
    .Mass = Mass_Min+(i/Ubound(Ball))*(Mass_Max-Mass_Min)
    .Density = Density
    .Radius = ((.Mass/.Density)/((4/3)*pi))^(1/3)
    .Pos.x = .Radius+Rnd*(Scrn_Wid-.Radius)
    .pos.y = .Radius+Rnd*(scrn_hgt-.Radius)
    .vel.x = (Rnd-Rnd)*200
    .vel.y = (Rnd-Rnd)*200
    
  End With
  
Next

''  main program loop
Do
  
  ''  test all elastic balls against each other. calculate forces if they touch.
  For i = Lbound(Ball) To Ubound(Ball)-1
    
    For i2 = i+1 To Ubound(Ball)
      
      ''  distance between elastic balls (Pythagoras' theorem)
      Dst.X = Ball(i2).Pos.X-Ball(i).Pos.X
      Dst.Y = Ball(i2).Pos.Y-Ball(i).Pos.Y
      Ball_Distance = Sqr(Dst.X*Dst.X+Dst.Y*Dst.Y)
      Dist_Min = Ball(i).Radius+Ball(i2).Radius
      
      If Ball_Distance < Dist_Min Then
        
        ''  cosine and sine to the angle between Ball i and i2 (trigonometry)
        Cos_Angle = Dst.X/Ball_Distance
        Sin_Angle = Dst.Y/Ball_Distance
        
        ''  spring force (Hooke's law of elasticity)
        Force = -Spring_Stiffnes*(Ball_Distance-Dist_Min)
        
        ''  split spring force vector into horizontal and vertical component
        Ball(i).Frc.X -= Force*Cos_Angle
        Ball(i).Frc.Y -= Force*Sin_Angle
        Ball(i2).Frc.X += Force*Cos_Angle
        Ball(i2).Frc.Y += Force*Sin_Angle
        
      End If
      
    Next
    
  Next
  
  ''  update acceleration, velocity, and position of elastic balls
  ''  (using the Euler-Cromer 1st order integration algorithm)
  For i = Lbound(Ball) To Ubound(Ball)
    
    With Ball(i)
      
      ''  accelerate balls (acceleration = force / mass)
      .Acc.X = .Frc.X/.Mass
      .Acc.Y = .Frc.Y/.Mass
      
      ''  reset force vector
      .frc.x = 0
      .frc.y = 0
      
      ''  update velocity
      ''  delta velocity = acceleration * delta time
      ''  new velocity = old velocity + delta velocity
      .Vel.X += .Acc.X*Timestep
      .Vel.Y += .Acc.Y*Timestep
      
      ''  update position 
      ''  delta position = velocity * delta time
      ''  new position = old position + delta position
      .Pos.X += .Vel.X*Timestep
      .Pos.Y += .Vel.Y*Timestep
      
    End With
    
  Next
  
  ''  keep elastic balls within screen boundaries
  For i = Lbound(Ball) To Ubound(Ball)
    
    With Ball(i)
      
      ''  right
      If .Pos.x > Scrn_Wid-.Radius Then
        .Vel.X = -.vel.x
        .pos.x = Scrn_Wid-.Radius
      End If
      ''  left
      If .Pos.x < .Radius Then
        .Vel.X = -.vel.x
        .pos.x = .Radius
      End If
      ''  bottom
      If .Pos.Y > Scrn_Hgt-.Radius Then
        .vel.y = -.vel.y
        .pos.y = Scrn_Hgt-.Radius
      End If
      ''  top
      If .Pos.Y < .Radius Then
        .vel.y = -.vel.y
        .pos.y = .Radius
      End If
      
    End With
    
  Next
  
  ''  display elastic balls in graphics
  ScreenLock
    
    ''  clear screen
    Cls
    
    ''  draw elastic balls
    For i = Lbound(Ball) To Ubound(Ball)
      With Ball(i)
        Circle(.Pos.X, .Pos.Y), .Radius    , .Col,,, 1
        Circle(.Pos.X, .Pos.Y), .Radius-0.5, .Col,,, 1
        Circle(.Pos.X, .Pos.Y), .Radius-1.0, .Col,,, 1
      End With
    Next i
    
  ScreenUnlock
  
  Sleep 1, 1
  
Loop Until Multikey(1)

End

4. DAMPING

In the beginning of part 1 of this tutorial I mentioned that the correct technical term is mass-spring-damper system. We've been around masses and springs, but what about dampers?

In the previous code examples you've probably noticed how the springs and rubberballs keep on bouncing - apparently forever - once they've been prodded into movement. We don't see this behaviour in real life, and damping is simply a way of reducing this effect in a somewhat physically plausible way.

Damping a simulated spring works a little like the damper inside a car spring. Without damping the car would bounce up and down for quite a while every time it hit a bump on the road. On the other hand a car with oversized, hard dampers would be very rigid, because they counter the effect of the spring. If you get springs and dampers adjusted just right, you've got a spring that serves your need and at the same time stops bouncing very quickly - like a real world spring.

So, how do we simulate spring damping? Again we need to look at some vector math. First we need to look at the so called dot product of two vectors. Vectors are different from scalars in the way that they can be multiplied by each other in several different ways. The one we are interested in is called the dot product or scalar product, because the result is expressed as a "regular" scalar number.

Let's look at scalar s and 2d vectors A and B. This is how for example "s equals A dot B" is written on paper:

s = A . B (except the dot is a bit higher up)

And this is what it means:

s = A.x * B.x + A.y * B.y

Or in plane words... the product of the horizontal components is added to the product of the vertical components. The useful thing about the dot product is that

A . B = |A| * |B| * cos(v)

where A and B are 2d vectors, |A| and |B| is the magnitude of A and B, and v is the angle between the two vectors.

Next we need to look at how to do a so called scalar projection of one vector onto another vector, as shown in this illustration:

Image

Basically what we do is find out how long a line we would get if we flattened vector A out on vector B, or if we "cast the shadow" of A at a right angle onto B. The scalar projection of A onto B is equal to

|A| * cos(v)

where |A| is the magnitude of A, and v is the angle between vector A and B.

Since we already know that

A . B = |A| * |B| * cos(v)

then we can see that

A . B / |B| = |A| * cos(v)

Or in other words that the scalar projection of A onto B is equal to A dot B divided by the magnitude of B. Now look at this image where we find the scalar projection of the velocity vector on the spring vector:

Image

Notice the arrow called damping force? It is simply the negative value of the scalar projection of the velocity vector onto the spring. And it is the heart of our damping algorithm! This force has two very useful qualities:

-It always has the right direction. It always pushes the endpoint mass in the opposite direction of its velocity, thus slowing it down.

-It always has the right size. If the velocity vector is perpendicular to the spring, then the scalar projection and therefore the damping force is very small. If the velocity vector is almost parallel to the spring, then the scalar projection and damping force is big.

The damping force - which is a scalar - is subtracted from the spring force before this is split up into a horizontal and vertical component, so it will always work along the spring axis and never dampen the "sideways" or perpendicular movement of the spring masses.

Voila! And it actualy works. Take a look at this code sample. It is exactly like the first one, except damping has been added:

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's FreeBasic mass-spring-damper tutorial, feb 2008
'
'   Example # 2.3
'
'   Simple two dimensional spring simulation. This program simulates an 
'   elastic rope made of serially connected springs and point masses.
'   Damping is added, which makes the rope look more realistic. 
'   (Press esc to quit)
'
'******************************************************************************'

''  set constants. experiment with these and see how the simulation reacts
Const Pi                =     4*Atn(1)          ''  pi (better not change ^^)
Const Timestep          =     0.01              ''  timestep
Const Num_Masses        =     5             -1 ''  number of masses in rope
Const Point_Mass        =     10                 ''  mass of each point mass
Const Point_Density     =     0.01              ''  density of each point mass
Const Rope_Length       =     350               ''  rope length
Const Spring_Stiffnes   =     6000             ''  spring stiffnes
Const Spring_Damping    =     2e2              ''  spring damping
Const Grav_Acc          =     800              ''  gravitational acceleration
Const Rest_Length = Rope_Length/Num_Masses      ''  rest length of each spring

''  define types
Type Vector_Type
  As Single X, Y
End Type

Type Mass_Type
  As Single Mass, Density, Radius
  As Vector_Type Frc, Acc, Vel, Pos
End Type 

''  dimension variables
Dim As Vector_Type Lng, Norm_Lng, Vel
Dim As Mass_Type Mass(Num_Masses)
Dim As Single Force, Spring_Length
Dim As Integer i, i2, Scrn_Wid, Scrn_Hgt

''  set screen width, height, and bit depth
Scrn_Wid = 800
Scrn_Hgt = 600
ScreenRes Scrn_Wid, Scrn_hgt, 16

''  set startup conditions of masses
For i = Lbound(Mass) To Ubound(Mass)
  With Mass(i)
    .Mass = Point_Mass
    .Density = Point_Density
    .Radius = ((.Mass/.Density)/((4/3)*pi))^(1/3)
    .Pos.x = (Scrn_Wid\2)+(i*Rest_Length)
    .pos.y = 100
  End With
Next

''  main program loop
Do
  
  ''  calculate the forces acting upon each point mass in the rope
  For i = Lbound(Mass) To Ubound(Mass)-1
    
    i2 = i+1
    
    ''  spring length (Pythagoras' theorem)
    Lng.X = Mass(i2).Pos.X-Mass(i).Pos.X
    Lng.Y = Mass(i2).Pos.Y-Mass(i).Pos.Y
    Spring_Length = Sqr(Lng.X*Lng.X+Lng.Y*Lng.Y)
    
    ''  normalise the distance vector 
    ''  (same as finding cosine and sine to the angle)
    Norm_Lng.X = Lng.X/Spring_Length
    Norm_Lng.Y = Lng.Y/Spring_Length
    
    ''  relative velocity (needed for damping, see below)
    Vel.X = Mass(i2).Vel.X-Mass(i).Vel.X
    Vel.Y = Mass(i2).Vel.Y-Mass(i).Vel.Y
    
    ''  spring force (Hooke's law of elasticity)
    Force = -Spring_Stiffnes*(Spring_Length-Rest_Length)
    
    ''  spring damping 
    ''  (scalar projection of velocity vector onto spring)
    Force -= Spring_Damping*((Lng.X*Vel.X+Lng.Y*Vel.Y)/Spring_Length)
    
    ''  split spring force vector into horizontal and vertical component
    ''  (each mass moves in opposite directions, hence the - and + sign)
    Mass(i).Frc.X -= Force*Norm_Lng.X
    Mass(i).Frc.Y -= Force*Norm_Lng.Y
    Mass(i2).Frc.X += Force*Norm_Lng.X
    Mass(i2).Frc.Y += Force*Norm_Lng.Y
    
  Next
  
  ''  update acceleration, velocity, and position of point masses
  ''  (using the Euler-Cromer 1st order integration algorithm)
  For i = Lbound(Mass)+1 To Ubound(Mass)
    
    With Mass(i)
      
      ''  accelerate masses (acceleration = force / mass)
      .Acc.X = .Frc.X/.mass
      .Acc.Y = .Frc.Y/.mass
      
      ''  add gravity (downwards acceleration)
      .Acc.Y += Grav_Acc
      
      ''  update velocity
      ''  delta velocity = acceleration * delta time
      ''  new velocity = old velocity + delta velocity
      .Vel.X += .Acc.X*Timestep
      .Vel.Y += .Acc.Y*Timestep
      
      ''  update position 
      ''  delta position = velocity * delta time
      ''  new position = old position + delta position
      .Pos.X += .Vel.X*Timestep
      .Pos.Y += .Vel.Y*Timestep
      
      ''  reset force vector
      .frc.x = 0
      .frc.y = 0
      
    End With
    
  Next
  
  ''  display the rope in graphics
  ScreenLock
    
    ''  clear screen
    Cls
    
    ''  draw springs
    For i = Lbound(Mass) To Ubound(Mass)-1
      i2 = i+1
      Line (Mass(i).Pos.X, Mass(i).Pos.Y)-(Mass(i2).Pos.X, Mass(i2).Pos.Y), Rgb(64, 255, 64)
    Next i
    
    ''  draw point masses
    For i = Lbound(Mass) To Ubound(Mass)
      With Mass(i)
        Circle(.Pos.X, .Pos.Y), .Radius, Rgb(255, 64, 64),,, 1, f
      End With
    Next i
    
  ScreenUnlock
  
  Sleep 1, 1
  
Loop Until Multikey(1)

End
That's it, hope you had fun. In part 3 we will take a closer look at friction forces and go 3D. Homework until then is to add damping to the elastic ball code sample and experiment with it ^^

Cheers, Michael
Last edited by h4tt3n on Jan 26, 2009 11:20, edited 2 times in total.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Post by duke4e »

WOW, GREAT TUTORIAL!!!!
AWESOME AS PREVIOUS ONE!!!

Keep up the good work!
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

Michael "h4tt3n" Nissen's FreeBasic spring simulation tutorial
part 3 - angular springs, also appearing in QBExpress # 30...

1. INTRODUCTION

Hello everyone, and welcome to part three of my little spring tutorial. In part one we looked at all the basics needed to create a physics simulation, and we looked at a simple one dimensional code example of a mass-spring simulation. In part two we had a good look at the vector math needed in 2d spring physics, and finally we looked at a method for doing velocity dependent damping.

In this part we will stay in 2d and learn how to implement angular springs. In the previous two parts of this tutorial we've been working with what I prefer to call linear springs, which act on two point masses and tries to keep them at a fixed distance by applying a ”linear force” along the spring axis. Similarly, the type of angular spring we will be looking at acts on two linear springs and tries to keep them at a fixed angle by applying an ”angular force” - also called torque - perpendicular to the spring axis.

Image

Angular springs are a very useful tool for building many different types of deformable or flexible ”soft bodies” out of linear springs and point masses, like for instance hair, trees, water balloons, cloth, steel girders, and jelly. I'll cover this subject in greater detail in a later tutorial. Put very shortly, in this part we'll learn:

1.how to find the angle between two connected linear springs
2.how to apply an angular force to the springs in order to keep them at a fixed angle

Let's get started by looking at...


2. HOW TO FIND THE ANGLE BETWEEN TWO SPRINGS

In the previous part of this tutorial we learned a great deal of vector math, especially in the part about velocity dependent damping. One of the very important equations were about finding the scalar product or ”dot product” of two vectors:

A · B = |A| * |B| * Cos(v)

where A and B are 2d vectors, |A| and |B| is the magnitude of A and B, and v is the angle between the two vectors. If we rearrange the equation a bit we get

Cos(v) = A · B / |A| * |B|

Which can be boiled further down to just taking the dot product of the normalized vectors:

Cos(v) = (A / |A|) · (B / |B|) = A normalized · B normalized

And voila! Now we already know how to get cosine to the angle between the springs, even without calling any CPU-expensive trig functions. We've already found the magnitude of the springs and normalised them in order to calculate the linear spring force, so we just recycle that information here. This helps making the angular spring implementation computationally very cheap. But what about sine to the angle? Lucky for us there's a very similar function just for that. Remember the definition of the dot product:

A · B = A.x * B.x + A.y * B.y

There is a similar function called the perpendicular dot product, also called ”perproduct” or ”perpdot” for short, which looks like this:

A ”perpdot” B = A.x * B.y – A.y * B.x

If we use perpdot instead of dot in the above equation we get:

Sin(v) = A ”perpdot” B / |A| * |B|

Or boiled down like the cosine equation:

Sin(v) = A normalised ”perpdot” B normalised

Image

There we go. Now we know how to find sine and cosine of the angle between two springs. If we wanted to know the actual angle in radians we could use the ATan2() function, but as we shall see later there's really no need to. Here you can see what all the above equations look like when they're implemented in a working code sample:

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's mass-spring-damper tutorial, November 2008
'
'   Example # 3.1
'
'   Simple angle test program. In this example we use vector math to find sine
'		and cosine of the angle between the two springs. (Press esc to quit)
'
'******************************************************************************'

''	define constants
const pi            = 4*atn(1)		''	pi
Const rad2deg 			= 180/pi			''	radians to degrees
Const num_masses 		= 3						''	number of point masses
Const Rope_Length		=	350					''	total length of all springs
Const Scrn_Wid      = 600					''  screen width
Const Scrn_Hgt      = 600					''  screen height

''  define types
Type Vector_Type
  As Double X, Y
End Type

Type Mass_Type
  As Double Mass, dens, radius
  As Vector_Type Frc, Acc, Vel, Psn
End Type 

Type Spring_Type
  As Integer a, b
  As Double Length
  As Vector_Type Nor_Lng, Lng
End Type

Type Angular_Spring_Type
  As Integer a, b
  As Double sin_angle, cos_angle
End Type

''	dim arrays & variables
Dim As mass_type mass(1 to num_masses)
dim as spring_type spring(1 to num_masses-1)
dim as angular_spring_type angular_spring(1 to num_masses-2)
Dim As Integer mouse_x, mouse_y

''	initiate point masses
for i As Integer = lbound(mass) to Ubound(mass)
  with mass(i)
    .radius = 5
    .Psn.y = scrn_hgt*0.5-Rope_Length*0.5+(Rope_Length/(num_masses-1))*(i-1)
    .Psn.x = scrn_wid*0.5
  end with  
Next

''	initiate linear springs
''	(a and b stores the number of the masses it's attached to)
for i As Integer = lbound(spring) to ubound(spring)
  with spring(i)
    .a = i
    .b = i+1
  end with
next

''	initiate angular springs
''	(a and b stores the number of the linear springs it's attached to)
for i As Integer = lbound(angular_spring) to ubound(angular_spring)
  with angular_spring(i)
    .a = i
    .b = i+1
  end with
Next

''	initiate screen
ScreenRes Scrn_Wid, Scrn_hgt, 16
Color RGB(0, 0, 0), RGB(255, 255, 255)

Do
	
	''	use mouse to change spring angle
	GetMouse Mouse_x, Mouse_Y
	Mass(1).Psn.X = Mouse_X
	Mass(1).Psn.Y = Mouse_Y
	
	''	get linear spring data
	For i As Single = Lbound(Spring) to Ubound(Spring)
		With Spring(i)
			''	spring length vector
			.Lng.X = Mass(.b).Psn.X-Mass(.a).Psn.X
			.Lng.Y = Mass(.b).Psn.Y-Mass(.a).Psn.Y
			''	spring length scalar
			.Length = Sqr(.Lng.X*.Lng.X+.Lng.Y*.Lng.Y)
			''	normalised spring length vector
			.Nor_Lng.X = .Lng.X/.Length
			.Nor_Lng.Y = .Lng.Y/.Length
		End With
	Next
	
	''	get angular spring data
	For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
		With Angular_Spring(i)
			''	cosine of the angle - dot product of normalized spring vectors
			.cos_angle = Spring(.a).Nor_Lng.x*Spring(.b).Nor_Lng.x+Spring(.a).Nor_Lng.y*Spring(.b).Nor_Lng.y
			''	sine of the angle - perpendicular dot product of normalized spring vectors
			.sin_angle = Spring(.a).Nor_Lng.x*Spring(.b).Nor_Lng.y-Spring(.a).Nor_Lng.y*Spring(.b).Nor_Lng.x
		End With
	Next
	
	''	draw to screen
	ScreenLock
	
		''	clear screen
		Cls 
		
		''	print data
		With Angular_Spring(1)
			Locate 2, 2: Print Using "   Sine of spring angle: ####.### (-1, 1)"; .sin_Angle
			Locate 4, 2: Print Using " Cosine of spring angle: ####.### (-1, 1)"; .cos_Angle
			Locate 6, 2: Print Using "Spring angle in radians: ####.### (-pi, pi)"; ATan2(.Sin_Angle, .Cos_Angle)
			Locate 8, 2: Print Using "Spring angle in degrees: ####.### (-180, 180)"; ATan2(.Sin_Angle, .Cos_Angle)*rad2Deg
		End With
		
		''  draw springs
		For i As Integer = Lbound(Spring) To Ubound(Spring)
			With Spring(i)
				Line (Mass(.a).Psn.X, Mass(.a).Psn.Y)-(Mass(.b).Psn.X, Mass(.b).Psn.Y), Rgb(64, 160, 32)
			End With
		Next
		
		''  draw particles
		For i As Integer = Lbound(mass) To Ubound(mass)
			With mass(i)
				Circle(.Psn.X, .Psn.Y), .radius, rgb(0, 0, 0),,, 1, f
				Circle(.Psn.X, .Psn.Y), .radius-1, rgb(255, 32, 32),,, 1, f
			End With
		Next
		
	ScreenUnLock
	
	Sleep 1, 1
	
Loop Until MultiKey(1)

End
3. A LITTLE MATH AND PHYSICS ABOUT ROTATION

In the previous parts of this tutorial we learned how to apply a linear force along the spring axis. When dealing with angular springs we don't want to apply forces parallel to the spring axis but rather perpendicular to it. We also need to find a good angular spring equation that applies the right amount of ”angular force” in the right situation. To do this we first need to get familiar with two new concepts: normal vectors and torque.

A normal vector - not to be confused with a normalised vector! - is simply a different name for a perpendicular vector. Finding the normal of a vector is very easy - the normal of vector (x, y) is simply (y, -x). This vector always points to the right of the vector it is normal to and is hence called a right-hand normal vector. There also is a left-hand normal (-y, x), but the right-hand normal is most commonly used, so we'll stick to that. If we find the normal of a spring vector and then normalise it – in other words the normalised normal vector - and use this to define the angular spring force direction, it will always be perpendicular to the spring's axis, which is just what we want.

So much for the direction of the angular force vector, but what about its magnitude – how much force do we apply? This is where torque comes into the picture. When dealing with physics, torque is simply the rotational equivalent of force. When you apply force to an object it moves along a straight line, and when you apply torque to an object it rotates. Within the boundaries of our simple 2d simulation, the definition of torque is

T = r * F

Which also means that

F = T / r

Where T is torque, r is distance to center of rotation, and F is force. Please note that this definition is only valid if the force vector is at a right angle to the distance vector, but since this is always the case in our simulation we'll just keep it at that and skip the more complex definitions of torque for now.

Image

So, why do we need to know about torque? Let's go back and have a look at Newton's 3rd law, stating that for every force (action) there is an equal and opposite force (reaction). For linear springs this means that when applying spring force F to one of the endpoint particles, we have to apply exactly -F to the other particle, otherwise the simulation won't behave physically correct. We say that the net force on the system must be zero. This law is also valid for torque, which means that the torque we apply on one of the linear springs must be ”canceled out” by the torque we apply on the other - the net torque on the system must be zero. Please note that in this case force doesn't have to zero out, but force times distance to rotational center – or torque - does! And this means trouble when the springs are of different length, which they are more or less all the time. We somehow need to scale the forces, so they result in equal, opposite torques.

Besides all this there is another problem with the equation F = T / r, namely that force increases towards infinity as spring length decreases towards zero and vice versa. Practically speaking this means that two identical angular springs behave differently depending on wether they're applied to long or short linear springs. If one or both of the linear springs has near-zero length, then the applied force will be near-infinite, which usually results in simulation instability or even program crash. So we need to scale the forces, so they're always proportional to the length of the springs they're applied to.

Long story short, I've found a very simple and stable solution to both of the above mentioned problems. The applied angular force on linear spring 1 and 2 is simply

Angular Force 1 = Spring Force * Angular Spring 2 Length
Angular Force 2 = -Spring Force * Angular Spring 1 Length

Ok, now we've been around normals and torque. Next thing is to find a suitable angular spring force equation, and then we're ready to rock. Let's look at a practical example: We'll use the damped rope simulation (example 2.3) from the previous part of this tutorial and apply an angular spring between each pair of linear springs in the rope. Let's say we want to make all the springs stay parallel and form a sort of ”soft rod”, essentially behaving like a hair or a flag pole. In order to do this we first need an angular spring force equation. Remember the linear spring force equation, Hooke's law of elasticity?

Linear Force = -Stiffnes * (Length – Rest length)

This equation always produces a scalar spring force value of the right magnitude and with the right sign (+ or -). What we need to do is write a similar equation for angular springs. It's really tempting and intuitive to assume that we'll be all right by just translating the linear spring equation into

Angular Force = -Stiffnes * (Angle – Rest angle)

I've tested a few implementations of this method and discarded them for several reasons. One of the main reasons is that you need to call the cpu-expensive ATan2() function in order to get the angle between the springs. Besides, when bending the angular spring past the ”180 degree mark” away from it's rest angle, the spring force – which is at maximum - suddenly changes direction, often resulting in chaotic movement and simulation instability. Of course we could implement it if we really wanted to, but I just don't think it's worth the effort.

Instead, let me introduce you to a more stable and computationally much cheaper method. In our case with the ”soft rod” simulation all we need is sine of the angle between the springs in order to calculate the angular force: If the angle between two lines is 0 or 180 degrees – in other words they're parallel – then sine of the angle is zero. If for instance the angle is a little bigger than – say - 0 degrees, then sine of the angle is also slightly higher than zero and vice versa. This is a really good foundation for a spring equation. In our case the angular spring force equation is simply

Force = -Stiffnes * Sine angle

This equation will produce a force value that increases as the angle diverges more and more from 0 and which is signed either positive or negative depending on what direction the spring is bent. In other words we have an angular spring with a rest angle of 0 degrees.

Image

If the angular spring is bent +/-90 degrees away from it's rest angle – which means sine to the angle is 1 – the spring force reaches its maximum. If the spring is bent even further the spring force decreases and reaches 0 when the spring is bent +/-180 degrees, where sine becomes 0 again. This sounds like a drawback, since the angular spring force will be 0 when the spring is bent furthest away from its rest angle, but practically it works quite well, and it avoids the above mentioned ”force-suddently-changes-direction” problem in the angle based equation.

If we used cosine instead of sine in our equation we'd get an angular spring that would be at rest at +/- 90 degrees. Combining sine and cosine in one spring equation allows us to pick several different rest angles. Here's an equation template ranging from -180 to +180 degrees at 45 degree intervals:

+/-180 degrees: Force = -Angular Stiffnes*(-Sin Angle)
+135 degrees: Force = -Angular Stiffnes*(-Sin Angle-Cos Angle)
+90 degrees: Force = -Angular Stiffnes*(-Cos Angle)
+45 degrees: Force = -Angular Stiffnes*(Sin Angle-Cos Angle)
0 degrees: Force = -Angular Stiffnes*(Sin Angle)
-45 degrees: Force = -Angular Stiffnes*(Sin Angle+Cos Angle)
-90 degrees: Force = -Angular Stiffnes*(Cos Angle)
-135 degrees: Force = -Angular Stiffnes*(-Sin Angle+Cos Angle)

If we want to simulate a soft body composed of simple geometric shapes like rods, boxes, triangles, octagons and the like, these equations will be adequate. Let's look at a code sample and see if all this works, shall we?

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's mass-spring-damper tutorial, November 2008
'
'   Example # 3.2
'
'		Simple angular spring simulation. This code sample is showing how linear 
'		and angular springs can be connected into a "soft rod" of different shapes. 
'		One of the point masses has been glued stuck and acts as a pivot for the
'		rest of the body. Scroll down a bit and test the different angular spring 
'		force equations. 
'
'  	(Press esc to quit)
'
'******************************************************************************'

''	define constants. Play around with these and see what happens!
const pi            		=	4*atn(1)		''	pi (better not change ^^)
Const rad2deg 					= 180/pi			''	radians to degrees
Const deg2rad						= Pi/180			''	degrees to radians
Const Grav_acc		  		= 2000				''	gravitational acceleration
const num_masses    		= 8						''	number of masses in rod
Const Rod_Length   			= 300					''	total spring body length
Const Linear_Stiffnes  	= 10000       ''  linear spring stiffnes
Const Linear_Damping    = 40         	''  linear spring damping
Const Angular_Stiffnes	= 10000      	''  angular spring stiffnes
Const Scrn_Wid      		= 600         ''  screen width
Const Scrn_Hgt      		= 600         ''  screen height
Const dt            		= 0.0025			''  timestep, delta time
Const	air_friction			= 0.002				''	air friction (to calm springs down)

''  define types
Type Vector_Type
  As Double X, Y
End Type

Type Mass_Type
  As Double Mass, radius
  As Vector_Type Frc, Acc, Vel, Psn
End Type 

Type Linear_Spring_Type
  As Integer a, b
  As Double Length, Rest_Length
  As Vector_Type Nor_Lng, Vel, Lng
End Type

Type Angular_Spring_Type
  As Integer a, b
  As Double sin_angle, cos_angle
End Type

''	dim variables
Dim As mass_type mass(1 to num_masses)
dim as Linear_spring_type Linear_spring(1 to num_masses-1)
dim as angular_spring_type angular_spring(1 to num_masses-2)
Dim As Integer mouse_x, mouse_y

Randomize timer
 
''	initiate point masses
for i As Integer = lbound(mass) to Ubound(mass)
  with mass(i)
    .mass = 1
    .radius = 4
    .Psn.x = scrn_wid*0.5+(Rod_Length/(num_masses-1))*(i-2)
    .Psn.y = scrn_hgt*0.5
  end with  
Next

''	initiate linear springs
for i As Integer = lbound(Linear_spring) to ubound(Linear_spring)
  with Linear_spring(i)
    .a = i
    .b = i+1
    .Lng.X = Mass(.b).Psn.X-Mass(.a).Psn.X
    .Lng.Y = Mass(.b).Psn.Y-Mass(.a).Psn.Y
    .Length = Sqr(.Lng.X*.Lng.X+.Lng.Y*.Lng.Y)
    .rest_length = .length
  end With
next

''	initiate angular springs
for i As Integer = lbound(angular_spring) to ubound(angular_spring)
  with angular_spring(i)
    .a = i
    .b = i+1
  end with
Next

''	initiate screen
ScreenRes Scrn_Wid, Scrn_hgt, 16
Color RGB(0, 0, 0), RGB(255, 255, 255)

Mouse_X = Mass(1).Psn.X
Mouse_Y = Mass(1).Psn.Y
SetMouse Mouse_x, Mouse_Y

Do
	
	''	use mouse to change spring angle
	GetMouse Mouse_x, Mouse_Y
	Mass(1).Psn.X = Mouse_X
	Mass(1).Psn.Y = Mouse_Y
	
	''	get linear spring data
	For i As Integer = Lbound(Linear_Spring) to Ubound(Linear_Spring)
		With Linear_Spring(i)
			''	spring length vector
			.Lng.X = Mass(.b).Psn.X-Mass(.a).Psn.X
			.Lng.Y = Mass(.b).Psn.Y-Mass(.a).Psn.Y
			''	spring length scalar
			.Length = Sqr(.Lng.X*.Lng.X+.Lng.Y*.Lng.Y)
			''	normalised spring length vector
			.Nor_Lng.X = .Lng.X/.Length
			.Nor_Lng.Y = .Lng.Y/.Length
			''	velocity (used for damping)
			.Vel.X = Mass(.b).Vel.X-Mass(.a).Vel.X
			.Vel.Y = Mass(.b).Vel.Y-Mass(.a).Vel.Y
		End With
	Next
	
	''	get angular spring data
	For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
		With Angular_Spring(i)
			''	sine of the angle = perpendicular dot product of normalized spring vectors
			.sin_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.y-Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.x
			''	cosine of the angle = dot product of normalized spring vectors
			.cos_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.x+Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.y
		End With
	Next
	
	''  apply linear spring force
	For i As Integer = Lbound(Linear_Spring) to Ubound(Linear_Spring)
		With Linear_Spring(i)
			Dim As Single Force = -Linear_Stiffnes*(.Length-.Rest_Length)-Linear_Damping*(.Nor_Lng.X*.Vel.X+.Nor_Lng.Y*.Vel.Y)
			Dim As Vector_Type Frc
			''	force vector = scalar force * normalised spring vector
			frc.x = Force*.Nor_Lng.X
			frc.y = Force*.Nor_Lng.Y
			Mass(.a).Frc.X -= Frc.X
			Mass(.a).Frc.Y -= Frc.Y
			Mass(.b).Frc.X += Frc.X
			Mass(.b).Frc.Y += Frc.Y
		End With
	Next
	
	''  apply angular spring force
	For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
		With Angular_Spring(i)
			Dim As Vector_Type Frc
			
      ''  angular spring force templates ranging from -180 to 180 degrees in 45 degree intervals
      ''	try these out and see how the "soft rod" stabilizes in different shapes!
      
      'Dim As Single Force = -Angular_Stiffnes*(-.Sin_Angle)							'' 	+/-180 deg
      'Dim As Single Force = -Angular_Stiffnes*(-.Sin_Angle-.Cos_Angle)		'' 		+135 deg
      'Dim As Single Force = -Angular_Stiffnes*(-.Cos_Angle)							''  	 +90 deg
      'Dim As Single Force = -Angular_Stiffnes*(.Sin_Angle-.Cos_Angle)		''  	 +45 deg
      Dim As Single Force = -Angular_Stiffnes*(.Sin_Angle)								''		   0 deg
      'Dim As Single Force = -Angular_Stiffnes*(.Sin_Angle+.Cos_Angle)		''  	 -45 deg
      'Dim As Single Force = -Angular_Stiffnes*(.Cos_Angle)								''  	 -90 deg
      'Dim As Single Force = -Angular_Stiffnes*(-.Sin_Angle+.Cos_Angle)		'' 		-135 deg
      
      
      ''	weighted average force keeps net torque zero and scales 
      ''	torque so it is always proportional to spring length.
      Dim As Single Force_a =  Force*Linear_spring(.b).Length
      Dim As Single Force_b = -Force*Linear_spring(.a).Length
      
			With Linear_spring(.a)
				''	force vector
				frc.x = Force_a* .Nor_Lng.y
				frc.y = Force_a*-.Nor_Lng.x
				''	apply force to point masses
				Mass(.a).Frc.x -= Frc.x
				Mass(.a).Frc.y -= Frc.y
				Mass(.b).Frc.x += Frc.x
				Mass(.b).Frc.y += Frc.y
			End with
			With Linear_spring(.b)
				''	force vector
				frc.x = Force_b* .Nor_Lng.y
				frc.y = Force_b*-.Nor_Lng.x
				''	apply force to point masses
				Mass(.a).Frc.x -= Frc.x
				Mass(.a).Frc.y -= Frc.y
				Mass(.b).Frc.x += Frc.x
				Mass(.b).Frc.y += Frc.y
			End with
		End With
	Next
	
	''	apply velocity squared air friction
	For i As Integer = lbound(mass) to ubound(mass)
		With mass(i)
			.Frc.x -= .Vel.X*Abs(.Vel.X)*Air_Friction
			.Frc.y -= .Vel.y*Abs(.Vel.y)*Air_Friction
		End With
	Next
	
	''	integrate (symplectic Euler 1st order)
	For i As Integer = lbound(mass)+2 to ubound(mass)
		With mass(i)
			.acc.x = (.frc.x/.mass)
			.acc.y = (.frc.y/.mass) + grav_acc
			.frc.x = 0
			.frc.y = 0
			.vel.x += .acc.x*dt
			.vel.y += .acc.y*dt
			.Psn.x += .vel.x*dt
			.Psn.y += .vel.y*dt
		End with
	Next
	
	''	draw to screen
	ScreenLock
	
		''	clear screen
		Cls 
		
		''	print data
		With Angular_Spring(1)
			Locate 2, 2: Print Using "Angular spring 1 angle: ####.## (-180, 180)"; ATan2(.Sin_Angle, .Cos_Angle)*rad2Deg
		End With
		
		''  draw springs
		For i As Integer = Lbound(Linear_Spring) To Ubound(Linear_Spring)
			With Linear_Spring(i)
				Line (Mass(.a).Psn.X, Mass(.a).Psn.Y)-(Mass(.b).Psn.X, Mass(.b).Psn.Y), Rgb(64, 160, 32)
			End With
		Next
		
		''  draw particles
		For i As Integer = Lbound(mass) To Ubound(mass)
			With mass(i)
				Circle(.Psn.X, .Psn.Y), .radius, rgb(0, 0, 0),,, 1, f
				Circle(.Psn.X, .Psn.Y), .radius-1, rgb(255, 32, 32),,, 1, f
			End With
		Next
		
	ScreenUnLock
	
	Sleep 1, 1
	
Loop Until Multikey(1)

End
But what if we want to assign a rest angle that isn't listed above? What we need is a generic all-purpose equation that will allow us to assign any rest angle and even change it during run-time. And here it is without further discussion:

Fa = -Ka * (cos (Vrest) * sin (V) – sin (Vr est)* cos (V))

Where Fa is angular force, Ka is angular spring stiffnes, V is current angle, and Vrest is rest angle. Notice how very similar this is to the perpdot equation? I don't know what it means, but it's cool! Here's a code sample showing the equation in action:

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's mass-spring-damper tutorial, November 2008
'
'   Example # 3.3
'
'		Angular spring simulation. This code sample is showing how linear 
'		and angular springs can form a string shaped soft body. 
'		Use mouse to move the soft body around and arrow keys to change
'		the size and shape of the body.
'
'  	(Press esc to quit)
'
'******************************************************************************'

''	define constants. Play around with these and see what happens!
const pi            		=	4*atn(1)		''	pi (better not change ^^)
Const rad2deg 					= 180/pi			''	radians to degrees
Const deg2rad						= Pi/180			''	degrees to radians
Const Grav_acc		  		= 2000				''	gravitational acceleration
const num_masses    		= 12					''	number of masses in rod (3 or higher)
Const Rod_Length   			= 300					''	total spring body length
Const Linear_Stiffnes   = 10000       ''  linear spring stiffnes
Const Linear_Damping    = 40         	''  linear spring damping
Const Angular_Stiffnes	= 10000      	''  angular spring stiffnes
Const Scrn_Wid      		= 600         ''  screen width
Const Scrn_Hgt      		= 600         ''  screen height
Const dt            		= 0.0025			''  timestep, delta time
Const	air_friction			= 0.002				''	air friction (to calm springs down)
Const Delta_Angle				= 0.1					''	spring angle step
Const Delta_Length			= 1.002				''	spring length step

''  define types
Type Vector_Type
  As Double X, Y
End Type

Type Mass_Type
  As Double Mass, radius
  As Vector_Type Frc, Acc, Vel, Psn
End Type 

Type Linear_Spring_Type
  As Integer a, b
  As Double Length, Rest_Length, Mass
  As Vector_Type Nor_Lng, Vel, Lng
End Type

Type Angular_Spring_Type
  As Integer a, b
  As Double sin_angle, cos_angle, sin_rest_angle, cos_rest_angle, Rest_Angle
End Type

''	dim variables
Dim As mass_type mass(1 to num_masses)
dim as Linear_spring_type Linear_spring(1 to num_masses-1)
dim as angular_spring_type angular_spring(1 to num_masses-2)
Dim As Integer mouse_x, mouse_y

Randomize timer
 
''	initiate point masses
for i As Integer = lbound(mass) to Ubound(mass)
  with mass(i)
    .mass = 1
    .radius = 4
    .Psn.x = (scrn_wid\2)
    .Psn.y = (scrn_hgt\4) + (Rod_Length/(num_masses-1))*i-1
  end with  
Next

''	initiate linear springs
for i As Integer = lbound(Linear_spring) to ubound(Linear_spring)
  with Linear_spring(i)
    .a = i
    .b = i+1
    .Lng.X = Mass(.b).Psn.X-Mass(.a).Psn.X
    .Lng.Y = Mass(.b).Psn.Y-Mass(.a).Psn.Y
    .Length = Sqr(.Lng.X*.Lng.X+.Lng.Y*.Lng.Y)
    .Nor_Lng.X = .Lng.X/.Length
    .Nor_Lng.Y = .Lng.Y/.Length
    .rest_length = .length
    .Mass = mass(.a).mass+mass(.b).mass
  end With
next

''	initiate angular springs
for i As Integer = lbound(angular_spring) to ubound(angular_spring)
  with angular_spring(i)
    .a = i
    .b = i+1
		.sin_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.y-Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.x
		.cos_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.x+Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.y
		.Sin_Rest_Angle = .Sin_Angle
		.Cos_Rest_Angle = .Cos_Angle
		.rest_angle = ATan2(.Sin_Rest_Angle, .Cos_Rest_Angle)
  end with
Next

''	initiate screen
ScreenRes Scrn_Wid, Scrn_hgt, 16
Color RGB(0, 0, 0), RGB(255, 255, 255)

Mouse_X = Mass(1).Psn.X
Mouse_Y = Mass(1).Psn.Y
SetMouse Mouse_x, Mouse_Y

Do
	
	''	use mouse to change spring angle
	GetMouse Mouse_x, Mouse_Y
	Mass(1).Psn.X = Mouse_X
	Mass(1).Psn.Y = Mouse_Y
	
	''  up arrow  ->  increase spring length
  If Multikey(&h48) Then
		For i As Integer = Lbound(Linear_Spring) To Ubound(Linear_Spring)
			With Linear_Spring(i)
        .Rest_Length *= Delta_Length
      End With
    Next
  End If
  
  ''  down arrow  -> decrease spring length
  If Multikey(&h50) Then
		For i As Integer = Lbound(Linear_Spring) To Ubound(Linear_Spring)
			With Linear_Spring(i)
        .Rest_Length /= Delta_Length
      End With
    Next
  End If
  
	''  Right arrow  ->  increase angular spring angle
  If Multikey(&h4d) Then
		For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
			With Angular_Spring(i)
        .Rest_Angle += Delta_Angle*deg2rad
        If .rest_angle < -Pi Then .rest_angle += 2*Pi
        If .rest_angle >  Pi Then .rest_angle -= 2*Pi
        .Sin_Rest_angle = Sin(.Rest_Angle)
        .Cos_Rest_angle = Cos(.Rest_Angle)
      End With
    Next
  End If
  
  ''  Left arrow  ->  decrease angular spring angle
  If Multikey(&h4b) Then
		For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
			With Angular_Spring(i)
        .Rest_Angle -= Delta_Angle*deg2rad
        If .rest_angle < -Pi Then .rest_angle += 2*Pi
        If .rest_angle >  Pi Then .rest_angle -= 2*Pi
        .Sin_Rest_angle = Sin(.Rest_Angle)
        .Cos_Rest_angle = Cos(.Rest_Angle)
      End With
    Next
  End If
	
	''	get linear spring data
	For i As Integer = Lbound(Linear_Spring) to Ubound(Linear_Spring)
		With Linear_Spring(i)
			''	spring length vector
			.Lng.X = Mass(.b).Psn.X-Mass(.a).Psn.X
			.Lng.Y = Mass(.b).Psn.Y-Mass(.a).Psn.Y
			''	spring length scalar
			.Length = Sqr(.Lng.X*.Lng.X+.Lng.Y*.Lng.Y)
			''	normalised spring length vector
			.Nor_Lng.X = .Lng.X/.Length
			.Nor_Lng.Y = .Lng.Y/.Length
			''	velocity (used for damping)
			.Vel.X = Mass(.b).Vel.X-Mass(.a).Vel.X
			.Vel.Y = Mass(.b).Vel.Y-Mass(.a).Vel.Y
		End With
	Next
	
	''	get angular spring data
	For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
		With Angular_Spring(i)
			''	sine of the angle = perpendicular dot product of normalized spring vectors
			.sin_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.y-Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.x
			''	cosine of the angle = dot product of normalized spring vectors
			.cos_angle = Linear_Spring(.a).Nor_Lng.x*Linear_Spring(.b).Nor_Lng.x+Linear_Spring(.a).Nor_Lng.y*Linear_Spring(.b).Nor_Lng.y
		End With
	Next
	
	''  apply linear spring force
	For i As Integer = Lbound(Linear_Spring) to Ubound(Linear_Spring)
		With Linear_Spring(i)
			Dim As Single Force = -Linear_Stiffnes*(.Length-.Rest_Length)-Linear_Damping*(.Nor_Lng.X*.Vel.X+.Nor_Lng.Y*.Vel.Y)
			Dim As Vector_Type Frc
			''	force vector = scalar force * normalised spring vector
			frc.x = Force*.Nor_Lng.X
			frc.y = Force*.Nor_Lng.Y
			Mass(.a).Frc.X -= Frc.X
			Mass(.a).Frc.Y -= Frc.Y
			Mass(.b).Frc.X += Frc.X
			Mass(.b).Frc.Y += Frc.Y
		End With
	Next
	
	''  apply angular spring force
	For i As Integer = Lbound(Angular_Spring) To Ubound(Angular_Spring)
		With Angular_Spring(i)
			Dim As Vector_Type Frc
			
			''	generic spring force equation
      Dim As Single Force = -Angular_Stiffnes*(.Cos_Rest_angle*.Sin_Angle-.Sin_Rest_angle*.Cos_angle)

      ''	weighted average force keeps net torque zero and scales 
      ''	torque so it is always proportional to spring length.
      Dim As Single Force_a =  Force*Linear_spring(.b).Length
      Dim As Single Force_b = -Force*Linear_spring(.a).Length
      
			With Linear_spring(.a)
				''	force vector = scalar force * normalised normal vector
				frc.x = Force_a*.Nor_Lng.y
				frc.y = Force_a*-.Nor_Lng.x
				''	apply force to point masses
				Mass(.a).Frc.x -= Frc.x
				Mass(.a).Frc.y -= Frc.y
				Mass(.b).Frc.x += Frc.x
				Mass(.b).Frc.y += Frc.y
			End with
			With Linear_spring(.b)
				''	force vector as above
				frc.x = Force_b*.Nor_Lng.y
				frc.y = Force_b*-.Nor_Lng.x
				''	apply force to point masses
				Mass(.a).Frc.x -= Frc.x
				Mass(.a).Frc.y -= Frc.y
				Mass(.b).Frc.x += Frc.x
				Mass(.b).Frc.y += Frc.y
			End with
		End With
	Next
	
	''	apply velocity squared air friction
	For i As Integer = lbound(mass) to ubound(mass)
		With mass(i)
			.Frc.x -= .Vel.X*Abs(.Vel.X)*Air_Friction
			.Frc.y -= .Vel.y*Abs(.Vel.y)*Air_Friction
		End With
	Next
	
	''	integrate (symplectic Euler 1st order)
	For i As Integer = lbound(mass)+1 to ubound(mass)
		With mass(i)
			.acc.x = (.frc.x/.mass)
			.acc.y = (.frc.y/.mass) + grav_acc
			.frc.x = 0
			.frc.y = 0
			.vel.x += .acc.x*dt
			.vel.y += .acc.y*dt
			.Psn.x += .vel.x*dt
			.Psn.y += .vel.y*dt
		End with
	Next
	
	''	draw to screen
	ScreenLock
	
		''	clear screen
		Cls 
		
		''	print data
		With Angular_Spring(1)
			Locate 2, 2: Print Using "     Angular spring 1 angle: ####.##"; ATan2(.Sin_Angle, .Cos_Angle)*rad2Deg
			Locate 4, 2: Print Using "Angular spring 1 rest angle: ####.##"; .rest_angle*rad2Deg
		End With
		
		''  draw springs
		For i As Integer = Lbound(Linear_Spring) To Ubound(Linear_Spring)
			With Linear_Spring(i)
				Line (Mass(.a).Psn.X, Mass(.a).Psn.Y)-(Mass(.b).Psn.X, Mass(.b).Psn.Y), Rgb(64, 160, 32)
			End With
		Next
		
		''  draw particles
		For i As Integer = Lbound(mass) To Ubound(mass)
			With mass(i)
				Circle(.Psn.X, .Psn.Y), .radius, rgb(0, 0, 0),,, 1, f
				Circle(.Psn.X, .Psn.Y), .radius-1, rgb(255, 32, 32),,, 1, f
			End With
		Next
		
	ScreenUnLock
	
	Sleep 1, 1
	
Loop Until Multikey(1)

End
That's all for this time, hope you enjoyed it.

Cheers,
Michael
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

This is really nice work. Thanks for this.

rb
veggie
Posts: 75
Joined: May 17, 2009 12:52

Great stuff!

Post by veggie »

Loved this tutorial series, pity it did not continue further... :(
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

Talking about a really weird coincidence :-)

Just today I finally cracked open the whole spring concept and figured out how to control it on a higher level. I've figured out how to implement unconditionally stable (damped or undamped) springs, and I'm able to explain exactly how it works. I've also expanded the spring concept into a really nice and simple type of constraint. I'll write a tutorial on this as soon as I've finished going through the maths and made a few proof-of-concept code samples.

Cheers,
Mike
veggie
Posts: 75
Joined: May 17, 2009 12:52

Post by veggie »

Good news Mike, its would be great to see some more of these excellent tutorials on what is usually a pretty dry subject, one which you manage to convey almost in the form of light reading... no mean feat!

Can't wait to see what you come up with. ;)
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

Post by h4tt3n »

veggie wrote:Good news Mike, its would be great to see some more of these excellent tutorials on what is usually a pretty dry subject, one which you manage to convey almost in the form of light reading... no mean feat!
Thanks man, that's exactly what I'm trying to do :-)
h4tt3n
Posts: 698
Joined: Oct 22, 2005 21:12
Location: Denmark

small teaser

Post by h4tt3n »

As promised, a new spring tutorial is coming up, this time on how to make the perfect spring, no less :-)

While I'm still struggling with some of the math, I decided to post a small teaser code sample. This spring example returns to it's rest length in ONE simulation step, which is as good as it gets. Besides, it is unconditionally stable, which means that it will not blow up, no matter what startup values you give it - timestep, mass, and startup velocity and position can have any value, and it will still behave perfectly.

The tutorial will revolve around how to derive the optimal spring stiffnes and damping values in order to achieve this. I don't think anyone has done this before :-)

Cheers,
Mike

Code: Select all

'******************************************************************************'
'
'   Michael "h4tt3n" Nissen's mass-spring-damper tutorial, December 2009
'
'   Example # 4.4
'
'		This damped spring example shows a "perfect" spring, which is unconditionally 
'		stable, and which returns to it's rest length in ONE simulation step, no matter
'		how distorted it is! Try setting the sleep_time constant to 1000, and you'll 
'		be able to see it happen.
'		Press left mouse to place particle. 
'		Press right mouse to change rest length.
'
'		(Press esc to quit)
'
'******************************************************************************'

''	define constants. Play around with these and see what happens!
const as integer 	screen_x								= 500					''	screen width
const as integer 	screen_y								= 500					''	screen height
const as integer 	center_x 								= screen_x\2	''	horizontal screen center
const as integer 	center_y	 							= screen_y\2	''	verticlal screen center
const as integer 	sleep_time	 						= 10					''	sleep time
const as single		pi											= 4*atn(1)		''	pi
const as single		particle_startup_angle	= 3/4 * pi		''	particle startup angle (0 - 2*pi)
const as single 	gravitational_acc		  	= 400					''	gravitational acceleration
const as integer 	spring_rest_length  		= 100					''	spring rest length
const as single 	dt            					= 0.01				''  timestep, delta time
const as single 	spring_stiffnes 				= 1						''  spring stiffnes (0 - 1)
const as single 	spring_damping 					= 1						''  spring damping  (0 - 1)

''  define types
type Vector_Type
  as single X, Y
end type

type Particle_Type
  as single Mass, radius
  as Vector_Type Frc, Acc, Vel, Psn
end type 

''	dim variables
dim as Particle_type Particle
dim as Vector_Type dst, nrm_dst
dim as integer mouse_x, mouse_y, mouse_b, frame
dim as single spring_length, spring_force, rest_distance, rest_length, normal_velocity

with Particle
  .mass 	= 1
  .radius = ((.Mass/0.001)/((4/3)*pi))^(1/3)
  .Psn.x 	= center_x + cos(0.5*pi + particle_startup_angle) * spring_rest_length
  .Psn.y 	= center_y + sin(0.5*pi + particle_startup_angle) * spring_rest_length
end with  

rest_length = spring_rest_length

''	initiate screen
screenres screen_x, screen_y, 16
color rgb(0, 0, 0), rgb(255, 255, 255)

do
	
	frame += 1
	
	getmouse mouse_x, mouse_y,, mouse_b
	
	''	on right klick set new rest length
	if mouse_b = 2 then
		dst.X = mouse_X - center_x:	dst.Y = mouse_Y - center_y
		rest_length = sqr(dst.X*dst.X+dst.Y*dst.Y)
	end if
	
	''	on leftclick set particle to mouse position
	if mouse_b = 1 then
		with Particle
			.Frc.X = 0			:	.Frc.y = 0
			.vel.X = 0			:	.vel.y = 0
			.psn.x = mouse_x:	.psn.y = mouse_y
		end with
	end if
	
	''	gravitational acceleration
	Particle.vel.y += gravitational_acc*dt
	
	''	distance vector
	dst.X = Particle.Psn.X - center_x
	dst.Y = Particle.Psn.Y - center_y
	
	''	distance scalar
	spring_length = sqr(dst.X*dst.X+dst.Y*dst.Y)
	
	''	normalise distance vector
	if spring_length > 0 then
		nrm_dst.x = dst.x/spring_length
		nrm_dst.y = dst.y/spring_length
	else
		nrm_dst.x = 0
		nrm_dst.y = 0
	end if
	
	''	particle velocity along spring axis
	normal_velocity = nrm_dst.x*particle.vel.x+nrm_dst.y*particle.vel.y
	
	''	difference between spring length and rest length
	rest_distance = spring_length - rest_length
	
	''	Damped spring equation f = -k * x - c * v (expressed in different ways below)
	''	(spring stiffnes coefficient k is automatically set to m / dt^2)
	''	(spring damping coefficient c is automatically set to m / dt)
	'spring_force = -spring_stiffnes*Particle.Mass/(dt*dt)*rest_distance - spring_damping*Particle.mass/dt*normal_velocity
	spring_force = -(rest_distance*spring_stiffnes/(dt*dt) + normal_velocity*spring_damping/dt)*Particle.mass
	
	''	apply force vector
	Particle.Frc.X = spring_force*nrm_Dst.x
	Particle.Frc.y = spring_force*nrm_Dst.y
	
	''	draw to screen
	screenlock
	
		''	clear screen
		cls
		
		''	print data
		locate 2, 2:	print using "rest length: ###.####"; rest_length
		locate 4, 2:	print using "real length: ###.####"; spring_length
		
		locate 2, screen_x\8 - 12: print using "frame: ######"; frame
		
		''	draw fixed point in center of screen
		circle(center_x, center_y), 4, rgb(128, 128, 128),,, 1, f
		
				''  draw free moving particle
		circle(Particle.Psn.X, Particle.Psn.Y), Particle.radius, rgb(255, 0, 0),,, 1, f
		
		''	draw circle showing actual particle distance
		circle(center_x, center_y), spring_length, rgb(192, 192, 192),,, 1
		
		''	draw circle showing all valid particle positions for the given rest length
		circle(center_x, center_y), rest_length, rgb(192, 192, 192),,, 1
		
		''	draw spring
		line(center_x, center_y)-(Particle.Psn.X, Particle.Psn.Y), rgb(0, 192, 0)
		
	screenunlock
	
	''	integrate (symplectic Euler 1st order)
	with Particle
		.acc.x  = .frc.x/.mass:	.acc.y  = .frc.y/.mass
		.vel.x += .acc.x*dt		:	.vel.y += .acc.y*dt
		.Psn.x += .vel.x*dt		:	.Psn.y += .vel.y*dt
	end with
	
	sleep sleep_time, 1
	
loop until multikey(1)

end
veggie
Posts: 75
Joined: May 17, 2009 12:52

Post by veggie »

Bah, where's my "Notify me when a reply is posted"? Damn spam filter!

Anyway', lookin real good Mike, the subject matter, "derive the optimal spring stiffnes and damping values" is a very important spring physics aspect to me personally as its one of the features I was (not very successfuly :P) trying to figure out.

Will keep a closer eye on this thread for updates (via an old fashioned link! :))
Post Reply