--Render and simulate 1D wave equation.
local love = love
local step = assert( step )

local N = 25

--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) / 13
    radii = {},
    --TIME derivative of radius
    vrad = {},

    --SPACE DFT of radius function (which is periodic)
    dftre = {},
    dftim = {},
  }

  for i = 1, N do
    t.radii[i] = 1.0
    t.vrad[i] = 0.0
  end
  DFT( t )

  return setmetatable(t, mt)
end


local old = Wave() --State at beginning of tick.
local cur = Wave() --State at end of tick.


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.5 )
      th = ( i - 1 + k ) * 2.0 * math.pi / N
      r = cur:Interpolate( th )
      x, y = r * math.cos( th ), r * math.sin( th )
      love.graphics.circle( "fill", x, y, 0.01 )
      
      --First derivative.
      love.graphics.setColor( 0, 1.0, 0, 0.5 )
      r = 1.0 + cur:Derivative( th )
      x, y = r * math.cos( th ), r * math.sin( th )
      love.graphics.circle( "fill", x, y, 0.01 )
      
      --Second derivative.
      love.graphics.setColor( 0, 0, 1.0, 0.5 )
      r = 1.0 + cur:SecondDerivative( th )
      x, y = r * math.cos( th ), r * math.sin( th )
      love.graphics.circle( "fill", x, y, 0.01 )
    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()
  local t = love.timer.getTime()
  for i = 1, N do
    local rxx = old:SecondDerivative( 2.0 * math.pi * (i - 1) / N )
    cur.vrad[i] = 0.7 * old.vrad[i] + 0.5 * step * rxx
    cur.radii[i] = old.radii[i] + step * cur.vrad[i]
  end
  
  cur:DFT( )
  
  if love.math.random( 120 ) < 3 then
    cur:Impulse( 0, 0.2 )
  end
  
  --Deep copy of current state to old state.
  for name, t in pairs( cur ) do
    for i = 1, 13 do
      old[name][i] = t[i]
    end
  end

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()
end

Reset()

return {
  Reset = Reset,
  Update = Update,
  DetectCollision = DetectCollision,
  Draw = Draw,
}