Collision impulses

.gitignore

@@ -1 +1,2 @@
 build/
+cover.png

collision.lua

@@ -37,9 +37,12 @@ local function DetectCollision( curMarble, newMarble, curWave, newWave )
   local iWaveRadius = curWave:Interpolate( iTheta )
   local fWaveRadius = curWave:Interpolate( fTheta )
   
+  local speed = newMarble.dx * newMarble.dx + newMarble.dy * newMarble.dy
+
   --Case 1: marble is already outside the curve.
   if iRadius > iWaveRadius then
     return {
+      speed = speed,
       dt = newMarble.t - curMarble.t,
       t = curMarble.t,
       r = iWaveRadius,
@@ -51,6 +54,7 @@ local function DetectCollision( curMarble, newMarble, curWave, newWave )
   if fRadius > fWaveRadius then
     --This isn't right, fix later.
     return {
+      speed = speed,
       dt = 0,
       t = newMarble.t,
       r = fWaveRadius,

main.lua

@@ -1,5 +1,5 @@
 local love = love
-local step = 1.0 / 120.0
+local step = 1.0 / 120
 local sitelenpona
 local text
 local marble
@@ -107,7 +107,8 @@ local function BeatScore( t )
 
   --Safety: avoid a dbz.
   if mu < 0.001 or dt < 0.001 then
-    error( "DBZ! Beat length too small." )
+    return 0.0
+    --error( "DBZ! Beat length too small." )
   end
 
   --Calculate beat score.
@@ -171,6 +172,31 @@ function love.draw()
   love.graphics.push( "transform" )
   love.graphics.applyTransform( transform )
   wave.Draw()
+  
+  if debugRenderImpact then
+  
+  love.graphics.setLineWidth( 0.01 )
+  love.graphics.setColor( 1, 0, 0, 0.5 ) --Red: Incoming
+  love.graphics.line(
+    debugRenderImpact.xi,
+    debugRenderImpact.yi,
+    debugRenderImpact.xf,
+    debugRenderImpact.yf)
+  love.graphics.setColor( 0, 1, 0, 0.5 ) --Green: Normal
+  love.graphics.line(
+    debugRenderImpact.xi,
+    debugRenderImpact.yi,
+    debugRenderImpact.xn,
+    debugRenderImpact.yn)
+  love.graphics.setColor( 0, 0, 1, 0.5 ) -- Blue: Outgoing
+    love.graphics.line(
+    debugRenderImpact.xi,
+    debugRenderImpact.yi,
+    debugRenderImpact.vxout,
+    debugRenderImpact.vyout)
+  
+  end
+  
   love.graphics.pop()
 
   sitelenpona.Draw( text.tok[state.currentBeat] )

marble.lua

@@ -26,7 +26,47 @@ end
 local function OnImpact( impact )
   --Adjust current trajectory according to collision.
   if not impact.dt then return end
-  return Integrate( impact.dt ) --Hmm! Maybe this should be a fixed step instead for stability's sake.
+  
+  local x, y = impact.r * math.cos( impact.th ), impact.r * math.sin( impact.th )
+  local vx, vy = newState.dx, newState.dy --Velocity of particle going into collision.
+  local unx, uny = math.cos( impact.normal ), math.sin( impact.normal ) --Outward-facing normal of wave.
+  local uvx, uvy --Unit vector velocity of particle.
+  local speed = math.sqrt( vx * vx + vy * vy )
+  if speed < 0.01 then
+    --If the marble is motionless, there is no angular velocity wrt 0,
+    --so the wave is headed directly inward.
+    --We handle the collision as if the marble is headed directly outward.
+    uvx, uvy = unx, uny
+    vx, vy = uvx, uvy
+  else
+    uvx, uvy = vx / speed , vy / speed
+  end
+  
+  
+  
+  --Get signed angle between normal and incoming velocity (both unit vectors)
+  local dot = unx * uvy - uny * uvx
+  
+  --Fudge factor: apply an impulse inward so that you don't stick or slide on the wave.
+  local inward = ( dot > 0 ) and dot or -dot
+  inward = inward * inward * inward
+  
+  --Calculate the rotation matrix:
+  --counterclockwise rotation by  2 * pi - 2 * arccos( n dot v )
+  local c, s = 1 - 2 * dot * dot, - 2 * dot * math.sqrt( 1 - dot * dot )
+  --Apply:
+  local vxout, vyout =
+  inward * (- math.cos(impact.th) ) - (1.0 - inward) * ( vx * c - vy * s ),
+  inward * (- math.sin(impact.th) ) - (1.0 - inward) * ( x * s + vy * c )
+  
+  curState.x, curState.y = x, y
+  curState.dx, curState.dy = 0.5 * vxout, 0.5 * vyout
+  
+  debugRenderImpact = { xi = x, yi = y, xf = x - vx, yf = y - vy,
+    xn = 0.2 * unx + x, yn = 0.2 * uny + y,
+    vxout = x + vxout, vyout = y + vyout}
+  
+  return Integrate( math.max( impact.dt , 1 / 60 ) ) --Hmm! Maybe this should be a fixed step instead for stability's sake.
 end
 
 local function Update()

wave.lua

@@ -1,15 +1,23 @@
 --Render and simulate 1D wave equation.
 local love = love
-local old, cur, new --States add beginning of last tick, current tick, and next tick respectively.
+local N = 15
+local SOUNDSPEED = 0.5
+local IMPULSESIZE = 20
+local DAMPING = 0.01
 
-local N = 17
+
-local SOUNDSPEED = 2.0
+local old, cur, new --States add beginning of last tick, current tick, and next tick respectively.
 
 local function Current() return cur end
 local function Next() return new end
+local Integrate
+local Interpolate
+local Derivative
+local SecondDerivative
+local DFT
+
 
 --Calculate discrete fourier transform of radius function.
-local DFT
 do
   local twiddlec = {}
   local twiddles = {}
@@ -43,7 +51,7 @@ do
 end
 
 --Minimal-oscillation interpolant of a real function from its discrete Fourier coefficients.
-local Interpolate = function( wave, x )
+Interpolate = function( wave, x )
   local re, im = wave.dftre, wave.dftim
   local y = re[1]
 
@@ -61,7 +69,7 @@ local Interpolate = function( wave, x )
 end
 
 --Derivative of the interpolation.
-local Derivative = function( wave, x )
+Derivative = function( wave, x )
   local re, im = wave.dftre, wave.dftim
   local y = 0
 
@@ -79,7 +87,7 @@ local Derivative = function( wave, x )
 end
 
 --Second derivative of the interpolation.
-local SecondDerivative = function( wave, x )
+SecondDerivative = function( wave, x )
   local re, im = wave.dftre, wave.dftim
   local y = 0
 
@@ -103,13 +111,20 @@ local function AliasedSinc( theta, x )
 end
 
 --Apply bandlimited impulse to wave.
-local Impulse = function( wave, location, magnitude )
+local function OnImpact( impact )
-  local speed = wave.vrad
   
+  local r = cur.radii
+  local theta = impact.th
+  local magnitude = IMPULSESIZE * impact.speed
+  local dt = math.max( impact.dt, 1 / 120.0 )
   for i = 0, N - 1 do
-    speed[ i + 1 ] = speed[ i + 1 ] + magnitude * AliasedSinc( 2.0 * math.pi * i, location )
+    r[ i + 1 ] = r[ i + 1 ] + dt * magnitude * AliasedSinc(  theta, 2.0 * math.pi * i / N )
   end 
   
+  --We've updated the positions, now we need to take a DFT
+  --in order to get the bandlimited second spatial derivative.
+  cur:DFT() 
+  return Integrate( dt )
 end
 
 local mt = { __index = {
@@ -134,7 +149,7 @@ local function Wave( )
   }
 
   for i = 1, N do
-    t.radii[i] = 1.0 + 0.05 * math.sin( i * 2.0 * math.pi / N  )
+    t.radii[i] = 1.0 + 0.05 * ( i/N + math.sin( i * 2.0 * math.pi / N  ))
     t.vrad[i] = 0.0
   end
   DFT( t )
@@ -166,11 +181,11 @@ local function Draw()
       th = ( i - 1 + k ) * 2.0 * math.pi / N
       local r = cur:Interpolate( th )
       local x, y = r * math.cos( th ), r * math.sin( th )
-      --love.graphics.circle( "fill", x, y, 0.01 )
+      love.graphics.circle( "fill", x, y, 0.01 )
-      love.graphics.circle( "fill", th / math.pi - 1.0 , r, 0.01)
+      --love.graphics.circle( "fill", th / math.pi - 1.0 , r, 0.01)
       
       --First derivative.
-      love.graphics.setColor( 0, 1.0, 0, 0.7 )
+      --[[love.graphics.setColor( 0, 1.0, 0, 0.7 )
       r = cur:Derivative( th )
       x, y = r * math.cos( th ), r * math.sin( th )
       --love.graphics.circle( "fill", x, y, 0.01 )
@@ -184,7 +199,7 @@ local function Draw()
       love.graphics.circle( "fill", th / math.pi - 1.0, r, 0.01)
       
       love.graphics.setColor( 1.0, 1.0, 1.0, 0.2 )
-      love.graphics.circle( "fill", 2.0 * ( i + k ) / N - 1.2, 0, 0.02 )
+      love.graphics.circle( "fill", 2.0 * ( i + k ) / N - 1.2, 0, 0.02 )]]
       
     end
   end
@@ -195,9 +210,6 @@ local function Draw()
   love.graphics.circle( "fill", x, y, 0.02 )
 end
 
-local function OnImpact( impact )
-  
-end
 
 local function Update()
 
@@ -216,12 +228,15 @@ local function Update()
   end
 end
 
-local function Integrate( step )
+Integrate = function( step )
   
   for i = 1, N do
-    --new.vrad[i] = 
     local rxx = cur:SecondDerivative( math.pi * 2.0 * ( i - 1 ) / N )
-    new.radii[i] = 2.0 * cur.radii[i] - old.radii[i] + step * step * SOUNDSPEED * rxx
+    local r = ( 1.0 - DAMPING ) * ( 2.0 * cur.radii[i] - old.radii[i] + step * step * SOUNDSPEED * rxx ) --Verlet
+     + DAMPING --Damping: oscillate toward 1.
+    if r > 1.5 then r = 1.5 end
+    if r < 0.5 then r = 0.5 end
+    new.radii[i] = r
   end
   
   new:DFT( )