your-own-drum/wave.lua

249 lines
5.7 KiB
Lua

--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 = 13
local function Current() return cur end
local function Next() return new end
--Calculate discrete fourier transform of radius function.
local DFT
do
local twiddlec = {}
local twiddles = {}
for i = 1, N do
twiddlec[i], twiddles[i] = math.cos( - 2.0 * ( i - 1 ) * math.pi / N ), math.sin( - 2.0 * (i - 1) * math.pi / N )
end
local function Twiddle( j )
j = 1 + j % N
return twiddlec[j], twiddles[j]
end
--Slow Discrete Fourier Transform of a real sequence.
DFT = function( wave )
local seq = wave.radii
local dre, dim = wave.dftre, wave.dftim
for k = 0, N - 1 do
local x, y = 0, 0 --Fourier coefficients in bin k.
for n = 0, N - 1 do
local cos, sin = Twiddle( n * k )
x = x + seq[n + 1] * cos
y = y + seq[n + 1] * sin
end
dre[k + 1], dim[k + 1] = x / N , y / N
end
end
end
--Minimal-oscillation interpolant of a real function from its discrete Fourier coefficients.
local Interpolate = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = re[1]
for k = 1, N / 2 do
local c, s = math.cos( x * k ), math.sin( x * k )
y = y
+ c * re[k + 1]
- s * im[k + 1]
+ c * re[N - k + 1]
+ s * im[N - k + 1]
end
return y
end
--Derivative of the interpolation.
local Derivative = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = 0
for k = 1, N / 2 do
local c, s = k * math.cos( x * k ), k * math.sin( x * k )
y = y
- c * im[k + 1]
- s * re[k + 1]
+ c * im[N - k + 1]
- s * re[N - k + 1]
end
return y
end
--Second derivative of the interpolation.
local SecondDerivative = function( wave, x )
local re, im = wave.dftre, wave.dftim
local y = 0
for k = 1, N / 2 do
local c, s = -k * k * math.cos( x * k ), -k * k * math.sin( x * k )
y = y
+ c * re[k + 1]
- s * im[k + 1]
+ c * re[N - k + 1]
+ s * im[N - k + 1]
end
return y
end
--Bandlimited impulse located at angle theta.
local function AliasedSinc( theta, x )
local n, d = math.sin( N * 0.5 * ( x - theta ) ), N * math.sin( 0.5 * ( x - theta ) )
if d < 0.0001 and d > -0.0001 then return 1.0 end
return n / d
end
--Apply bandlimited impulse to wave.
local Impulse = function( wave, location, magnitude )
local speed = wave.vrad
for i = 0, N - 1 do
speed[ i + 1 ] = speed[ i + 1 ] + magnitude * AliasedSinc( 2.0 * math.pi * i, location )
end
end
local mt = { __index = {
DFT = DFT,
Interpolate = Interpolate,
Derivative = Derivative,
SecondDerivative = SecondDerivative,
Impulse = Impulse } }
local function Wave( )
local t = {
--radii[k] = radius of point on curve at angle (k - 1) / N
radii = {},
--TIME derivative of radius
vrad = {},
--SPACE DFT of radius function (which is periodic)
dftre = {},
dftim = {},
}
for i = 1, N do
t.radii[i] = 0.3 * math.sin( i * 2.0 * math.pi / N )
t.vrad[i] = 0.0
end
DFT( t )
return setmetatable(t, mt)
end
local function Draw()
-- Blue circle.
love.graphics.setColor( 91 / 255, 206 / 255, 250 / 255 )
love.graphics.circle("fill", 0, 0, 1)
local t = love.timer.getTime()
-- Debug dots.
for i = 1, N do
local th = ( i - 1 ) * 2.0 * math.pi / N
local cx, cy = cur.radii[i] * math.cos( th ), cur.radii[i] * math.sin( th )
love.graphics.setCanvas()
love.graphics.setColor( 0, 0, 0, 0.5 )
love.graphics.circle( "fill", cx, cy, 0.02 )
for k = 0.1, 1.0, 0.1 do
--Interpolant.
love.graphics.setColor( 1.0, 0, 0, 0.7 )
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", th / math.pi - 1.0 , r, 0.01)
--First derivative.
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 )
love.graphics.circle( "fill", th / math.pi - 1.0, r, 0.01)
--Second derivative.
love.graphics.setColor( 0, 0, 1.0, 0.7 )
r = cur:SecondDerivative( th )
x, y = r * math.cos( th ), r * math.sin( th )
--love.graphics.circle( "fill", x, y, 0.01 )
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 )
end
end
love.graphics.setColor( 1, 1, 1, 0.5 )
local r = cur:Interpolate( t )
local x, y = r * math.cos( t ), r * math.sin( t )
love.graphics.circle( "fill", x, y, 0.02 )
end
local function OnImpact( impact )
end
local function Update()
--Deep copy of current state to old state.
for name, t in pairs( cur ) do
for i = 1, N do
old[name][i] = t[i]
end
end
--Deep copy of new state to current state.
for name, t in pairs( new ) do
for i = 1, N do
cur[name][i] = t[i]
end
end
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 ) )
end
cur:DFT( )
end
local function DetectCollision( px, py, vpx, vpy )
local xi, xf = px + vpx * step, py + vpy * step
end
local function Reset()
old = Wave()
cur = Wave()
new = Wave()
end
Reset()
return {
Reset = Reset,
Update = Update,
Integrate = Integrate,
OnImpact = OnImpact,
DetectCollision = DetectCollision,
Draw = Draw,
Current = Current,
Next = Next,
}