'******************************************************************************'
'
'   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