インフレーションルールに従ったPenroseタイルの生成

使用するパッケージ;

Package.ji

using Plots, LinearAlgebra, DataFrames

Penrose tilingのインフレーションに関しては
https://scipython.com/blog/penrose-tiling-1/
などが詳しい。
要点としては、2種類の三角形に関して特定の条件で辺を抜き出し、その線分を黄金比を用いて内分するといったものである。

まず、初期条件となる三角形の集合を生成する

initial-triangle.jl

function get_init()
    triangle = [] #格子点のデータを初期化
    for i in 1:10
        θi =-2π/5 + (2i-1)*π/10
        θi1 =-2π/5 + (2i+1)*π/10
        B = [cos(θi), sin(θi)]
        C = [cos(θi1), sin(θi1)]
        if i%2 == 0
            push!(triangle, [0, [0.0,0.0], B, C]) #[col::Int, A::Vect, B::Vect, C::Vect]
        else
            push!(triangle, [0, [0.0,0.0], C, B])
        end
    end
    return triangle
end

これにより次のような三角形の集合が形成される

この図は次の関数を用いて出力可能である。

Plot-triangle.jl

function plot_triangle(triangle,name)
    plot()
    for T in triangle
        col = T[1]
        list = [T[2], T[3], T[4]]
        sort!(list, by = x -> x[1])
        
        A = list[1]
        B = list[3]
        C = list[2]
        
        y1 = range(A[2], stop=C[2], length=10)
        x1 = range(A[1], stop=C[1], length=10)
        y2 = range(C[2], stop=B[2], length=10)
        x2 = range(C[1], stop=B[1], length=10)

        # 三角形の頂点をプロット
        x = [A[1], B[1], C[1], A[1]]
        y = [A[2], B[2], C[2], A[2]]
        
        if col == 0
            plot!(x, y, seriestype = :shape, fillcolor = :red, label = :none)
        else
            plot!(x, y, seriestype = :shape, fillcolor = :blue, label = :none)
        end
    end
    savefig(name)
end

次にこの初期条件を用いて、インフレーションルールに基づいて三角形の分割を行う

Spilit-triangle.jl

function split_triangle(triangle)
    new_triangles = [] #新たな格子点データを格納するリストを初期化
    goldenRatio = (1 + sqrt(5))/2
    for T in triangle
        col = T[1]
        A = T[2]
        B = T[3]
        C = T[4]
        if col == 0 
            P = A + (B - A) / goldenRatio #三角形を内分
            push!(new_triangles, [0, C, P, B])
            push!(new_triangles, [1, P, C, A])
        else 
            Q = B + (A - B) / goldenRatio #三角形を内分
            R = B + (C - B) / goldenRatio #三角形を内分
            push!(new_triangles, [1, Q, R, B])
            push!(new_triangles, [0, R, Q, A])
            push!(new_triangles, [1, R, C, A])
        end
    end
    return new_triangles
end

一回の分割操作によって次のような図形が生成される

この操作を5回繰り返した図形を例として示す

以上により、インフレーションルールに従って格子点のデータを生成することに成功した。
以下では、このデータを物性理論などの様々な解析に使うため、最近接サイトなどのデータをわかりやすい形式で格納し、出力する方法を示す。
以下はhttps://qiita.com/cometscome_phys/items/659a569497eef3a85d86 を参考にして作成したものである。

Make-data.jl

function is_inside_pentagon(x, y, R) #点(x,y)が正五角形領域に含まれるかを判定(無くても良い)
    θ = [π/2, π/2 + 2π/5, π/2 + 4π/5, π/2 + 6π/5, π/2 + 8π/5]
    vertices = [(R*cos(θ[i]), R*sin(θ[i])) for i in 1:5]

    function sign(p1, p2, p3)
        return (p1[1] - p3[1]) * (p2[2] - p3[2]) - (p2[1] - p3[1]) * (p1[2] - p3[2])
    end

    function point_in_triangle(pt, v1, v2, v3)
        d1 = sign(pt, v1, v2)
        d2 = sign(pt, v2, v3)
        d3 = sign(pt, v3, v1)

        has_neg = (d1 < 0) || (d2 < 0) || (d3 < 0)
        has_pos = (d1 > 0) || (d2 > 0) || (d3 > 0)

        return !(has_neg && has_pos)
    end

    for i in 1:5
        if point_in_triangle((x, y), vertices[i], vertices[mod1(i+1, 5)], (0, 0))
            return true
        end
    end

    return false
end

function make_data(triangle) #格子点の原点からの距離のデータを生成し、格納
    #作成した格子点を境界条件に従ってカットするときなどに使用できる
    Ax,Ay = [triangle[i][2][1] for i in 1:length(triangle)], [triangle[i][2][2] for i in 1:length(triangle)]
    Bx,By = [triangle[i][3][1] for i in 1:length(triangle)], [triangle[i][3][2] for i in 1:length(triangle)]
    Cx,Cy = [triangle[i][4][1] for i in 1:length(triangle)], [triangle[i][4][2] for i in 1:length(triangle)]

    dfA = DataFrame([1:length(triangle),Ax,Ay,zeros(Int64,length(triangle))], [:id, :Ax, :Ay,:name]) #A; name = 0
    dfB = DataFrame([1:length(triangle),Bx,By,zeros(Int64,length(triangle)).+1 ], [:id, :Bx, :By,:name]) #B; name = 1
    dfC = DataFrame([1:length(triangle),Cx,Cy,zeros(Int64,length(triangle)).+2 ], [:id, :Cx, :Cy,:name]) #C; name = 2

    df2 = rename!(dfA[:,[:id,:Ax,:Ay,:name]], :Ax => :x, :Ay => :y)
    df2 = append!(df2, rename!(dfB[:,[:id,:Bx,:By,:name]], :Bx => :x, :By => :y))
    df2 = append!(df2,rename!(dfC[:,[:id,:Cx,:Cy,:name]], :Cx => :x, :Cy => :y))
    df2.Abs2 = df2.x.^2 + df2.y.^2

    return sort!(df2, :Abs2)
end


function make_points_df(df::DataFrame)

    ids = [] #格子点のid
    names = [] #格子点が三角形のどの頂点に対応するのか
    points = [] #格子点のベクトルデータ
    num_hop = [] #ある格子点が接続する格子点の数(coodination number)
    hop_ind = [] #ある格子点が接続する格子点のid (length(hop_ind[i]) == num_hop[i])
    
    df.Abs2 = round.(df.Abs2, digits=5)
    df.x = round.(df.x, digits=5)
    df.y = round.(df.y, digits=5)
    for abs in unique(df.Abs2)
        df2 = df[df.Abs2 .== abs, :]
        
        for (j, id) in enumerate(df2.id)
            if [df2.x[j], df2.y[j]] ∉ points && is_inside_pentagon(df2.x[j], df2.y[j], 1.05)
                push!(points, [df2.x[j], df2.y[j]])
                push!(num_hop, 0) #num_hopを初期化 (0でfill)
                push!(hop_ind, [-1])   #hop_indを初期化 (ホッピングなし → [-1])
            end
        end

        push!(ids, df2.id)
        push!(names, df2.name)
    end

    df_points = DataFrame([1:length(points), points,hop_ind,num_hop], [:id, :point,:hop_ind,:num_hop])
    
    df_points.id = 1:size(df_points)[1]

    for (i,p) in enumerate(df_points.point)
        df_same = df[df.x .== p[1], :] #同じ三角形に属する格子点のみを抜き出す
        df_same = df_same[df_same.y .== p[2], :]

        for (j, id) in enumerate(df_same.id)
            df_triangle = df[df.id .== id, :]
            if df_same[j,:].name[1] == 0
                df_B = df_triangle[df_triangle.name.==1,:]
                df_C = df_triangle[df_triangle.name.==2,:]

                B = [df_B.x[1], df_B.y[1]]
                C = [df_C.x[1], df_C.y[1]]

                df_connected_B = df_points[[ p == B for p in df_points.point ],:]
                df_connected_C = df_points[[ p == C for p in df_points.point ],:]

                df_connected = append!(df_connected_B, df_connected_C) #接続される格子点のデータ
                
                for id_connec in df_connected.id 
                    if df_points[i,:num_hop] == 0
                        df_points[i,:hop_ind] = [id_connec]
                        df_points[i,:num_hop] += 1
                    elseif id_connec ∉ df_points[i,:hop_ind]
                        push!(df_points[i,:hop_ind], id_connec)
                        df_points[i,:num_hop] += 1
                    end
                end
            else 
                df_A = df_triangle[df_triangle.name.==0,:]
                
                A = [df_A.x[1], df_A.y[1]]
                df_connected = df_points[[ p == A for p in df_points.point ],:]

                for id_connec in df_connected.id 
                    if df_points[i,:num_hop] == 0
                        df_points[i,:hop_ind] = [id_connec]
                        df_points[i,:num_hop] += 1
                    elseif id_connec ∉ df_points[i,:hop_ind]
                        push!(df_points[i,:hop_ind], id_connec)
                        df_points[i,:num_hop] += 1
                    end
                end
            end
            
        end
    end

    return df_points
end


function save_data(df_points::DataFrame) #作成した格子点データを出力
    file_dot = open("dotxy_"*string(size(df_points)[1])*".in","w")
    file_hop = open("hopindices_"*string(size(df_points)[1])*".in","w")
    file_hopnum = open("hopnum_"*string(size(df_points)[1])*".in","w")
    for i in 1:size(df_points)[1]
        println(file_dot,df_points[i,:point][1]," ",df_points[i,:point][2])
        println(file_hopnum,df_points[i,:num_hop])
        for hop in 1:df_points[i,:num_hop]
            print(file_hop,df_points[i,:hop_ind][hop]," ")
        end
        println(file_hop)
    end
    close(file_dot)
    close(file_hop)
    close(file_hopnum)
end

実行例;

test.jl

triangle = get_init()
for i in 1:5
    triangle = split_triangle(triangle)
end

df_m = make_data(triangle)
df_points = make_points_df(df_m)
save_data(df_points)


plt = plot(legend = false)
for i in 1:size(df_points)[1]
    for hop in 1:df_points[i,:num_hop]
        j = df_points[i,:hop_ind][hop]

        plot!([df_points[i,:point][1],df_points[j,:point][1]], [df_points[i,:point][2],df_points[j,:point][2]], color="red", lw=0.5,label="")
    end
end
plt