--Render and simulate 1D wave equation. local love = love local N = 33 local SOUNDSPEED local IMPULSESIZE local DAMPING 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 local shader = love.graphics.newShader([[ uniform float re[33]; uniform float im[33]; uniform float score; //Slow IDFT float r( float x ) { float r = re[0]; for( int k = 1; k < 17; k++ ) { float c = cos( x * float( k ) ); float s = sin( x * float( k ) ); r += + c * re[k] - s * im[k] + c * re[33 - k] + s * im[33 - k]; } return r; } vec4 effect(vec4 color, Image tex, vec2 texture_coords, vec2 screen_coords) { vec2 p = (3.0 * screen_coords - 1.5 * love_ScreenSize.xy ) / love_ScreenSize.y; p.y = -p.y; float r = r( atan(p.y, p.x) ) - length( p ); float q = float( r < 0.05 ) * clamp( 0.5 - score, 0.0, 1.0 ) ; return vec4( q + (1.0 + clamp( score, 0.0, 1.0 ) * r * r * 0.2) * color.xyz, float(r > 0.0) ) ; } ]]) --Calculate discrete fourier transform of radius function. 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. 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. 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. 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, adjust free parameters according to game state. local function OnImpact( impact, level ) IMPULSESIZE = 10.0 * ( level - 2.0) / 120.0 SOUNDSPEED = 25 - 10 * level / 120 DAMPING = 0.01 * ( 1.0 - 0.4 * level / 120 ) --Apply bandlimited impulse. 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 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 = { 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 = {}, --SPACE DFT of radius function (which is periodic) dftre = {}, dftim = {}, } for i = 1, N do t.radii[i] = 1.0 end DFT( t ) return setmetatable(t, mt) end local function Draw( score ) -- Blue circle. love.graphics.setColor( 91 / 255, 206 / 255, 250 / 255 ) shader:send( "re", unpack( cur.dftre ) ) shader:send( "im", unpack( cur.dftim ) ) shader:send( "score", score ) love.graphics.setShader( shader ) love.graphics.circle("fill", 0, 0, 1.5) local t = love.timer.getTime() love.graphics.setShader( ) -- 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 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 Integrate = function( step ) for i = 1, N do local rxx = cur:SecondDerivative( math.pi * 2.0 * ( i - 1 ) / N ) local r = ( 1.0 - DAMPING ) * ( 2.0 * cur.radii[i] - old.radii[i] + step * step * SOUNDSPEED * rxx ) --Verlet + DAMPING --Damping: oscillate toward 1. --Avoid explosions. r = math.max( 0.5, math.min( r, 4.0 ) ) new.radii[i] = r end new:DFT( ) end local function DetectCollision( px, py, vpx, vpy ) local xi, xf = px + vpx * step, py + vpy * step end local function Reset() IMPULSESIZE = 1 / 10.0 SOUNDSPEED = 5.0 DAMPING = 0.1 / 1 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, }