--Render and simulate 1D wave equation.
local love = love
local N = 33
local SOUNDSPEED = 0.5
local IMPULSESIZE = 5
local DAMPING = 0.01


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 );

    return vec4( (1.0 + float(score > 1.0) * r * 0.1) * color.xyz, float(r > -0.01)) ;
}
  ]])


--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.
local function OnImpact( impact )

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

end

local function DetectCollision( px, py, vpx, vpy )
  local xi, xf = px + vpx * step, py + vpy * step
end

local function Reset()
  
  SOUNDSPEED = 0.5
  IMPULSESIZE = 5
  DAMPING = 0.02

  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,
}