Julia で青海波を描く

Julia で青海波を描く

plotter.jl を include
https://blog.goo.ne.jp/r-de-r/e/bd71a52a09801335d56f7c47d879bfe3

include("plotter.jl")

function seigaiha(nx=6, ny=6; r=1, color=:deepskyblue2, width=720, height=360)
    function unit(x, y)
        for i = 6r:-r:1
            col = i % 2 == 1 ? :white : color
            plotcircle(x, y, i, startangle=0, endangle=180, col=col, fcol=col)
        end
    end
    plotbegin(w=width, h=height)
    for j = 2ny:-1:0
        for i = 0:nx
            unit(i*12r + (j % 2)*6r, j*3r)
        end
    end
    x1, y1 = 0, 3r
    x2, y2 = 12nx*r, (6ny+3)r
    plotlimit(x1, y1, x2, y2)
    println("(width, height) = ($(x2 - x1), $(Int(y2 - y1)))")
    plotend()
    savefig("parts.pdf")
end

seigaiha()
savefig("fig1.png")

Plots で使える描画関数に糖衣をかぶせた自作描画関数群

【描画関数を組み合わせて「麻の葉格子」を描く】で使った自作描画関数群

# ファイル名 plotter.jl

using Plots

# 開始
function plotbegin(; w = 600, h = 600)
    pyplot(grid=false, aspect_ratio=1, size=(w, h),
        xlabel="", ylabel="", label="")
    plt = plot(showaxis=false)
end

# プロット領域（(x1, y1), (x2, y2) を対角座標とする矩形内）の設定
function plotlimit(x1, y1, x2, y2)
    plot!(xlims=(x1, x2), ylims=(y1, y2))
end

# 終了
function plotend()
    plot!()
end

# 座標ベクトル x, y による折線
function plotline(x::Vector{Float64}, y::Vector{Float64};
        col=:black, lty=:solid, lwd=1)
    plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd)
end

# （x1, y1) - (x2, y2) の線分
function plotline(x1, y1, x2, y2; col=:black, lty=:solid, lwd=1)
    plot!([x1, x2], [y1, y2], linecolor=col, linestyle=lty, linewidth=lwd)
end

# （x1, y) - (x2, y) の水平線分
function plothline(x1, x2, y; col=:black, lty=:solid, lwd=1)
    plot!([x1, x2], [y, y], linecolor=col, linestyle=lty, linewidth=lwd)
end

# (x, y1) - (x, y2) の垂直線分
function plotvline(x, y1, y2; col=:black, lty=:solid, lwd=1)
    plot!([x, x], [y1, y2], linecolor=col, linestyle=lty, linewidth=lwd)
end

# (x1, y1), (x2, y2) を対角座標とする矩形
function plotbox(x1, y1, x2, y2; col=:black, lty=:solid, lwd=1, fcol="")
    if fcol == ""
        plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], linecolor=col,
              linestyle=lty, linewidth=lwd)
    else
        plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], linecolor=col,
              linestyle=lty, seriestype=:shape, linewidth=lwd, fillcolor=fcol)
    end
end

# (ox, oy) を原点とする半径 r の円（円弧）。fcol を指定すると塗りつぶし。
function plotcircle(ox, oy, r; startangle=0, endangle=360,
                    col=:black, lty=:solid, lwd=1, fcol="")
    plotellipse(ox, oy, r, r, φ=0, startangle=startangle,
                endangle=endangle, col=col, lty=lty, lwd=lwd, fcol=fcol)
end

# 度をラジアンに変換する
radian(degree) = degree / 180 * π

# (ox, oy) を中心，ra, rb を半径とする楕円（楕円弧）。fcol を指定すると塗りつぶし。
function plotellipse(ox, oy, ra, rb; φ=0, startangle=0, endangle=360,
                     length=100, col=:black, lty=:solid, lwd=1, fcol="")
    θ = vcat(range(radian(startangle), radian(endangle),
              length=length), radian(endangle))
    if φ == 0
        if fcol == ""
            plot!(ra .* cos.(θ) .+ ox, rb .* sin.(θ) .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd)
        else
            plot!(ra .* cos.(θ) .+ ox, rb .* sin.(θ) .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd,
                  seriestype=:shape, fillcolor=fcol)
        end
    else
        x = ra .* cos.(θ)
        y = rb .* sin.(θ)
        φ = radian(φ)
        cosine = cos(φ)
        sine = sin(φ)
        if fcol == ""
            plot!(cosine .* x .- sine .* y .+ ox,
                  sine .* x .+ cosine .* y .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd)
        else
            plot!(cosine .* x .- sine .* y .+ ox,
                  sine .* x .+ cosine .* y .+ oy,
                  linecolor=col, linestyle=lty, linewidth=lwd,
                  seriestype=:shape, fillcolor=fcol)
        end
    end
end

# 任意の多角形。fcol を指定すると塗りつぶし。
function plotpolygon(x, y; col=:black, lty=:solid, lwd=1, fcol="")
    x = vcat(x, x[1])
    y = vcat(y, y[1])
    if fcol == ""
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd)
    else
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd,
              seriestype=:shape, fillcolor=fcol)
    end
end

# (x1, y1) を描き始め，一辺の長さ l の正 n 多角形。fcol を指定すると塗りつぶし。
function plotpolygon1(x1, y1, l, n; col=:black, lty=:solid, lwd=1, fcol="")
    θ = range(0, 2π, length=n+1)
    x = fill(float(x1), n+1)
    y = fill(float(y1), n+1)
    for i = 2:n
        x[i] = x[i-1]+l*cos(θ[i])
        y[i] = y[i-1]+l*sin(θ[i])
    end
    if fcol == ""
        plot!(x, y, linecolor=col, linestyle=lty, linewidth=lwd)
    else
        plot!(x, y, linecolor=col, linestyle=lty,
              seriestype=:shape, fillcolor=fcol, linewidth=lwd)
    end
end

# (ox, oy) を中心，r を半径とする円に内接する正 n 多角形。fcol を指定すると塗りつぶし。
function plotpolygon2(ox, oy, r, n; φ=90, col=:black, lty=:solid, lwd=1, fcol="")
    θ = range(0, 2π, length=n+1) .+ radian(φ)
    if fcol == ""
        plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
              linecolor=col, linestyle=lty, linewidth=lwd)
    else
        plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
              linecolor=col, linestyle=lty, linewidth=lwd,
              seriestype=:shape, fillcolor=fcol)
    end
end

# (x1,y1) - (x2, y2) の矩形範囲に格子を描く。間隔は wx, wy。
function plotgrid(x1, y1, x2, y2, wx; wy=wx, col=:gray, lty=:dot, lwd=1)
    X1 = min(x1, x2)
    X2 = max(x1, x2)
    Y1 = min(y1, y2)
    Y2 = max(y1, y2)
    for i = 0:fld(abs(X2-X1), wx)
        plot!([X1+wx*i, X1+wx*i], [Y1, Y2], linecolor=col,
              linestyle=lty, linewidth=lwd)
    end
    for i = 0:fld(abs(y2-y1), wy)
        plot!([X1, X2], [Y1+wy*i, Y1+wy*i], linecolor=col,
              linestyle=lty, linewidth=lwd)
    end
end

# 二点 (x1, y1), (x2, y2) を通る半径 r の円弧（短い方）を描く。
# 短い方の円弧も，円の中心が 2 つあるので，x1 = x2 の場合は，右か左いずれの中心を使うかを :right, :left で示す。
# x1 ≠ x2 の場合は，上か下いずれの中心を使うかを :upper, :lower で示す。
function plotarc(x1, y1, x2, y2, r; udrl="", col=:black, lwd=1)
    function adjust(φ, x1, y1, x3, y3)
        if x1 >= x3 && y1 >= y3 # 第1象限
            return φ
        elseif x1 >= x3 && y1 < y3 # 第4象限
            return 360-φ
        elseif x1 < x3 && y1 >= y3 # 第2象限
            return 180-φ
        elseif x1 < x3 && y1 < y3 # 第3象限
            return 180+φ
        end
    end
    if y1 > y2
        x1, y1, x2, y2 = x2, y2, x1, y1
    end
    x31 = x32 = y31 = y32 = NaN
    if x1 == x2
        w = r^2 - (y2 - y1)^2 / 4
        if w >= 0
            t0 = sqrt(w)
            x31 = x1 + t0
            x32 = x1 - t0
            y31 = y32 = (y1 + y2) / 2
            θ = asind(abs(y2 - y32) / r)
            if udrl == :right
                plotcircle(x31, y31, r, startangle=180 - θ, endangle=180 + θ, col=col, lwd=lwd)
            elseif udrl == :left
                plotcircle(x32, y32, r, startangle=-θ, endangle=θ, col=col, lwd=lwd)
            else
                error("if x1 == x2, udrl must be ':right' or ':left'")
            end
        end
    else
        t1 = x1^2 - x2^2 + y1^2 - y2^2
        t5 = x1 - x2
        t2 = y2 - y1
        t3 = t5^2 + t2^2
        t4 = t3 - 4r^2
        w = -t3*t4
        if w > 0
            t6 = y1/2 + y2/2
            t7 = sqrt(w)*t5/2t3
            y31 = t6 - t7
            y32 = t6 + t7
            if y31 < y32
                y31, y32 = y32, y31
            end
            x31 = (t1 + 2t2*y31)/2t5
            x32 = (t1 + 2t2*y32)/2t5
            posy1 = y31 - (y2 - y1)/(x2 - x1)*(x31 - x1) - y1
            posy2 = y32 - (y2 - y1)/(x2 - x1)*(x32 - x1) - y1
            if udrl == :upper
                x3, y3 = x31, y31
            else#if udrl == :lower
                x3, y3 = x32, y32
            end
            #println("(x1, y1) = ($x1, $y1),  (x2, y2) = ($x2, $y2),  (x3, y3) = ($x3, $y3)")
            φ = adjust(asind(abs(y1 - y3) / r), x1, y1, x3, y3)
            θ = adjust(asind(abs(y2 - y3) / r), x2, y2, x3, y3)
            θ - φ > 180 && (θ -= 360)
            φ - θ > 180 && (φ -= 360)
            #println("posy1 = $posy1,  posy2 = $posy2,  θ = $θ,  φ = $φ")
            plotcircle(x3, y3, r, startangle=φ, endangle=θ, col=col, lwd=lwd)
        end
    end
    x31, y31, x32, y32
end

描画関数を組み合わせて「麻の葉格子」を描く

日本古来の「麻の葉格子」を描いてみる。

using Plots で使える関数を直接使っても当然描けるが，使い勝手のよいように糖衣をかぶせた関数群を用意している。以前このブログにも掲載したが，その後少々手直しをして使っている。この関数群は plotter.jl というファイルに用意しているので，使うためには include する。using で使えるようにしてもよいのだが，かえって面倒くさいようなので，このようにしている。

include("plotter.jl")

参照：Plots で使える描画関数に糖衣をかぶせた自作描画関数群

なお，以下の出画例は png ファイルであるが，pdf ファイルで保存するとジャギーでない綺麗な（拡大縮小が自在な）図形が描かれると思う。

麻の葉格子は何通りかあるが，一番オーソドックスなのが以下のもの。繰り返しの基本パターンを unit() で定義している。それを二重の for 文で領域に敷き詰めていく。

function asanoha(nx=6, ny=5; r=1, width=600, height=400)
    a = r * sqrt(3) / 2
    b = r / 2
    c = a / 3
    function unit(x, y)
        plotline(x .+ [-a, -a, -c, 0, -a+c, -a, 0, 0], y .+ [0,  b,  b, 0, 0,  b, 0,  b])
        plotline(x .+ [ a,  a,  c, 0,  a-c,  a, 0, 0], y .+ [0,  b,  b, 0, 0,  b, 0,  b])
        plotline(x .+ [-a, -a, -c, 0, -a+c, -a, 0, 0], y .+ [0, -b, -b, 0, 0, -b, 0, -b])
        plotline(x .+ [ a,  a,  c, 0,  a-c,  a, 0, 0], y .+ [0, -b, -b, 0, 0, -b, 0, -b])
    end
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = a, b, 2(nx + 0.5) * a, (2ny + 1) * b
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2)
    for x = 1:nx
        for y = 1:ny
            unit(2a * x, 2b * y)
        end
    end
    plotend()
end

asanoha(4, 3, width=500, height=216)
savefig("fig.png")

2 つ目は，よく似ているが，三菱マークが挿入されているもの。これも基本パターンを unit2() で定義しているのを，二重の for ループで，領域内に敷き詰めていく。

function asanoha2(nx=6, ny=5; r=1, width=600, height=400)
    function unit2(x, y)
        xs = [b, c, 0, a, c]
        ys = [r2, r2, 0, 0, r2]
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:5
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2)
        end
    end
    sqrt3 = sqrt(3)
    a = r / sqrt3
    b = a / 2
    c = 3b
    r2 = r / 2
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = 1.5r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2)
    for y = 1:ny
        for x = 1:nx + 1
          unit2(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * r * 1.5)
        end
    end
    plotend()
end

asanoha2(5, 4, width=500, height=318)
savefig("fig2.png")

3 つ目は，2 つ目の麻の葉格子に円を加えたもの。ちょっとした違いだが，印象はかなり違ってくる。

function asanoha3(nx=6, ny=5; r=1, width=600, height=400)
    function unit2(x, y)
        xs = [b, c, 0, a, c]
        ys = [r2, r2, 0, 0, r2]
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:5
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2)
            plotcircle(x, y, r)
        end
    end
    sqrt3 = sqrt(3)
    a = r / sqrt3
    b = a / 2
    c = 3b
    r2 = r / 2
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = 1.5r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2)
    for y = 1:ny
        for x = 1:nx + 1
          unit2(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * r * 1.5)
        end
    end
    plotend()
end

asanoha3(5, 4, width=500, height=318)
savefig("fig3.png")

R の gmp, Rmpfr が動かない

R で ペル方程式，Pell's equation

にも書いたのだけど，また gmp, Rmpfr が動かなくなっていた。

で，同じように，

install.packages("gmp", type = "source")
install.packages("Rmpfr", type = "source")

で，動くようになった。