using Plots;
using QuadGK;
using DelimitedFiles;
using Interpolations;

struct FourierSeries
    f::Function
    a0::Float64
    ans::Array{Float64,1}
    bns::Array{Float64,1}
end

"""
    FourierSeries(f; N=10)

Generate a FourierSeries object containing the the first `N` (excluding a0)
Fourier coefficients of the function `f` over the interval (-π, π).
"""
function FourierSeries(f::Function; N::Int = 10)
    # Compute ans
    a0::Float64 = quadgk(t -> f(t), -π, π)[1] / π
    ans = Array{Float64,1}(undef, N)
    bns = Array{Float64,1}(undef, N)
    for n = 1:N
        ans[n] = quadgk(t -> f(t) * cos(n * t), -π, π)[1] / π
        bns[n] = quadgk(t -> f(t) * sin(n * t), -π, π)[1] / π
    end
    FourierSeries(f, a0, ans, bns)
end

"""
  eval(t::Float64, fs::FourierSeries)

Evaluate `FourierSeries` object at time `t`.
"""
function eval(t::Float64, fs::FourierSeries)
    fs.a0 / 2 +
    sum(fs.ans[n] * cos(n * t) + fs.bns[n] * sin(n * t) for n = 1:length(fs.ans))
end

"""
  eval(ts::AbstractRange, fs::FourierSeries)

Evaluate `FourierSeries` object at times `ts`.
"""
function eval(ts::AbstractRange, fs::FourierSeries)
    [eval(t, fs) for t in ts]
end

"""
  circle(x,y,r[; npts=500])

Compute x and y coordinates of a circle of radius `r` centered at (`x`, `y`).
"""
function circle(x, y, r; npts=500)
    θ = LinRange(0.0, 2π, npts)
    x .+ r * sin.(θ), y .+ r * cos.(θ)
end

"""
  saw_tooth(t)

Returns values of the saw tooth waveform
"""
saw_tooth(t) = t

"""
  square_wave(t)

Returns values of the square_wave waveform
"""
square_wave(t) = t < 0 ? 0.0 : π

"""
Make the plot!
"""
fs = FourierSeries(saw_tooth)
@gif for i = 1:1500
    shift = i / 100
    plot(
        ts,
        [eval(t + shift, fs) for t in ts],
        ylims = (-1.5, π + 1.5),
        legend = false,
        xlims = (-π, 4π),
        framestyle = :box,
        aspect_ratio = 1
    )
    x = 2.5π
    y = fs.a0 / 2
    r = sqrt(fs.ans[1]^2 + fs.bns[1]^2)
    ϕ = -atan(fs.bns[1], fs.ans[1])
    plot!(
        circle(x, y, r),
        seriestype = [:shape,],
        linecolor = :black,
        c = :white,
        legend = false
    )
    for n = 2:length(fs.ans)
        x += r * sin((n - 1) * (ts[end] + shift) + ϕ)
        y += r * cos((n - 1) * (ts[end] + shift) + ϕ)
        r = sqrt(fs.ans[n]^2 + fs.bns[n]^2)
        ϕ = -atan(fs.bns[n], fs.ans[n])
        plot!(
            circle(x, y, r),
            seriestype = [:shape,],
            linecolor = :black,
            c = :white,
            legend = false,
            fillalpha = 0.0
        )
        scatter!([x], [y], markersize = 0.3)
    end
    x += r * sin(length(fs.ans) * (ts[end] + shift) + ϕ)
    y += r * cos(length(fs.ans) * (ts[end] + shift) + ϕ)
    scatter!([x], [y], markersize = 0.1)
    interps = range(ts[end], stop = x, length = 50)
    plot!(interps, [y for _ in interps], color = :purple)
end every 10