算額（その1522）

算額（その1522）

四十 岩手県一関市牧沢 牧沢八幡神社 明治5年(1872)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03060
https://w.atwiki.jp/sangaku/pages/308.html

キーワード：円5個，斜線2本
#Julia, #SymPy, #算額, #和算

大円の中に等円を 3 個入れ，大円に外接する外円と等円の共通接線を引く。大円，外円の直径がそれぞれ 7 寸，2 寸のとき，等円の直径はいかほどか。

山村の図は不適切である。「今有如図」の図に基づいて解を求める。

大円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r1, (x1, y1), (r1 - R, 0)
外円の半径と中心座標を r2, (R + r2, 0)
共通接線と外円の接点座標を (x0, sqrt(r2^2 - (R + r2 - x0)^2))
とおき，以下の連立方程式を解く。

SymPy の性能上，一度には解けないので，逐次解いていく。

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

using SymPy
@syms R::positive, r1::positive, x1::positive, y1::positive,
      r2::positive, x0::positive
y0 = sqrt(r2^2 - (R + r2 - x0)^2)
eq1 = (x1 + (R - r1))^2 + y1^2 - 4r1^2
eq2 = x1^2 + y1^2 - (R - r1)^2
eq3 = dist2(-R, 0, x0, y0, x1, y1, r1)
eq4 = dist2(-R, 0, x0, y0, R + r2, 0, r2);
#res = solve([eq1, eq2, eq3, eq4], (r1, x1, y1, x0))

まず，eq1, eq2, eq4 の連立方程式を解いて，x1, y1, x0 を求める。

res = solve([eq1, eq2, eq4], (x1, y1, x0))[1]

    (-(-R^2 + 2*R*r1 + r1^2)/(-R + r1), -2*sqrt(R)*r1*sqrt(R - 2*r1)/(-R + r1), R*(2*R + 3*r2)/(2*R + r2))

得られた解を eq3 に代入し，新たな方程式 eq13 とする。

eq13 = eq3(x1 => res[1], y1 => res[2], x0 => res[3]) |> simplify |> numerator;

方程式を解き，r1 の解を求める。
r1 = 2sqrt(R + r2)*sqrt(3R + 2r2) - 3(R + r2) である。

res2 = solve(eq13, r1)[2]
res2 |> println

    -3*R - 3*r2 + 2*sqrt(R + r2)*sqrt(3*R + 2*r2)

R = 7/2 寸, r2 = 2/2 寸のとき，r1 = 2sqrt(R + r2)*sqrt(3R + 2r2) - 3(R + r2) = 1.5 となる。等円の直径は 3 寸である。

res2(R => 7//2, r2 => 2//2) |> println
res2(R => 7//2, r2 => 2//2).evalf() |> println

    3/2
    1.50000000000000

function draw(R, r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = -3*R - 3*r2 + 2*sqrt(R + r2)*sqrt(3*R + 2*r2)
    (x1, y1, x0) = (-(-R^2 + 2*R*r1 + r1^2)/(-R + r1), -2*sqrt(R)*r1*sqrt(R - 2*r1)/(-R + r1), R*(2*R + 3*r2)/(2*R + r2))
    y0 = sqrt(r2^2 - (R + r2 - x0)^2)
    plot()
    circle(0, 0, R, :green)
    circle(r1 - R, 0, r1, :blue)
    circle22(x1, y1, r1, :blue)
    circle(R + r2, 0, r2)
    segment(-R, 0, x0, y0)
    segment(-R, 0, x0, -y0)
    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, R, "R", :green, :center, :bottom, delta=delta)
        point(R + r2, 0, "外円:r2\n(R+r2,0)", :red, :center, delta=-delta)
        point(x0, y0, "(x0,y0)", :black, :center, :bottom, delta=delta)
        point(r1 - R, 0, "等円:r1,(r1-R,0)", :blue, :center, :bottom, delta=3delta)
        point(x1, y1, "等円:r1,(x1,y1)", :blue, :center, :bottom, delta=delta)
    end
end;

draw(7/2, 2/2, true)