## h4tt3n's spring physics simulation tutorial

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

### h4tt3n's spring physics simulation tutorial

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: 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: 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 reactsconst Pi = 4*Atn(1)             ''  Pi (better not change this ^^)const Mass_1_mass = 100         ''  set mass of the spring's endpointsconst mass_2_mass = 20          ''  set mass of the spring's endpointsconst spring_distortion = 50    ''  how stretched the spring is at startupconst spring_rest_length = 200  ''  spring neutral length or rest lengthconst spring_stiffnes = 50      ''  spring stiffnes. Larger equals stifferconst timestep = 0.01           ''  delta time. Smaller is slower but more accurate''  dimension variablesDim 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 massesMass_1_Position = 300-(spring_rest_length+spring_distortion)/2Mass_2_Position = Mass_1_Position + spring_rest_length + spring_distortion''  set mass radius using the equation of a sphere's volumeMass_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 pixelScreenres 600, 400, 16''  main program loopDo    ''  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 valuesEnd`

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:
Great... will this be in QBE as well?
h4tt3n
Posts: 681
Joined: Oct 22, 2005 21:12
Location: Denmark
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:
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:
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: 681
Joined: Oct 22, 2005 21:12
Location: Denmark
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: 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 reactsConst Pi                =     4*Atn(1)          ''  pi (better not change ^^)Const Timestep          =     0.01              ''  timestepConst Num_Masses        =     5              -1 ''  number of masses in ropeConst Point_Mass        =     10                ''  mass of each point massConst Point_Density     =     0.01              ''  density of each point massConst Rope_Length       =     350               ''  rope lengthConst Spring_Stiffnes   =     6000              ''  spring stiffnesConst Grav_Acc          =     800               ''  gravitational accelerationConst Rest_Length = Rope_Length/Num_Masses      ''  rest length of each spring''  define types''  two dimensional vector typeType Vector_Type  As Single X, YEnd Type''  mass type including force-, acceleration-, velocity-, and position vectorsType Mass_Type  As Single Mass, Density, Radius  As Vector_Type Frc, Acc, Vel, PosEnd Type ''  dimension variablesDim As Vector_Type Length, Normalised_Length, VelDim As Mass_Type Mass(Num_Masses)Dim As Single Force, Spring_LengthDim As Integer i, i2, Scrn_Wid, Scrn_Hgt''  set screen width, height, and bit depthScrn_Wid = 800Scrn_Hgt = 600ScreenRes Scrn_Wid, Scrn_hgt, 16''  set startup conditions of massesFor 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 loopDo    ''  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: 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 ballsConst Mass_Min          =     400           ''  smallest elastic ball massConst Mass_Max          =     8000         ''  biggest elastic ball massConst Density           =     0.01          ''  elastic ball density Const Spring_Stiffnes   =     5e5           ''  elastic stiffnes''  define typesType Vector_Type  As Single X, YEnd TypeType Ball_Type  As Uinteger Col  As Single Mass, Density, Radius  As Vector_Type Frc, Acc, Vel, PosEnd Type ''  dimension variablesDim As Vector_Type Dst, VelDim As Ball_Type Ball(Num_Balls)Dim As Single Force, Ball_Distance, Cos_Angle, Sin_Angle, Dist_MinDim As Integer i, i2, Scrn_Wid, Scrn_Hgt''  set screen width, height, and bit depthScrn_Wid = 800Scrn_Hgt = 600ScreenRes Scrn_Wid, Scrn_hgt, 16''  use timer to generate random numbersRandomize Timer''  set startup conditions of elastic ballsFor 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 loopDo    ''  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: 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: 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 reactsConst Pi                =     4*Atn(1)          ''  pi (better not change ^^)Const Timestep          =     0.01              ''  timestepConst Num_Masses        =     5             -1 ''  number of masses in ropeConst Point_Mass        =     10                 ''  mass of each point massConst Point_Density     =     0.01              ''  density of each point massConst Rope_Length       =     350               ''  rope lengthConst Spring_Stiffnes   =     6000             ''  spring stiffnesConst Spring_Damping    =     2e2              ''  spring dampingConst Grav_Acc          =     800              ''  gravitational accelerationConst Rest_Length = Rope_Length/Num_Masses      ''  rest length of each spring''  define typesType Vector_Type  As Single X, YEnd TypeType Mass_Type  As Single Mass, Density, Radius  As Vector_Type Frc, Acc, Vel, PosEnd Type ''  dimension variablesDim As Vector_Type Lng, Norm_Lng, VelDim As Mass_Type Mass(Num_Masses)Dim As Single Force, Spring_LengthDim As Integer i, i2, Scrn_Wid, Scrn_Hgt''  set screen width, height, and bit depthScrn_Wid = 800Scrn_Hgt = 600ScreenRes Scrn_Wid, Scrn_hgt, 16''  set startup conditions of massesFor 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 WithNext''  main program loopDo    ''  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:
WOW, GREAT TUTORIAL!!!!
AWESOME AS PREVIOUS ONE!!!

Keep up the good work!
h4tt3n
Posts: 681
Joined: Oct 22, 2005 21:12
Location: Denmark
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. 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 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 constantsconst pi            = 4*atn(1)      ''   piConst rad2deg          = 180/pi         ''   radians to degreesConst num_masses       = 3                  ''   number of point massesConst Rope_Length      =   350               ''   total length of all springsConst Scrn_Wid      = 600               ''  screen widthConst Scrn_Hgt      = 600               ''  screen height''  define typesType Vector_Type  As Double X, YEnd TypeType Mass_Type  As Double Mass, dens, radius  As Vector_Type Frc, Acc, Vel, PsnEnd Type Type Spring_Type  As Integer a, b  As Double Length  As Vector_Type Nor_Lng, LngEnd TypeType Angular_Spring_Type  As Integer a, b  As Double sin_angle, cos_angleEnd Type''   dim arrays & variablesDim 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 massesfor 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 withnext''   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 withNext''   initiate screenScreenRes Scrn_Wid, Scrn_hgt, 16Color 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. 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. 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 degreesConst deg2rad                  = Pi/180         ''   degrees to radiansConst Grav_acc              = 2000            ''   gravitational accelerationconst num_masses          = 8                  ''   number of masses in rodConst Rod_Length            = 300               ''   total spring body lengthConst Linear_Stiffnes     = 10000       ''  linear spring stiffnesConst Linear_Damping    = 40            ''  linear spring dampingConst Angular_Stiffnes   = 10000         ''  angular spring stiffnesConst Scrn_Wid            = 600         ''  screen widthConst Scrn_Hgt            = 600         ''  screen heightConst dt                  = 0.0025         ''  timestep, delta timeConst   air_friction         = 0.002            ''   air friction (to calm springs down)''  define typesType Vector_Type  As Double X, YEnd TypeType Mass_Type  As Double Mass, radius  As Vector_Type Frc, Acc, Vel, PsnEnd Type Type Linear_Spring_Type  As Integer a, b  As Double Length, Rest_Length  As Vector_Type Nor_Lng, Vel, LngEnd TypeType Angular_Spring_Type  As Integer a, b  As Double sin_angle, cos_angleEnd Type''   dim variablesDim 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_yRandomize timer ''   initiate point massesfor 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 springsfor 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 Withnext''   initiate angular springsfor i As Integer = lbound(angular_spring) to ubound(angular_spring)  with angular_spring(i)    .a = i    .b = i+1  end withNext''   initiate screenScreenRes Scrn_Wid, Scrn_hgt, 16Color RGB(0, 0, 0), RGB(255, 255, 255)Mouse_X = Mass(1).Psn.XMouse_Y = Mass(1).Psn.YSetMouse Mouse_x, Mouse_YDo      ''   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 degreesConst deg2rad                  = Pi/180         ''   degrees to radiansConst Grav_acc              = 2000            ''   gravitational accelerationconst num_masses          = 12               ''   number of masses in rod (3 or higher)Const Rod_Length            = 300               ''   total spring body lengthConst Linear_Stiffnes   = 10000       ''  linear spring stiffnesConst Linear_Damping    = 40            ''  linear spring dampingConst Angular_Stiffnes   = 10000         ''  angular spring stiffnesConst Scrn_Wid            = 600         ''  screen widthConst Scrn_Hgt            = 600         ''  screen heightConst dt                  = 0.0025         ''  timestep, delta timeConst   air_friction         = 0.002            ''   air friction (to calm springs down)Const Delta_Angle            = 0.1               ''   spring angle stepConst Delta_Length         = 1.002            ''   spring length step''  define typesType Vector_Type  As Double X, YEnd TypeType Mass_Type  As Double Mass, radius  As Vector_Type Frc, Acc, Vel, PsnEnd Type Type Linear_Spring_Type  As Integer a, b  As Double Length, Rest_Length, Mass  As Vector_Type Nor_Lng, Vel, LngEnd TypeType Angular_Spring_Type  As Integer a, b  As Double sin_angle, cos_angle, sin_rest_angle, cos_rest_angle, Rest_AngleEnd Type''   dim variablesDim 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_yRandomize timer ''   initiate point massesfor 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 springsfor 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 Withnext''   initiate angular springsfor 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 withNext''   initiate screenScreenRes Scrn_Wid, Scrn_hgt, 16Color RGB(0, 0, 0), RGB(255, 255, 255)Mouse_X = Mass(1).Psn.XMouse_Y = Mass(1).Psn.YSetMouse Mouse_x, Mouse_YDo      ''   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
This is really nice work. Thanks for this.

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

### Great stuff!

Loved this tutorial series, pity it did not continue further... :(
h4tt3n
Posts: 681
Joined: Oct 22, 2005 21:12
Location: Denmark
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
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: 681
Joined: Oct 22, 2005 21:12
Location: Denmark
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: 681
Joined: Oct 22, 2005 21:12
Location: Denmark

### small teaser

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 widthconst as integer    screen_y                        = 500               ''   screen heightconst as integer    center_x                         = screen_x\2   ''   horizontal screen centerconst as integer    center_y                         = screen_y\2   ''   verticlal screen centerconst as integer    sleep_time                      = 10               ''   sleep timeconst as single      pi                                 = 4*atn(1)      ''   piconst as single      particle_startup_angle   = 3/4 * pi      ''   particle startup angle (0 - 2*pi)const as single    gravitational_acc           = 400               ''   gravitational accelerationconst as integer    spring_rest_length        = 100               ''   spring rest lengthconst as single    dt                           = 0.01            ''  timestep, delta timeconst as single    spring_stiffnes             = 1                  ''  spring stiffnes (0 - 1)const as single    spring_damping                = 1                  ''  spring damping  (0 - 1)''  define typestype Vector_Type  as single X, Yend typetype Particle_Type  as single Mass, radius  as Vector_Type Frc, Acc, Vel, Psnend type ''   dim variablesdim as Particle_type Particledim as Vector_Type dst, nrm_dstdim as integer mouse_x, mouse_y, mouse_b, framedim as single spring_length, spring_force, rest_distance, rest_length, normal_velocitywith 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_lengthend with  rest_length = spring_rest_length''   initiate screenscreenres screen_x, screen_y, 16color 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
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! :))