裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia で円をモチーフにした格子模様を描く(1)

2021年08月13日 | ブログラミング

Julia で円をモチーフにした格子模様を描く(1)

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

lwd を適当(適切)に変えてください。

include("plotter.jl")

function en(nx=6, ny=5; r=1, lwd=2, width=600, height=360)
    x1, y1, x2, y2 = 4r, r, 5r + (2nx+2)r -3r, (2ny+1)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:2ny+1
        for x = 1:nx+1
          plotcircle((2x+2 + (y % 2))r, y*r, 0.94r, lwd=lwd, col=:cornsilk)
        end
    end
    plotend()
end

en(10, 6, lwd=3, width=600, height=360)
savefig("en1.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で正方形基準の格子模様を描く(3)

2021年08月13日 | ブログラミング

Julia で正方形基準の格子模様を描く(3)

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

include("plotter.jl")

function asanoha7(nx=6, ny=5; a=1, width=600, height=400)
    function unit(x, y, a)
        xy = [0 0; 4a 0; 4a 4a; 0 4a; 0 0; 3a a; 4a 4a; a 3a; 0 0; 4a 4a; NaN NaN;
                  0 4a; a 3a; NaN NaN; 4a 0; 3a a]
        xs = vcat(xy[:, 1], xy[:, 1], -xy[:, 1], -xy[:, 1])
        ys = vcat(xy[:, 2], -xy[:, 2], xy[:, 2], -xy[:, 2])
        plotline(4a*x - 2a .+ xs, 4a*y - 2a .+ ys, lwd=2, col=:bisque)
    end
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = 0, 0, 4a*nx, 4a*ny
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for x = 1:nx
        for y = 1:ny
            unit(x, y, a)
        end
    end
    plotend()
end

asanoha7(6, 4, width=480, height=320)
savefig("fig7.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で正方形基準の格子模様を描く(2)

2021年08月13日 | ブログラミング

Julia で正方形基準の格子模様を描く(2)

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

include("plotter.jl")

function asanoha6(nx=6, ny=5; a=1, width=600, height=400)
    function unit(x, y, a)
        xy = [0 0; 4a 0; 4a 4a; 0 4a; 0 0; 3a a; 4a 4a; a 3a; 0 0; 4a 4a; NaN NaN;
                  0 4a; a 3a; NaN NaN; 4a 0; 3a a]
        xs = vcat(xy[:, 1], xy[:, 1], -xy[:, 1], -xy[:, 1])
        ys = vcat(xy[:, 2], -xy[:, 2], xy[:, 2], -xy[:, 2])
        plotline(8a*x -4a .+ xs, 8a*y - 4a .+ ys, lwd=2, col=:bisque)
    end
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = 0, 0, 8a*nx, 8a*ny
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for x = 1:nx
        for y = 1:ny
            unit(x, y, a)
        end
    end
    plotend()
end

asanoha6(6, 4, width=480, height=320)
savefig("fig6.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で正方形基準の格子模様を描く(1)

2021年08月13日 | ブログラミング

Julia で正方形基準の格子模様を描く(1)

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

include("plotter.jl")

function asanoha5(nx=6, ny=5; a=1, width=600, height=400)
    function unit(x, y, a)
        xy = [0 0; 4a 0; 4a 4a; 0 4a; 0 0]
        if x % 2 != y % 2
            xy = vcat(xy, [3a a; 4a 4a; a 3a; 0 0; 4a 4a; NaN NaN;
                  0 4a; a 3a; NaN NaN; 4a 0; 3a a])
        end
        xs = vcat(xy[:, 1], xy[:, 1], -xy[:, 1], -xy[:, 1])
        ys = vcat(xy[:, 2], -xy[:, 2], xy[:, 2], -xy[:, 2])
        plotline(8a*x-4a .+ xs, 8a*y-4a .+ ys, lwd=2, col=:bisque)
    end
    plotbegin(w=width, h=height)
    x1, y1, x2, y2 = 0, 0, 8a*nx, 8a*ny
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for x = 1:nx
        for y = 1:ny
            unit(x, y, a)
        end
    end
    plotend()
end

asanoha5(6, 4, width=480, height=320)
savefig("fig5.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で麻の葉格子模様を描く(2)

2021年08月13日 | ブログラミング

Julia で麻の葉格子模様を描く(2)

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

include("plotter.js")

function asanoha4(nx=6, ny=5; r=1, width=600, height=400)

    a = r * sqrt(3) / 2
    b = r / 2
    c = a / 3
    d = (a - c) / 2
    e = d/sqrt(3)
    f = a - d
    g = b-e
    h = b - 2e
    function unit(x, y)
        xs = [-a, -a, -c, 0, -a+c, -a, 0, 0, NaN,
            a,  a,  c, 0,  a-c,  a, 0, 0, NaN,
            -a, -a, -c, 0, -a+c, -a, 0, 0, NaN,
            a,  a,  c, 0,  a-c,  a, 0, 0, NaN,
            d, 0,  -d, -d, 0, d, d, NaN,
            f, f, a, NaN, -f, -f, -a, NaN, f, f, a, NaN, -f, -f, -a]
        ys = [0,  b,  b, 0, 0,  b, 0,  b, NaN,
            0,  b,  b, 0, 0,  b, 0,  b, NaN,
            0, -b, -b, 0, 0, -b, 0, -b, NaN,
            0, -b, -b, 0, 0, -b, 0, -b, NaN,
            e, 2e, e, -e, -2e, -e, e, NaN,
            b, g, h, NaN, b, g, h, NaN, -b, -g, -h, NaN, -b, -g, -h]
        plotline(x .+ xs, y .+ ys, lwd=2, col=:bisque)
        plotline(x .+ [-a, a, NaN, -a, a], y .+ [b, -b, NaN, -b, b], lwd=4, col=:bisque)
    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, col=:gray, fcol=:gray)
    for x = 1:nx
        for y = 1:ny
            unit(2a * x, 2b * y)
        end
    end
    plotend()
end

asanoha4(4, 3, width=693, height=300)
savefig("fig4.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で亀甲ベースの格子模様を描く(5)

2021年08月12日 | ブログラミング

Julia で亀甲ベースの格子模様を描く(5)

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

include("plotter.jl")

function kikkou5(nx=6, ny=5; r=1, width=600, height=400)
    function unit(x, y)
        sqrt32 = sqrt(3)/2
        ax, ay = 0, 0
        bx, by = sqrt32*r, r/2
        cx, cy = sqrt32*r, -r/2
        dx, dy = sqrt32/3*r, 0
        ex, ey = sqrt32/2*r, r/4
        fx, fy = 5sqrt32/6*r, r/4
        gx, gy = sqrt32*r, 0
        hx, hy = 5sqrt32/6*r, -r/4
        ix, iy = sqrt32/2*r, -r/4
        xs = [ax, bx, cx, ax, dx, dx, ex, fx, gx, hx, ix, dx, NaN, bx, fx, NaN, cx, hx];
        ys = [ay, by, cy, ay, dy, dy, ey, fy, gy, hy, iy, dy, NaN, by, fy, NaN, cy, hy];
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:length(xs)
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2, lwd=2, col=:cornsilk)
        end
    end
    sqrt3 = sqrt(3)
    x1, y1, x2, y2 = r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:ny
        for x = 1:nx + 1
          unit(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * 1.5r)
        end
    end
    plotend()
end

kikkou5(4, 4, width=780, height=550)
savefig("fig5.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で亀甲ベースの格子模様を描く(4)

2021年08月12日 | ブログラミング

Julia で亀甲ベースの格子模様を描く(4)

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

include("plotter.jl")

function kikkou4(nx=6, ny=5; r=1, p=0.15r, width=600, height=400)
    function unit(x, y)
        sqrt32 = sqrt(3)/2
        ax, ay = 0, 0
        bx, by = sqrt32*r, r/2
        cx, cy = sqrt32*r, -r/2
        dx, dy = 2/3*sqrt32*r - p/2, sqrt32*p
        ex, ey = 2/3*sqrt32*r + p, 0
        fx, fy = 2/3*sqrt32*r - p/2, -sqrt32*p
        x0, y0 = 2/3*sqrt32*r, 0
        xs = [ax, bx, cx, ax, dx, bx, ex, cx, fx, ax, NaN, dx, 2x0-dx, NaN, ex, 2x0-ex, NaN, fx, 2x0-fx];
        ys = [ay, by, cy, ay, dy, by, ey, cy, fy, ay, NaN, dy, -dy, NaN, ey, -ey, NaN, fy, -fy];
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:length(xs)
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2[1:10], y .+ ys2[1:10], lwd=2, col=:cornsilk)
            plotline(x .+ xs2[12:end], y .+ ys2[12:end], lwd=3, col=:snow)
        end
    end
    sqrt3 = sqrt(3)
    x1, y1, x2, y2 = r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:ny
        for x = 1:nx + 1
          unit(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * 1.5r)
        end
    end
    plotend()
end

kikkou4(4, 4, width=780, height=550)
savefig("fig4.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で亀甲ベースの格子模様を描く(3)

2021年08月12日 | ブログラミング

Julia で亀甲ベースの格子模様を描く(3)

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

include("plotter.jl")

function kikkou3(nx=6, ny=5; r=1, p=r/6, width=600, height=400)
    function unit(x, y)
        sqrt32 = sqrt(3)/2
        ax, ay = 0, 0
        bx, by = sqrt32*r, r/2
        cx, cy = sqrt32*r, -r/2
        dx, dy = 2/3*sqrt32*r - p/2, sqrt32*p
        ex, ey = 2/3*sqrt32*r + p, 0
        fx, fy = 2/3*sqrt32*r - p/2, -sqrt32*p
        x0, y0 = 2/3*sqrt32*r, 0
        xs = [ax, bx, cx, ax, dx, bx, ex, cx, fx, ax, NaN, dx, cx, NaN, ex, ax, NaN, fx, bx];
        ys = [ay, by, cy, ay, dy, by, ey, cy, fy, ay, NaN, dy, cy, NaN, ey, ay, NaN, fy, by];
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:length(xs)
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2, lwd=2, col=:cornsilk)
        end
    end
    sqrt3 = sqrt(3)
    p = 0.15r
    x1, y1, x2, y2 = r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:ny
        for x = 1:nx + 1
          unit(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * 1.5r)
        end
    end
    plotend()
end

kikkou3(4, 4, width=780, height=550)
savefig("fig3.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で亀甲ベースの格子模様を描く(2)

2021年08月12日 | ブログラミング

Julia で亀甲ベースの格子模様を描く(2)

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

include("plotter.jl")

function kikkou2(nx=6, ny=5; r=1, p=r/6, width=600, height=400)
    function unit(x, y)
        xs = [0, r, r, NaN,      r/3, r, NaN,          r/3, r, NaN,   2r/3, 2r/3] .* sqrt3/2;
        ys = [0, r/2, -r/2, NaN, -r/6, r/6, NaN, r/6, -r/6, NaN, r/3, -r/3];
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:length(xs)
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2, lwd=5, col=:cornsilk)
        end
    end
    sqrt3 = sqrt(3)
    x1, y1, x2, y2 = r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:ny
        for x = 1:nx + 1
          unit(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * 1.5r)
        end
    end
    plotend()
end

kikkou2(4, 4, width=780, height=550)
savefig("fig2.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で亀甲ベースの格子模様を描く(1)

2021年08月12日 | ブログラミング

Julia で亀甲ベースの格子模様を描く(1)

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

include("plotter.jl")

function kikkou1(nx=6, ny=5; r=1, p=r/6, width=600, height=400)
    function unit(x, y)
        xs = [0, r, r, NaN,      p, r, NaN,          p, r, NaN,           r-p, r-p] .* sqrt3/2;
        ys = [0, r/2, -r/2, NaN, p/2, (2p-r)/2, NaN, -p/2, (r-2p)/2, NaN, (r-p)/2, (p-r)/2];
        xs2 = copy(xs)
        ys2 = copy(ys)
        for θ = 0:60:300 # degree
            for i = 1:length(xs)
                xs2[i], ys2[i] = xs[i] * cosd(θ) + ys[i] * sind(θ), xs[i] * sind(θ) - ys[i] * cosd(θ)
            end
            plotline(x .+ xs2, y .+ ys2, lwd=5, col=:cornsilk)
        end
    end
    sqrt3 = sqrt(3)
    x1, y1, x2, y2 = r * sqrt3, r, (nx + 1.5) * r * sqrt3, (1.5ny + 0.5)r
    println("(width, height) = ($(x2 - x1), $(y2 - y1))")
    plotbegin(w=width, h=height)
    plotlimit(x1, y1, x2, y2)
    plotbox(x1, y1, x2, y2, col=:gray, fcol=:gray)
    for y = 1:ny
        for x = 1:nx + 1
          unit(x * sqrt3 * r + (y % 2) * sqrt3 * r / 2, y * 1.5r)
        end
    end
    plotend()
end

kikkou1(4, 4, width=780, height=550)
savefig("fig1.png")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で青海波を描く

2021年08月10日 | ブログラミング

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")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

2021年08月10日 | ブログラミング

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

# ファイル名 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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

2021年08月10日 | ブログラミング

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

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")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R の gmp, Rmpfr が動かない

2021年08月10日 | ブログラミング

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

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

で,同じように,

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

で,動くようになった。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

数学検定 過去問題 1級(大学程度・一般) 2次:数理技能検定 問題4. 統計学検定 statistical test

2021年08月01日 | ブログラミング

二次は数理技能検定というだけあって,証明問題とか出てくる。

まあ,計算問題もそうであるが,答えがあることがわかっていて,その答えをあなたは求めることが出来ますか?ということである。

答えがあるかどうかわからないものに一生掛けて取り組むというような観点からみると,くだらないと言えばくだらないんだろうなあ。

 

数学検定 過去問題 1級(大学程度・一般) 2次:数理技能検定 問題4. 統計学検定
https://www.su-gaku.net/suken/support/past_questions/

1級(大学程度・一般)

〔1級〕2次:数理技能検定

問題4. A 大学の男子学生の身長について,男子学生全体の平均は 170 cm,A 大学 B 学部の男子学生全体の標準偏差は 5.0 cm であることがわかっています。
 このとき,A 大学 B 学部の男子学生の平均身長 m(cm) について,次の方法で検定を行います。

・ 学部の男子学生 100 人を無作為抽出し,身長を測定する。なお,この 100 は十分大きな値として扱ってよいものとする。このときの B 学部の男子学生 100 人の平均身長を X(cm) とする。
・ 帰無仮説 H0 を m = 170,対立仮説 H1 を m ≠ 170 とする。
・ 大学B学部の男子学生全体の身長は,正規分布に従うと仮定する。

実際に無作為抽出された 100 人の身長を測定したところ,X = 170.9 でした。これについて,1-2-6ページの正規分布表の値を用いて,次の問いに答えなさい。

(1) 有意水準を 0.05 とするとき,帰無仮説 H0 を棄却することができるかどうかを確かめなさい。
(2) 帰無仮説 H0 を棄却するためには,有意水準をいくら以上の値にする必要があるかを求めなさい。答えは小数第 4 位を切り上げて小数第 3 位まで答えなさい。

解答

平均値は A 大学の男子学生全体について,標準偏差は B 学部の男子学生全体についてわかっているとの状況設定は,回答者を惑わすのが目的か?と思われるような不自然な状況設定である。

要するに,必要なのは「母集団の平均値と標準偏差が 170 cm,5.0 cm」としてよいというだけのこと。

サンプルサイズは 100 で,標本平均は 170.9 cm であったと。

平均値の差の検定を行うが,サンプルサイズは十分大きいとして,標準正規分布を使って両側検定をやれと。

標本平均と母平均の差の標準誤差が 5.0 / sqrt(100) であることが要点。

Z = (170.9 - 170) / (5.0 / sqrt(100))
  = 1.8

標準正規分布表を見るまでもなく, Z < 1.96 ゆえ,帰無仮説は棄却できない。

「帰無仮説 H0 を棄却するためには,有意水準をいくら以上の値にする必要があるか」というのも,トンデモ設問だ。

まあしかたない。

2 * P(Z > 1.8) = 2 * (0.5 - P(u < 1.8)) = 2 *(0.5 - 0.46407) = 0.07186

なので,0.072

有意水準 0.072 で帰無仮説を棄却できます(有意です)」なんて言うか?

このことを統計学では,「有意確率 (p value) は 0.072です」って言うんだョ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村