算額（その1485）

算額（その1485）

四十四 岩手県一関市真滝 熊野白山滝神社 文久元年(1861)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，二等辺三角形，弦2本
#Julia, #SymPy, #算額, #和算

問題文，答，術ともに不明であるが，以下のようなものでもあろう。

外円の中に二等辺三角形と 2 本の弦を容れ，隙間に大円 1 個，中円 1 個，小円 2 個を容れる。外円の直径が与えられたとき，大円，中円，小円の直径を求める術を述べよ。

注：これだけの条件では大円の大きさは不定である（中円・小円とは無関係に，どのような大きさにもなれる）。これを解消するための制約条件はいくつか考えられるが，ここでは，外円の円周上にある ５ 点が，外円に内接する正五角形の頂点であるということにする。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
中円の半径と中心座標を r2, (0, R*sind(18) + r2)
小円の半径と中心座標を r3, (R*sind(18) + r3)
とおき，以下の連立方程式を解く。
なお，∠EAD = 54°，∠BAC = 18°，∠ADE = 36° である。

include("julia-source.txt");

using SymPy

@syms R, r1, r2, r3, x3;
# eq1 = dist2(R*sind(Sym(36)), -R*cosd(Sym(36)), 0, R, 0, r1 - R, r1)
# eq2 = dist2(R*cosd(Sym(18)), R*sind(Sym(18)), 0, R, 0, R*sind(Sym(18)) + r2, r2)
# eq3 = dist2(R*cosd(Sym(18)), R*sind(Sym(18)), 0, R, x3, R*sind(Sym(18)) + r3, r3)
# eq4 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2

eq11 = r1/(2R - r1) - sind(Sym(18));
eq12 = r2/(R - (r2 + R*sind(Sym(18)))) - sind(Sym(54)) |> simplify;
eq13 = r2*(R*cosd(Sym(18)) - x3) - r3*(R*cosd(Sym(18)));
eq14 = x3 - 2sqrt(r2*r3);
res = solve([eq11, eq12, eq13, eq14], (r1, r2, r3, x3))[2]

    (-4*R + 2*sqrt(5)*R, R*(-1 + sqrt(5))/4, R*(-35 - 6*sqrt(30 - 10*sqrt(5)) + 10*sqrt(6 - 2*sqrt(5)) + 19*sqrt(5))/20, R*((-5*sqrt(2) + 2*sqrt(10))*sqrt(sqrt(5) + 5) + 4*sqrt(25 - 10*sqrt(5)))/10)

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

    2*R*(-2 + sqrt(5))

# r2
res[2] |> println

    R*(-1 + sqrt(5))/4

# r3
res[3] |> sympy.sqrtdenest |> factor |> println

    R*(-15 + 7*sqrt(5))/4

# x3
res[4] |>sympy.sqrtdenest |> simplify |> println

    R*(-5*sqrt(2*sqrt(5) + 10) + 4*sqrt(25 - 10*sqrt(5)) + 2*sqrt(10*sqrt(5) + 50))/10

大円，中円，小円の半径は，外円の半径の 2(√5 - 2) 倍，(√5 - 1)/4 倍，(7√5 - 15)/4 倍である。

外円の直径が 1 のとき，大円，中円，小円の直径は 0.4721359549995796, 0.30901699437494745, 0.16311896062463216 である。

function draw(R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, r3, x3) = (-4*R + 2*sqrt(5)*R, R*(-1 + sqrt(5))/4, R*(-35 - 6*sqrt(30 - 10*sqrt(5)) + 10*sqrt(6 - 2*sqrt(5)) + 19*sqrt(5))/20, R*((-5*sqrt(2) + 2*sqrt(10))*sqrt(sqrt(5) + 5) + 4*sqrt(25 - 10*sqrt(5)))/10)
    @printf("外円の直径が %g のとき，大円，中円，小円の直径は %g, %g, %g である。\n", 2R, 2r1, 2r2, 2r3)
    r1 = 2*R*(-2 + sqrt(5))
    r2 = R*(-1 + sqrt(5))/4
    r3 = R*(-15 + 7*sqrt(5))/4
    x3 = R*(-5*sqrt(2*sqrt(5) + 10) + 4*sqrt(25 - 10*sqrt(5)) + 2*sqrt(10*sqrt(5) + 50))/10
    y = R*sind(18)
    plot([0, -R*cosd(18), R*cosd(18), 0], [R, y, y, R], color=:green, lw=0.5)
    circle(0, 0, R, :blue)
    circle(0, r1 - R, r1)
    circle(0, R*sind(18) + r2, r2, :magenta)
    circle2(x3, R*sind(18) + r3, r3, :orange)
    segment(0, R, R*sind(36), -R*cosd(36))
    segment(0, R, -R*sind(36), -R*cosd(36))
    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(0, r1 - R, "大円:r1,(0,r1-R)", :red, :center, delta=-delta)
        point(0, R - r2, "中円:r2,(0,R-r2)", :magenta, :center, delta=-delta)
        point(x3, R*sind(18) + r3, "小円:r3,(x3,R*sind(18)+r3)", :black, :center, delta=-delta)
        point(0, R, "A", :blue, :center, :bottom, delta=delta/2)
        point(0, -R, "B", :blue, :center, delta=-delta)
        point(R*sind(36), -R*cosd(36), " C", :blue, :left)
        point(R*cosd(18), R*sind(18), "  D", :blue, :left, :vcenter)
        point(0, R*sind(18), "E ", :blue, :right, :vcenter)
    end
end;

draw(1/2, true)