Juliaで点群からポリゴンを作る：フェルミ面を描く

ある三次元の配列があるとします。つまり、x,y,zを指定すると一つの値f(x,y,z)があるような状況です。この時、fがある値となるような点の集合を描きたいとします。
これは固体物理学では「フェルミ面を描く」ことに対応しています。Juliaの色々なライブラリを使うことでできることがわかりましたので、記録として残しておきます。

# バージョン
Julia 1.8.1

# インストール
必要なパッケージは以下の通りです。addでインストールしてください。
- Interpolations
- Meshing
- FileIO
- MeshIO
- GeometryBasics

# サンプルデータの生成
サンプルデータとして、
```math
E(k_x,k_y,k_z) = -(\cos(k_x) + \cos(k_y) + \cos(k_z))
```
とします。この関数がゼロとなるような点の集合を集めると、面になります。固体物理学ではこれはフェルミ面と呼ばれるものです。

この関数は解析的な形で求まっていますが、一般には数値的な値しかないかもしれません。そこで、この関数から適当に生成したサンプルデータを作成します。

```julia
using Interpolations
using Meshing
using FileIO #MeshIOもインストールしないといけない
using GeometryBasics
using LinearAlgebra

function make_data(nkx,nky,nkz)
    kxs = range(-pi,pi,length=nkx)
    kys = range(-pi,pi,length=nky)
    kzs = range(-pi,pi,length=nkz)
    energy= zeros(Float64,nkx,nky,nkz)
    for iz=1:nkz
        kz = kzs[iz]
        for iy=1:nky
            ky = kys[iy]
            for ix=1:nkx
                kx = kxs[ix]
                energy[ix,iy,iz] = -(cos(kx) + cos(ky) + cos(kz))
            end
        end
    end
    return energy
end
```

# データの補間
次に、メッシュで切ったデータ点以外の値を計算できるようにするために、スプライン補間のパッケージInterpolationsを用います。
```julia
function energy_dispersion(energymesh)
    nka,nkb,nkc = size(energymesh)
    kas = range(0,2pi,length=nka)
    kbs = range(0,2pi,length=nkb)
    kcs = range(0,2pi,length=nkc)
    interp_cubic = cubic_spline_interpolation((kas, kbs,kcs),  energymesh)

    function evaluate(ks)
        #println(typeof(ks))
        vecka = [ks[1],ks[2],ks[3]]
        for i=1:3
            while vecka[i] < 0
                vecka[i] += 2pi
            end
            while vecka[i] > 2pi
                vecka[i] -= 2pi
            end
        end
        #println(vecka)
        return interp_cubic(vecka...)
    end

    return ks -> evaluate(ks),interp_cubic 

end
```
この関数は、データのメッシュをもらうと、関数```evaluate(ks)```と補間用データを返してくれます。
ですので、
```julia
    nkx = 21
    nky = 21
    nkz = 21
    energy_mesh = make_data(nkx,nky,nkz)
    evaluate,interp_cubic = energy_dispersion(energy_mesh)

```
のようにすると、データのメッシュを21x21x21で作成した後に、補間用の関数```evaluate```を作成したことになります。


# ポリゴンデータの作成
最後に、ポリゴンを作成します。これは、
```julia
function main()
    nkx = 21
    nky = 21
    nkz = 21
    energy_mesh = make_data(nkx,nky,nkz)
    evaluate,interp_cubic = energy_dispersion(energy_mesh)

    μ = 0.0
    println(evaluate([0.2,0.1,0.3]))

    km = pi
    cosmesh = Mesh(evaluate, Rect(Vec(-km,-km,-km),Vec(2km,2km,2km)),
                        MarchingCubes(iso=-μ), samples=(50,50,50))
    save("cos.ply", cosmesh)
    a = zeros(Float64,3)
    b = zero(a)
    c = zero(a)
    g = zero(a)

    fp = open("fermi.txt","w")

    for tri in cosmesh
        kcount = 0
        for r in tri
            kcount += 1
            if kcount == 1
                a[:] = r[:]
            elseif kcount == 2
                b[:] = r[:]
            else
                c[:] = r[:]
            end
        end
        g[:] = (a+b+c)/3 #三角形の中心座標

        a[:] = a-c
        b[:] = b-c
        if norm(a-b) < 1e-10
            continue
        end
        s = sqrt(norm(a)^2*norm(b)^2 - dot(a,b)^2)/2 #三角形の面積

        println(fp,g[1],"\t",g[2],"\t",g[3],"\t",s)
    end
    
    close(fp)


end

main()
```
で可能です。
ここで、
```julia
    km = pi
    cosmesh = Mesh(evaluate, Rect(Vec(-km,-km,-km),Vec(2km,2km,2km)),
                        MarchingCubes(iso=-μ), samples=(50,50,50))
```
がポリゴンデータを作成しているところです。```Rect(Vec(-km,-km,-km),Vec(2km,2km,2km))```は、始点が```(-km,-km,-km)```で、そこからベクトル```(2km,2km,2km)```を足した点を対角線の点とした時の直方体領域を指定しています。
このデータは
```julia
save("cos.ply", cosmesh)
```
でply形式で保存ができます。ply形式は、[MeshLab](http://www.pointcloud.jp/blog_n17/)などで読み込むことができます。
例えば、
![スクリーンショット 2022-12-27 21.38.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246113/6a5ea573-12ad-cce6-080b-48107d5e3fd7.png)
こんな感じです。
また、Makie.jlも対応していますので、綺麗なプロットを描くこともできます。

サンプルコードではデータとして書き出していますから、gnuplotなどでもプロット可能で、この場合には、

![スクリーンショット 2022-12-27 21.39.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246113/5c219580-101e-a4cb-dffb-8296c8e8f1b3.png)

みたいにプロットできます。

コードではポリゴンを構成する三角形の面積を求めて、それもファイルに書き出しています。これは、表面積分を行うときの面積素として使うことができますから、表面積の計算や表面での何らかの量の計算が可能となります。