using PyPlot
xs = -5:0.01:5
figure(figsize=(5,3))
plot(xs, sin.(xs))
grid(true)

# 簡単な数式処理をできる

using SymPy
using PyPlot
ae()=PyPlot.axes()[:set_aspect]("equal") #プロットのアスペクト比を1:1にする函数

# 楕円曲線暗号で有名なEdwards曲線をプロット

x, y = symbols("x y") # 数式処理用に変数 x,y を準備
edwards=x^2+y^2-1+20x^2*y^2 # Edwards 曲線の定義式
display(edwards) # display函数で綺麗に数式を表示
plot_implicit(edwards,(x,-1.5,1.5),(y,-1.5,1.5))
ae() # この行をコメントアウトするとプロットの縦横比が1:1でなくなる

# 変な曲線

x, y = symbols("x y")
n=5
m=n+2
f=x*prod([x^2-i^2 for i in 1:n])
f-=2*n*subs(f,x,y)
display(f)
plot_implicit(f,(x,-m,m),(y,-m,m))
ae()

# 不定積分

x, y = symbols("x y")
integrate(x^y,x)

# 不定積分

x = symbols("x")
integrate(1/sqrt(1-x^2),x)

# 定積分

x = symbols("x", real=true)
integrate(1/(1+x^2),(x,0,oo))

# ζ(12)

n = symbols("n", integer=true)
summation(1/n^12,(n,1,oo))

# 複素函数のベクトル場表示

using PyPlot
plotvf(z,w) = plt[:quiver](real(z),imag(z),real(w)./abs.(w),-imag(w)./abs.(w)) #複素函数z→wのベクトル場をプロット
plotstream(z,w) = plt[:streamplot](real(z),imag(z),real(w),-imag(w)) #複素函数z→wの流れをプロット

x = [-2:0.2:2;]' # xは縦ベクトルでなければいけない
y = x' # yは横ベクトル 
z = x .+ y.*im # imは虚数単位、このように足すと自動的に正方行列になる
w = (-1-im)./z # w=(-1-i)/z
plotvf(z,w) # ベクトル場の作図
plotstream(z,w) # 流れの作図
PyPlot.axes()[:set_aspect]("equal") # アスペクト比を1:1に

x = [-2.5:0.2:2.5;]' # xは縦ベクトルでなければいけない
y = [-1.5:0.2:1.5;]  # yは横ベクトル 
z = x .+ y.*im # imは虚数単位、このように足すと自動的に正方行列になる
w = @. 1/(z+1) - 1/(z-1) # z=-1に正の電荷、z=-1に負の電荷
plotvf(z,w) # ベクトル場の作図
plotstream(z,w) # 流れの作図
PyPlot.axes()[:set_aspect]("equal") # アスペクト比を1:1に

# 3D plot

using PyPlot
using Distributions
linspace(start, stop, length) = range(start, stop=stop, length=length)

a=[
    2. -1.
    -1. 2.
]
n=101
x=linspace(-5,5,n) # -5:0.1:5 と同じ
y=linspace(-5,5,n)
z=zeros(n,n)
for i in 1:n
    for j in 1:n
        z[i,j]=pdf(MvNormal(a),[x[i];y[j]])
    end
end

surf(x, y, z, cmap=ColorMap("CMRmap_r"))

# GIGアニメーションファイル test.gif の作成
# Imagemagickのconvertコマンドを使う

using PyCall
#@pyimport matplotlib.animation as anim
using PyPlot
using Base64
using Printf

function showgif(filename)
    open(filename) do f
        base64_video = base64encode(f)
        display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
    end
end

fig = figure()
ax = PyPlot.axes()

x = [-1.5:0.05:1.0;]
y = [-1.0:0.05:1.0;]
z = x' .+ im*y
w0 = sqrt.(z.+1)./z

# plotw has period 1.
# plotw(t+1) == plotw(1)
function plotw(t)
    clf()
    title("\$t = $t\$")
    xlim(-1.4,0.9)
    ylim(-0.9,0.9)
    w = exp(pi*im*t)*w0
    streamplot(x,y,real(w),-imag(w))
    PyPlot.axes()[:set_aspect]("equal")
end

n=10 # frame/loop
l=1 # loop
interval=200 # milli seconds
filename = "test" # filename

withfig(fig) do
    for k in 1:n*l
        plotw(k/n)
        savefig(filename*@sprintf("_%04d",k), bbox_inches="tight")
    end
end

run(`convert -delay $(interval/10) -loop 0 $(filename)_\*.png $filename.gif`)
#run(`magick -delay $(interval/10) -loop 0 $(filename)_\*.png $filename.gif`)

showgif("$filename.gif")

repmat(a, m, n) = repeat(a, inner=(m,n))
using PyPlot
pygui(false)
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

r = 1
eps = 0.0001
h = 0.01
xx = [-r+eps:h:r+eps;]
yy = [-r+eps:h:r+eps;]
x, y = meshgrid(xx, yy)

z = @. x^2*y/(x^2+y^2)

# contour plot
figure(figsize=(6,5))
p1 = contour(x,y,z,[-0.5:0.05:0.5;])
title("Contour Plot")
xlabel("X")
ylabel("Y")
clabel(p1, inline=1, fontsize=6)
colorbar()
PyPlot.axes()[:set_aspect]("equal")

# 3D plot
figure(figsize=(6,6))
p2 = surf(xx,yy,z, cmap=ColorMap("rainbow"))
title("3D Plot")
xlabel("X")
ylabel("Y")

# Open interactive window

pygui(true)
figure(figsize=(6,6))
surf(xx,yy,z, cmap=ColorMap("rainbow"))

linspace(start, stop, length) = range(start, stop=stop, length=length)
repmat(a, m, n) = repeat(a, inner=(m,n))
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

r = 5.0
xx = linspace(-r,r,100)
yy = linspace(-r,r,100)
x, y = meshgrid(xx, yy)

a = 1
b = -20
c = 2

z = a*x.^2 + 2*b*x.*y + c*y.^2

pygui(false)
#pygui(true)

figure(figsize=(8,6))
imshow(z, extent=(-r,r,-r,r), cmap=ColorMap("rainbow"))
p3 = contour(x, y, z, cmap=ColorMap("gray"))
title("\$z=ax^2+2bxy+cy^2\$,   \$(a,b,c)=($a,$b,$c)\$")
clabel(p3, inline=1, fontsize=6)
colorbar()
PyPlot.axes()[:set_aspect]("equal")

figure(figsize=(8,7))
p1 = surf(x, y, z, cmap=ColorMap("rainbow"))
title("\$z=ax^2+2bxy+cy^2\$,   \$(a,b,c)=($a,$b,$c)\$")
xlabel("x")
ylabel("y")