算額（その1519）

算額（その1519）

百 大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

キーワード：円14個，外円
#Julia, #SymPy, #算額, #和算

大円の中に3個の中円を容れたもの 2 個と小円 6 個が図のように互いに接している。小円の直径が 1 寸のとき，甲円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (x1, 2r2/√3)
小円の半径と中心座標を r3, (0, r3), (r3, y3)
とおき，以下の連立方程式を解く。
甲円の半径 r1 は，（筆算により）中円の半径 r2 の (2/√3 + 1) 倍である。

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

using SymPy
@syms r1, x1, r2, y2, r3, y3
r1 = r2*(2/√Sym(3) + 1)
y2 = r1 - r2
eq1 = x1^2 + r3^2 - (r1 + r3)^2
eq2 = (x1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = r3^2 + (y3 - r3)^2 - 4r3^2
res = solve([eq1, eq2, eq3], (x1, r2, y3))[4]

    (r3*(sqrt(3) + 2), -2*sqrt(3)*r3 + 2*sqrt(3)*r3*sqrt(15*sqrt(3) + 26)/(4*sqrt(3) + 7) + 3*r3, r3*(1 + sqrt(3)))

中円の半径 r2 は，小円の半径 r3 の (3 + 3√2 - 2√3 - √6) 倍である。
小円の直径が 1 寸のとき，中円の直径は (3 + 3√2 - 2√3 - √6) = 1.3290493291983534 寸である。

「答」は「一寸三分二厘八毛余」と末尾に若干の誤差がある。

res[2] |> factor |> sympy.sqrtdenest |> simplify |> println

    -2*sqrt(3)*r3 - sqrt(6)*r3 + 3*r3 + 3*sqrt(2)*r3

3 + 3√2 - 2√3 - √6

    1.3290493291983534

function draw(r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (x1, r2, y3) = (r3*(sqrt(3) + 2), -2*sqrt(3)*r3 + 2*sqrt(3)*r3*sqrt(15*sqrt(3) + 26)/(4*sqrt(3) + 7) + 3*r3, r3*(1 + sqrt(3)))
    r1 = r2*(2/√3 + 1)
    #(r1, x1, r2, y3) = (-16*sqrt(3)*r3*sqrt(15*sqrt(3) + 26) - 8*sqrt(3)*r3*sqrt(45*sqrt(3) + 78) - r3 + 14*r3*sqrt(45*sqrt(3) + 78) + 28*r3*sqrt(15*sqrt(3) + 26), r3*(sqrt(3) + 2), -2*sqrt(3)*r3 + 2*sqrt(3)*r3*sqrt(15*sqrt(3) + 26)/(4*sqrt(3) + 7) + 3*r3, r3*(1 + sqrt(3)))
    println("x1 = $x1, r1 = $r1, 2r2 = $(2r2)")
    plot()
    circle(x1, 0, r1, :green)
    circle(-x1, 0, r1, :green)
    circle2(x1, 2r2/√3, r2, :blue)
    circle2(x1 + 2r2/√3*cosd(30), -2r2/√3*sind(30), r2, :blue)
    circle2(-x1 + 2r2/√3*cosd(30), -2r2/√3*sind(30), r2, :blue)
    circle22(0, r3, r3)
    circle4(r3, y3, r3)
    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)
        point(x1, 0, "大円:r1,(x1,r1)", :green, :center, delta=-delta)
        point(x1, 2r2/√3, "中円:r2,(x1,2r2/√3)", :blue, :center, delta=-delta)
        point(0, r3, "小円:r3\n(0,r3)", :red, :center, delta=-delta)
        point(r3, y3, "小円:r3\n(r3,y3)", :red, :center, delta=-delta)
    end
end;

draw(1/2, true)