using DifferentialEquations
using Plots

"""
Operador de diferenças finitas para trás, com condições de contorno periódicas
"""
function δ⁻(u::Vector{Float64}, h::Float64, ::Val{:per})
    du = zero(u)
    for j = 2:length(u)
        du[j] = (u[j] - u[j-1])/h
    end
    du[1] = (u[1] - u[end])/h
    return du
end

function dudt_advection⁻!(dudt::Vector{Float64}, u::Vector{Float64},
        p::Vector{Float64}, t::Float64)
    a, h = p
    du = δ⁻(u, h, Val(:per))
    dudt .= - a * du
    return nothing
end

a = 0.5 # # velocidade de propagação
L = 2π # comprimento do intervalo [0,L]
N = 60 # número de pontos da malha
h = L/(N-1) # comprimento de cada partição na malha
x = range(0.0, L, length=N) # discretização do intervalo [0,L] com N pontos, incluindo os extremos
u₀ = sin.(x) # condição inicial, com N pontos (incluindo os dois extremos)
p = [a, h] # parâmetros
Tf = L/a # tempo final (volta completa)
τ = 0.1 # intervalo de tempo para guardar a solução
tspan = (0.0,Tf) # intervalo de tempo
prob = ODEProblem(dudt_advection⁻!, u₀, tspan, p, saveat = τ)
nothing

sol = solve(prob, Tsit5())
sol.retcode

sol

plot(x, sol.u[end], ylims=(-1.1, 1.1), label="solução no tempo final $(sol.t[end])")

anim = @animate for (t,u) in zip(sol.t, sol.u)
    plot(x, u, ylims=(-1.1, 1.1), label="t=$(round(t,digits=2))",
        title="Uma solução da equação de adveção com a=$a", titlefont=10)
end
gif(anim, "../../../assets/attachments/img/anim_adveccao_dir.gif", fps = 20)
nothing

"""
Operador diferenças finitas para frente, com condições de contorno periódicas.
"""
function δ⁺(u, h, ::Val{:per})
    du = zero(u)
    for j = 1:length(u)-1
        du[j] = (u[j+1] - u[j])/h
    end
    du[end] = (u[1] - u[end])/h
    return du
end

function dudt_advection⁺!(dudt, u, p, t)
    a, h = p
    du = δ⁺(u, h, Val(:per))
    dudt .= - a * du
    return nothing
end

a = -0.5 # velocidade de propagação
p = [a, h] # parâmetros
prob = ODEProblem(dudt_advection⁺!, u₀, tspan, p, saveat = τ)
sol = solve(prob, Tsit5())
sol.retcode

anim = @animate for (t,u) in zip(sol.t, sol.u)
    plot(x, u, ylims=(-1.1, 1.1), label="t=$(round(t,digits=2))",
        title="Uma solução da equação de adveção com a=$a", titlefont=10)
end
gif(anim, "../../../assets/attachments/img/anim_adveccao_esq.gif", fps = 20)
nothing