算額（その1510）

算額（その1510）

九十九 岩手県江刺市 雨宝堂 現在中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03007
https://w.atwiki.jp/sangaku/pages/205.html

キーワード：円4個，直角三角形，斜線1本
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を引き，大円 1 個，中円 1 個，小円 2 個を容れる。鈎，股が 3 寸，4 寸のとき，中円の直径はいかほどか。

鈎，股をそのまま「鈎」，「股」
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (股 - r2, y2)
小円の半径と中心座標を r3, (x3, r3), (x32, y32)
斜線と斜辺の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

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, x3, x32, y32, x0, y0
eq1 = dist(0, 0, 股, 鈎, x1, r1) - r1^2
eq2 = dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2
eq3 = dist(0, 0, 股, 鈎, x32, y32) - r3^2
eq4 = dist(x0, y0, 股, 0, x1, r1) - r1^2
eq6 = dist(x0, y0, 股, 0, x32, y32) - r3^2
eq7 = dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2
eq8 = x0^2 + y0^2 - 股^2
eq9 = (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0
eq11 = r3/(股 - x3) - r1/(股 - x1)
eq12 = (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2;

function H(u)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = u
    return [
        dist(0, 0, 股, 鈎, x1, r1) - r1^2,
        dist(0, 0, 股, 鈎, 股 - r2, y2) - r2^2,
        dist(0, 0, 股, 鈎, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, x1, r1) - r1^2,
        dist(x0, y0, 股, 0, x32, y32) - r3^2,
        dist(x0, y0, 股, 0, 股 - r2, y2) - r2^2,
        x0^2 + y0^2 - 股^2,
        (2股 + sqrt((股 - x0)^2 + y0^2))*r1 - 股*y0,
        r3/(股 - x3) - r1/(股 - x1),
        (x3 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
    ]
end;
(鈎, 股) = (3, 4)
iniv = BigFloat[0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
res = nls(H, ini=iniv)

    ([0.9116963119775494, 2.7350889359326485, 0.3675444679663241, 2.2649110640673515, 0.2389194451114356, 3.668516977010949, 3.0781652486756204, 2.0099746301174206, 3.2, 2.4], true)

中円の直径は res[1][3]*2 = 0.7350889359326482 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (鈎, 股) = (3, 4)
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = [0.91753, 2.72942, 0.36849, 2.26201, 0.2366, 3.6581, 3.07768, 2.00873, 3.2, 2.4]
    (r1, x1, r2, y2, r3, x3, x32, y32, x0, y0) = res[1]
    plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
    circle(x1, r1, r1)
    circle(股 - r2, y2, r2, :blue)
    circle(x3, r3, r3, :green)
    circle(x32, y32, r3, :green)
    segment(股, 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(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(股 - r2, y2, "中円:r2\n(股-r2,y2)", :blue, :center, :bottom, delta=delta)
        point(x3, r3, "小円:r3,(x3,r3)", :green, :right, :vcenter, deltax=-8delta)
        point(x32, y32, "小円:r3,(x32,y32)", :green, :right, :bottom, delta=2delta, deltax=-6delta)
        point(股, 0, " 股", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
        point(股, 鈎, "(股,鈎)", :magenta, :left, :bottom, delta=delta, deltax=-delta/2)
    end
end;

draw(2, true)