using PyPlot
x = range(0,1,length=200)
plot(x, sin.(π*x), "-")
plot(x, sin.(2π*x), "-")
plot(x, sin.(3π*x), "-")
plot(x, sin.(4π*x), "-")
legend([L"\sin(\pi x)", L"\sin(2\pi x)", L"\sin(3\pi x)", L"\sin(4\pi x)"])
xlabel(L"x")
title("orthogonal sine functions")

using QuadGK

sinecoef(f, m) = 2 * quadgk(x -> f(x) * sin(m*π*x), 0,1, atol=1e-8 * sqrt(quadgk(x->abs2(f(x)),0,1)[1]))[1]

f(x) = 0.5 - abs(x - 0.5)
sinecoef.(f, 1:20)

using PyPlot
n = 1:2:99 # odd integers from 1 to 99
loglog(n, abs.(sinecoef.(f, n)), "bo")
xlabel(L"n")
ylabel(L"Fourier sine coefficient $|b_n|$")
title(L"Sine series convergence for $f(x) = 0.5 - |x - 0.5|$")
loglog(n, 1 ./ n.^2, "b--")
legend([L"|b_n|", L"1/n^2"])

# First, define a function to evaluate N terms of the sine series, given the coefficients b
function sinesum(b, x)
    f = 0.0
    for n = 1:length(b)
        f += b[n] * sin(n*π*x)
    end
    return f
end

#using Interact
fig = figure()
x = range(0,1, 1000)
#@manipulate for n=1:2:99
for n = 1:2:25
    display(
    withfig(fig) do
        plot(x, f.(x), "k--")
        b = sinecoef.(f, 1:n)
        plot(x, [sinesum(b, x) for x in x], "r-")
        xlabel(L"$x$")
        legend(["exact f", "$n-term sine series"])
    end
    )
end

fig = figure(figsize=(15,6))
g(x) = sin(sin(3π*x) + 5*sin(π*x))
# @manipulate for n=1:2:37
for n=1:2:17
    display(
    withfig(fig) do
        subplot(1,2,1)
        plot(x, g.(x), "k--")
        b = sinecoef.(g, 1:n)
        plot(x, [sinesum(b, x) for x in x], "r-")
        xlabel(L"x")
        legend(["exact g", "$n-term sine series"])
        
        subplot(1,2,2)
        semilogy(1:2:n, abs.(b[1:2:n]), "bo-")
        xlabel(L"n")
        ylabel(L"coefficient $|b_n|$")
    end
    )
end

using Interact

fig = figure()
@manipulate for t in slider(0:0.001:1, value=0.0, label="time t")
    withfig(fig) do
        b = @. sinecoef(f, 1:199) * exp(-((1:199)*π)^2 * t)
        plot(x, [sinesum(b, x) for x in x])
        xlabel(L"x")
        title("diffusion solution at time $t")
        ylim(0,0.5)
    end
end

fig = figure()
@manipulate for t in slider(0:0.01:10, value=0.0, label="time t")
    withfig(fig) do
        b = @. sinecoef(f, 1:199) * cos(((1:199)*π) * t)
        plot(x, [sinesum(b, x) for x in x])
        xlabel(L"x")
        title("wave solution at time $t")
        ylim(-0.5,0.5)
        grid()
    end
end