算額（その1514）

算額（その1514）

七十三 岩手県一関市川崎町門崎宮畑 伊吹神社 明治17年(1884)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03080
https://w.atwiki.jp/sangaku/pages/320.html

キーワード：球8個
#Julia, #SymPy, #算額, #和算

大球の中に小球 5 個，中球 2 個を容れる。小球は互いに外接しあい，すべてが中球とも外接し，大球に内接している。中球は大球に内接している（中球同士は外接していない）。小球の直径が与えられたとき，大球の直径はいかほどか。

大球の半径と中心座標を r1, (0, 0, 0)
中球の半径と中心座標を r2, (0, 0, z)
小球の半径と中心座標を r3, (l, 0, 0); l = r3/sin(π/10)
とおき，以下の方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1, r2, r3, z
l = r3/sind(Sym(36))
eq1 = z^2 + l^2 - (r2 + r3)^2
eq2 = z + r2 - r1
eq3 = l + r3 - r1;
res = solve([eq1, eq2, eq3], (r1, r2, z))[1]

    (r3 + 2*sqrt(2)*r3/sqrt(5 - sqrt(5)), r3*(-15*sqrt(2) - 10*sqrt(5 - sqrt(5)) + 2*sqrt(25 - 5*sqrt(5)) + 5*sqrt(10))/(5*(-3*sqrt(5 - sqrt(5)) - 3*sqrt(2) + sqrt(10) + sqrt(5)*sqrt(5 - sqrt(5)))), 2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30))

# r1
res[1] |> factor |> println

    r3*(sqrt(5 - sqrt(5)) + 2*sqrt(2))/sqrt(5 - sqrt(5))  

大球の半径 r1 は，小球の半径 r3 の (sqrt(5 - √5) + 2√2)/sqrt(5 - √5) = 2.70130161670408 倍である。

術では単に「3 倍」とあるが，理解不能である。

# r2
@syms d
apart(res[2]/r3, d) |> factor |> println

    -(-30 - 10*sqrt(5) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 3*sqrt(10)*sqrt(5 - sqrt(5)))/20

# z
res[3] |> simplify |> println

    2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30)

function draw(r3, more=false)
    pyplot(size=(500, 250), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, z) = (r3 + 2*sqrt(2)*r3/sqrt(5 - sqrt(5)), r3*(-15*sqrt(2) - 10*sqrt(5 - sqrt(5)) + 2*sqrt(25 - 5*sqrt(5)) + 5*sqrt(10))/(5*(-3*sqrt(5 - sqrt(5)) - 3*sqrt(2) + sqrt(10) + sqrt(5)*sqrt(5 - sqrt(5)))), 2*r3*(-5*sqrt(5) - sqrt(50 - 10*sqrt(5)) + 5*sqrt(10 - 2*sqrt(5)) + 15)/(-10*sqrt(5) - sqrt(10)*sqrt(5 - sqrt(5)) + 5*sqrt(2)*sqrt(5 - sqrt(5)) + 30))
    l = r3/sind(Sym(36))
    p1 = plot()
    circle(0, 0, r1, :green)
    circle(0, 0, r2, :blue)
    rotate(l, 0, r3, angle=72)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
    end
    p2 = plot()
    circle(0, 0, r1, :green)
    circle22(0, z, r2, :blue)
    circle(l, 0, r3)
    circle(l*cosd(72), 0, r3)
    circle(l*cosd(144), 0, r3)
    plot(p1, p2)
end;

draw(1, true)