using Plots
gr(); ENV["PLOTS_TEST"] = "true"
#clibrary(:colorcet)
clibrary(:misc)

function pngplot(P...; kwargs...)
    sleep(0.1)
    pngfile = tempname() * ".png"
    savefig(plot(P...; kwargs...), pngfile)
    showimg("image/png", pngfile)
end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))

showimg(mime, fn) = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:$mime;base64,$base64">""")
end

using SymPy
#sympy[:init_printing](order="lex") # default
#sympy[:init_printing](order="rev-lex")

using SpecialFunctions
using QuadGK

# e^{ipx} を sin(px) に置き換えた場合の図

a = 0.0
b = 10.0
f(p,x) = sin(p*x)
maxa(a,b,p) = a + 2π/abs(p)*fld(b-a, 2π/abs(p))
x = a:0.01:b
PP = []
for p in [5, 20]
    xx = maxa(a,b,p):0.001:b
    P = plot(title="p = $p", titlefontsize=10)
    plot!(legend=false, ylims=(-1.05,1.05))
    plot!(x, f.(p,x))
    hline!([0], color=:cyan, ls=:dot)
    plot!(xx, f.(p,xx), color=:red, fill=(0, 0.5, :red))
    push!(PP, P)
end
plot(PP..., size=(700, 250))

# ディリクレ核 D_N(x) のグラフ

D(N,x) = iszero(x) ? 2N : sin(2π*N*x)/(π*x)
PP = []
for N in [1,2,3,4,5,6]
    x = -6:0.01:6
    P = plot(x, D.(N,x), title="N=$N", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

plot(PP[4:6]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Δ_N(p) のグラフ
#
# N を大きくすると足し上げる p の範囲が拡大する.

Δ(N,p) = max(1-abs(p)/N, zero(p))
PP = []
for N in [1,3,8]
    p = -10:0.01:10
    P = plot(p, Δ.(N,p), title="N=$N", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Fejér核のグラフ
#
# Dirichlet核よりも x=0 から離れた部分の値が急速に 0 に収束することがわかる.

F(N,x) = iszero(x) ? N : (sin(π*N*x)/(π*x))^2/N
PP = []
for N in [1,2,3,4,5,6]
    x = -6:0.01:6
    P = plot(x, F.(N,x), title="N=$N", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

plot(PP[4:6]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Phat_t(p) のグラフ
#
# t > 0 を小さくすると足し上げる p の範囲が拡大する.

Phat(t,p) = exp(-2π*t*abs(p))
PP = []
for t in [0.5, 0.1, 0.05]
    p = -10:0.01:10
    P = plot(p, Phat.(t,p), title="t=$t", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Poisson核のグラフ
#
# t > 0 を小さくすると x 軸上での広がりの範囲が狭くなる.

Poisson(t,x) = t/(t^2+x^2)/π
PP = []
for t in [0.5, 0.1, 0.05]
    x = -1:0.001:1
    P = plot(x, Poisson.(t,x), title="t=$t", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

x = symbols("x", real=true)
p = symbols("p", real=true)
t = symbols("t", positive=true)
W(t,p) = simplify(integrate(exp(-x^2/(2t))*exp(-2*Sym(π)*im*p*x), (x,-oo,oo)))/sqrt(2*Sym(π)*t)
W(t,p)

G(t,x) = simplify(integrate(W(t,p)*exp(2*Sym(π)*im*p*x), (p,-oo,oo)))
G(t,x)

integrate(G(t,x), (x,-oo,oo))

[
    factor(diff(G(t,x), x, x))/2
    factor(diff(G(t,x), t))
]

# Ghat_t(p) のグラフ
#
# t > 0 を小さくすると足し上げる p の範囲が拡大する.

Ghat(t,p) = exp(-2t*π^2*p^2)
PP = []
for t in [0.1, 0.01, 0.004]
    p = -10:0.01:10
    P = plot(p, Ghat.(t,p), title="t=$t", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Gauss核のグラフ
#
# t > 0 を小さくすると x 軸上での広がりの範囲が狭くなる.

G(t,x) = exp(-x^2/(2t))/√(2π*t)
PP = []
for t in [0.1, 0.01, 0.004]
    x = -1:0.001:1
    P = plot(x, G.(t,x), title="t=$t", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# フーリエ級数のディリクレ核 D_N(x) のグラフ

D(N,x) = iszero(x) ? 2N+1 : sin(π*(2N+1)*x)/(sin(π*x))
PP = []
for N in [1:6; 10; 20; 30]
    x = -0.5:0.001:0.5
    P = plot(x, D.(N,x), title="N=$N", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

plot(PP[4:6]..., size=(750, 150), legend=false, layout=@layout([a b c]))

plot(PP[7:9]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# |sin πx| ≤ π|x| の確認

f(x) = abs(sin(π*x))
g(x) = π*abs(x)
x = -1/2:0.002:1/2
plot(size=(400, 250), legend=:top)
plot!(x, f.(x), label="|sin(pi x)|")
plot!(x, g.(x), label="pi |x|")

# |D_N(x)| の積分の漸近挙動の確認

D(N,x) = iszero(x) ? π*(2N+1) : sin(π*(2N+1)*x)/sin(π*x)
F(N) = quadgk(x-> abs(D(N,x)), 0, 1)[1]
G(N) = (2/π)^2*log(2N+1)+1
n = [1,2,3,6,10,20,30,60,100,200,300,600]
plot(size=(400,250), legend=:topleft, xlabel="N", xscale=:log)
plot!(n, F.(n), label="integal of |D_N(x)| from x=0 to 1", lw=2)
plot!(n, G.(n), label="(2/pi)^2*log(2N+1)+1", lw=2, ls=:dash)

# (1/π)∫_0^{aπ} sinc(x) dx のグラフ

f(a) = quadgk(x->sin(x)/x, eps(), a*π)[1]/π
a = 0.2:0.01:6
plot(size=(400,250), legend=false, xlims=(0,maximum(a)))
plot!(title="(1/pi) int_0^{a pi} sin(x)/x dx", titlefontsize=11)
plot!(xlabel="a")
plot!(a, f.(a))

# のこぎり波のグラフ

sawtooth(x) = 1/2 - (x - floor(x))
x = -2:0.01:1.99
plot(x, sawtooth.(x), legend=false, size=(400, 150))
hline!([0], ls=:dot, color=:black)

# のこぎり波のGibbs現象
#
# Fourier級数の部分和が不連続点で上下にオーバーシュートしている.

sawtooth(x) = 1/2 - (x - floor(x))
S(N,x) = sum(n->sin(2π*n*x)/(π*n), 1:N)
PP = []
for N in [10, 20, 50, 100]
    x = -1:0.0005:0.99
    P = plot(legend=false, ylims=(-0.61, 0.61))
    plot!(x, sawtooth.(x))
    plot!(x, S.(N, x))
    plot!(title="N=$N", titlefontsize=10)
    push!(PP, P)
    
    x = -0.01:0.0005:0.2
    P = plot(legend=false, ylims=(0.28, 0.61))
    plot!(x, sawtooth.(x))
    plot!(x, S.(N, x))
    plot!(title="N=$N", titlefontsize=10)
    push!(PP, P)
end

plot(PP[1:2]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[3:4]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[5:6]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[7:8]..., size=(500, 200), layout=@layout([a b{0.3w}]))

# 矩形波のGibbs現象
#
# Fourier級数の部分和が不連続点で上下にオーバーシュートしている.

squarewave(x) = (x - floor(x) ≤ 1/2) ? one(x) : -one(x)
S(K,x) = 4*sum(k->sin(2π*(2k-1)*x)/(π*(2k-1)), 1:K)
PP = []
for K in [5, 10, 25, 50]
    x = -1:0.0005:0.99
    P = plot(legend=false, ylims=(-1.22, 1.22))
    plot!(x, squarewave.(x))
    plot!(x, S.(K, x))
    plot!(title="K=$K", titlefontsize=10)
    push!(PP, P)
    
    x = -0.01:0.0005:0.2
    P = plot(legend=false, ylims=(0.7, 1.22))
    plot!(x, squarewave.(x))
    plot!(x, S.(K, x))
    plot!(title="K=$K", titlefontsize=10)
    push!(PP, P)
end

plot(PP[1:2]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[3:4]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[5:6]..., size=(500, 200), layout=@layout([a b{0.3w}]))

plot(PP[7:8]..., size=(500, 200), layout=@layout([a b{0.3w}]))

function f(x)
    xx = x - floor(x)
    xx == 0 && return 1/4
    0 < xx < 1/4 && return 1/2
    xx == 1/2 && return 0.0
    1/4 < xx < 1/2 && return -1/2
    xx == 1/2 && return -1/4
    1/2 < xx && return 0.0
end
S(M,N,x) = sum(m->(-1)^(m-1)*cos(2π*(2m-1)x)/(π*(2m-1)), 1:M) + sum(n->sin(2π*(4n-2)*x)/(π*(2n-1)), 1:N)

PP = []
for N in [15, 45]
    M = N
    x = -1.0:0.0005:0.9995
    P = plot(legend=false)#, ylims=(-0.61, 0.61))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)

    x = 0.6:0.0005:0.9
    P = plot(legend=false, ylims=(-0.2, 0.2))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)

    x = -0.01:0.0005:0.26
    P = plot(legend=false, ylims=(0.35, 0.61))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)
end

plot(PP[1:3]..., size=(640, 200), layout=@layout([a b{0.2w} c{0.3w}]))

plot(PP[4:6]..., size=(640, 200), layout=@layout([a b{0.2w} c{0.3w}]))

PP = []
for N in [10, 30]
    M = 2N
    x = -1.0:0.0005:0.9995
    P = plot(legend=false)#, ylims=(-0.61, 0.61))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)

    x = 0.6:0.0005:0.9
    P = plot(legend=false, ylims=(-0.2, 0.2))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)

    x = -0.01:0.0005:0.26
    P = plot(legend=false, ylims=(0.35, 0.61))
    plot!(x, f.(x))
    plot!(x, S.(M, N, x))
    plot!(title="(M,N)=($M,$N)", titlefontsize=10)
    push!(PP, P)
end

plot(PP[1:3]..., size=(640, 200), layout=@layout([a b{0.2w} c{0.3w}]))

plot(PP[4:6]..., size=(640, 200), layout=@layout([a b{0.2w} c{0.3w}]))

# Fourier級数のFejer核のプロット

Fejer(N,x) = 1/(N+1)*(sin(π*(N+1)*x)/sin(π*x))^2
PP = []
for N in [10; 20; 30]
    x = -0.5:0.001:0.5
    P = plot(x, Fejer.(N,x), title="N=$N", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

x = -1/2:0.005:1/2
f(x) = 1 - 8*x^2
P = plot(size=(400,230), legend=:bottom)
plot!(x, cos.(2π.*x), label="cos(2pi y)")
plot!(x, f.(x), label="1 - 8y^2")

Poisson(r,x) = (1-r^2)/(1-2r*cos(2π*x)+r^2)
g(r,x) = 1/8*(1-r)/x^2
x = -1/2:0.005:1/2
r = 0.8
plot(size=(400,230), ylim=(0, 1.01*(1+r)/(1-r)))
plot!(title="r = $r", titlefontsize=10)
plot!(x, Poisson.(r, x), label="P_r(x)")
plot!(x, g.(r,x), label="(1/8)(1-r)/x^2")

# Fourier級数のPoisson核のプロット

Poisson(r,x) = (1-r^2)/(1-2r*cos(2π*x)+r^2)
PP = []
for r in [0.3; 0.8; 0.9]
    x = -0.5:0.001:0.5
    P = plot(x, Poisson.(r,x), title="r = $r", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

# Fourier級数のGauss核のプロット

Gauss(t,x; N=100) = (1/√t)*sum(m->e^(-π*(x-m)^2/t), -N:N)
PP = []
for t in [0.1; 0.01; 0.004]
    x = -0.5:0.001:0.5
    P = plot(x, Gauss.(t,x), title="t = $t", titlefontsize=10)
    push!(PP, P)
end
plot(PP[1:3]..., size=(750, 150), legend=false, layout=@layout([a b c]))

n = big"1000000"
Float64(factorial(n)/(n^n*e^(-n)*√n)), √(2π)

u = symbols("u", positive=true)
[integrate(k/(1+u^k), (u, 0, oo)) for k in 2:6]

x = symbols("x")
f(x) = sin(π*x)/(π*x)
series(f(x), x, n=10)

x = symbols("x")
f(x) = sin(π*x)/(π*x)
series(f(x)*f(im*x), x, n=12)

x = symbols("x")
f(x) = sin(π*x)/(π*x)
ω = exp(im*Sym(π)/4)
series(f(x)*f(ω*x)*f(ω^2*x)*f(ω^3*x), x, n=24)

x = symbols("x")
f(x) = sin(π*x)/(π*x)
ω = exp(im*Sym(π)/8)
@time series(prod(f(ω^k*x) for k in 0:7), x, n=48)

BernoulliNumber(n) = sympy[:bernoulli](n)
[BernoulliNumber(n) for n in 0:8]

BernoulliNumber(n) = sympy[:bernoulli](n)
B = [BernoulliNumber(2k) for k in 1:8]

BernoulliNumber(n) = sympy[:bernoulli](n)
Z = [2^(2k-1)*(-1)^(k-1)*BernoulliNumber(2k)*PI^(2k)/factorial(2k) for k in 1:8]

[zeta(Sym(2k)) for k in 1:8]

Theta(t; N=20) = iszero(N) ? one(t) : 1 + 2*sum(n->e^(-π*n^2*t), 1:N)
t = 0.2:0.05:2.0
A = Theta.(t)
B = 1./sqrt.(t).*Theta.(1./t)
plot(size=(500, 350))
plot!(t, B, label="Theta(t)",               lw=2)
plot!(t, A, label="(1/sqrt(t)) Theta(1/t)", lw=2, ls=:dash)

Theta(t; N=20) = iszero(N) ? one(t) : 1 + 2*sum(n->e^(-π*n^2*t), 1:N)
t = 0.1
@show 1/√t*Theta(1/t; N=0)
@show 1/√t*Theta(1/t; N=1)
@show 1/√t*Theta(1/t; N=2)
@show Theta(t; N=9);
@show Theta(t; N=10);
@show Theta(t; N=11);