Juliaで回転体

関数をぐるぐる回してできる図形を回転体と呼びます。Juliaでこの回転体を描く方法について述べます。

環境

• Julia 1.9.3

コード

以下のコードを実行すると描くことができます。ここでは、z軸を中心に回転させています。

using Plots
function main()
    N = 100
    θ = range(0,2π,length=N)
    zs = range(0,1,length=N)
    f(z) = exp(-(z-0.2)^2/(0.4^2))+exp(-(z-1)^2/(0.2^2))*cos(z)
    z = []
    x = []
    y = []
    for zi in zs
        for θj in θ
            push!(z,zi)
            push!(x,f(zi)*cos(θj))
            push!(y,f(zi)*sin(θj))
        end
    end
    plot(x, y, z, label=nothing)
    x0,x1,y0,y1,z0,z1 = 0,0,0,0,0,1
    plot!([x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    x0,x1,y0,y1,z0,z1 = 0,f(1),0,0,1,1
    plot!([x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    x0,x1,y0,y1,z0,z1 = 0,f(0),0,0,0,0
    plot!([x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    plot!(f.(zs),zeros(N),zs,color = "red",lw = 1,label=nothing)

    savefig("revolution_z.png")
end
main()

x軸について回転させたければ、

using Plots
function main3()
    N = 100
    θ = range(0,2π,length=N)
    xs = range(0,1,length=N)
    f(x) = exp(-(x-0.2)^2/(0.4^2))+exp(-(x-1)^2/(0.1^2))*cos(x)^2
    z = []
    x = []
    y = []
    for xi in xs
        for θj in θ
            push!(x,xi)
            push!(z,f(xi)*cos(θj))
            push!(y,f(xi)*sin(θj))
        end
    end
    p = plot(dpi = 500)
    plot!(p,x, y, z, label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,0,0,0,0,1
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,f(1),0,0,1,1
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,f(0),0,0,0,0
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    plot!(p,xs,zeros(N),f.(xs),color = "red",lw = 1,label=nothing)

    savefig(p,"revolution_x.png")
end
main3()

とします。ここでは、pngファイルの解像度を指定してみました。

出力結果

出力結果は

となります。

x軸回転の方は

です。なお、Nを減らすと、

のように回転も減らせます。ただ、上のコードだと赤線もカクカクしてしまうので、赤線を滑らかに描くなら

function main3()
    N = 20
    N2 = 100
    θ = range(0,2π,length=N2)
    xs = range(0,1,length=N)
    xs2 = range(0,1,length=N2)
    f(x) = exp(-(x-0.2)^2/(0.4^2))+exp(-(x-1)^2/(0.1^2))*cos(x)^2
    z = []
    x = []
    y = []
    for xi in xs
        for θj in θ
            push!(x,xi)
            push!(z,f(xi)*cos(θj))
            push!(y,f(xi)*sin(θj))
        end
    end
    p = plot(dpi = 500)
    plot!(p,x, y, z, label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,0,0,0,0,1
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,f(1),0,0,1,1
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    z0,z1,y0,y1,x0,x1 = 0,f(0),0,0,0,0
    plot!(p,[x0,x1],[y0,y1],[z0,z1],color = "red",lw = 0.5,label=nothing)
    plot!(p,xs2,zeros(N2),f.(xs2),color = "red",lw = 1,label=nothing)

    savefig(p,"revolution_2x$N.png")
end
main3()

のようにすると、

が得られます。$\theta$方向の分割と赤線に関しては分割数としてN2 = 100を使っています。