STABLE! STABLE wave equation.

collision.lua

@@ -1,6 +1,67 @@
+local function Polar( x, y )
+  local arctan = math.atan2( y, x )
+  --maybe we don't need this line?
+  --if arctan < 0.0 then arctan = arctan + 2.0 * math.pi end
+  return math.sqrt( x * x + y * y ), arctan
+end
+
+local function tPolar( t )
+  return Polar( t.x, t.y )
+end
+
+local function LerpWave( cur, new, alpha, x )
+  local a = cur:Interpolate( x )
+  local b = new:Interpolate( x )
+  return alpha * b + ( 1.0 - alpha ) * a
+end
+
+--Peak of alpha * cur + (1 - alpha) * new on interval [xi, xf]
+local function GetMinimum( cur, new, alpha, xi, xf )
+  local ri, rf = LerpWave( cur, new, alpha, xi ), LerpWave( cur, new, alpha, xf )
+  local xm = 0.5 * ( xi + xf )
+  local rm = LerpWave( cur, new, alpha, 0.5 * ( xi + xf ) )
+  return math.min( ri, rf )
+end
+
+local function GetTimeFromAlpha( alpha, curM, newM )
+  return alpha * newM.t + ( 1.0 - alpha ) * curM.t
+end
+
 local function DetectCollision( curMarble, newMarble, curWave, newWave )
-  local t, r, th, normal
-  if t then return { t = t, r = r, th = th, normal = normal } end
+  local t, dt, r, th, normal
+
+
+  --Polar coordinates:
+  local iRadius, iTheta = tPolar( curMarble )
+  local fRadius, fTheta = tPolar( newMarble )
+  local iWaveRadius = curWave:Interpolate( iTheta )
+  local fWaveRadius = curWave:Interpolate( fTheta )
+
+  --Case 1: marble is already outside the curve.
+  if iRadius > iWaveRadius then
+    return {
+      dt = newMarble.t - curMarble.t,
+      t = curMarble.t,
+      r = iWaveRadius,
+      th = iTheta,
+      normal = iTheta - math.atan( curWave:Derivative( iTheta ) ), }
+  end
+  
+  --Case 2: marble is outside the curve by tick end.
+  if fRadius > fWaveRadius then
+    --This isn't right, fix later.
+    return {
+      dt = 0,
+      t = newMarble.t,
+      r = fWaveRadius,
+      th = fTheta,
+      normal = fTheta - math.atan( curWave:Derivative( iTheta ) ) }
+  end
+  
+  --Case 3: check whether marble passed through analog peak.
+  
+  --Case 4: marble must have remained inside curve.
+  return
 end
 
 return DetectCollision

marble.lua

@@ -1,6 +1,5 @@
 --Character controller. Renders point particle.
 local love = love
---local wave = assert( require "wave" )
 local oldBuffer = love.graphics.newCanvas()
 local newBuffer = love.graphics.newCanvas()
 local oldState, curState, newState
@@ -26,8 +25,8 @@ end
 
 local function OnImpact( impact )
   --Adjust current trajectory according to collision.
-  if not impact.tNext then return end
-  return Integrate( impact.tNext - impact.t ) --Hmm! Maybe this should be a fixed step instead for stability's sake.
+  if not impact.dt then return end
+  return Integrate( impact.dt ) --Hmm! Maybe this should be a fixed step instead for stability's sake.
 end
 
 local function Update()
@@ -79,9 +78,9 @@ local function Draw()
   --Render circle directly to screen.
   love.graphics.setCanvas()
   love.graphics.draw( newBuffer )
-  --love.graphics.circle( "fill", xf, yf, 5.0 )
   love.graphics.setColor( 1, 1, 1, 1.0 ) --White.
-  --love.graphics.circle( "line", xf, yf, 5.0 )
+  love.graphics.setLineWidth( 1.0 )
+  love.graphics.circle( "line", xn, yn, 4 )
   
   oldBuffer, newBuffer = newBuffer, oldBuffer
 

wave.lua

@@ -2,7 +2,8 @@
 local love = love
 local old, cur, new --States add beginning of last tick, current tick, and next tick respectively.
 
-local N = 13
+local N = 17
+local SOUNDSPEED = 2.0
 
 local function Current() return cur end
 local function Next() return new end
@@ -133,7 +134,7 @@ local function Wave( )
   }
 
   for i = 1, N do
-    t.radii[i] = 0.3 * math.sin( i * 2.0 * math.pi / N  )
+    t.radii[i] = 1.0 + 0.05 * math.sin( i * 2.0 * math.pi / N  )
     t.vrad[i] = 0.0
   end
   DFT( t )
@@ -216,12 +217,14 @@ local function Update()
 end
 
 local function Integrate( step )
-  local t = love.timer.getTime()
+  
   for i = 1, N do
-    new.radii[i] = 0.05 * math.cos( math.sin( t * 0.05 + i ) )
+    --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
   end
   
-  cur:DFT( )
+  new:DFT( )
 
 end