算額（その1561）
六十五 岩手県花泉町金沢 大門神社・大門観世音菩薩 明治33年(1900)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，外円，二等辺三角形
#Julia, #SymPy, #算額, #和算
甲円と交差する二等辺三角形を設け，隙間に乙円 3 個を容れる。乙円の直径が与えられたとき，甲円の直径はいかほどか。
甲円の半径と中心座標を R, (0, 0)
乙円の半径と中心座標を r, (0, -R - r), (x, y)
二等辺三角形の底辺の端点の座標を (-R - 2r, a)
とおき，以下の連立方程式を解く。
include("julia-source.txt")
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R, r, x, y, a
eq1 = r/(2R + r) - a/sqrt((2R + 2r)^2 + a^2)
eq2 = dist2(0, R, a, -R - 2r, x, y, r)
eq3 = x^2 + y^2 - (R - r)^2
eq4 = y/x - a/(2R + 2r);
一度に解けないので，eq1, eq3, eq4 から a, x, y を求める。
res = solve([eq1, eq3, eq4], (a, x, y))[2] # 2 of 4
(r*sqrt(R*(R^3 - R^2*r - R*r^2 + r^3)/(4*R^2 + 4*R*r + r^2))*(-2*R - r)/(R*(-R + r)), 2*sqrt(R*(R^3 - R^2*r - R*r^2 + r^3)/(4*R^2 + 4*R*r + r^2)), r*(R - r)/(2*R + r))
eq2 に代入し，eq12 式を求める。
eq12 = eq2(a => res[1], x => res[2], y => res[3]) |> simplify
eq12 |> println
4*R*(R^3 - R^2*r - 3*R*r^2 - r^3)
eq12 を解いて R を求める。
ans_R = solve(eq12, R)[4] # 4 of 4
ans_R |> println
r*(1 + sqrt(2))
甲円の直径は，乙円の直径の (1 + √2) 倍である。
乙円の半径が 1 のとき，全てのパラメータは以下のとおりである。
r = 1; R = 2.41421; a = 1.18921; x = 1.39324; y = 0.242641
function draw(r, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = r*(1 + √2)
(a, x, y) = (r*sqrt(R*(R^3 - R^2*r - R*r^2 + r^3)/(4R^2 + 4R*r + r^2))*(-2R - r)/(R*(-R + r)), 2sqrt(R*(R^3 - R^2*r - R*r^2 + r^3)/(4R^2 + 4R*r + r^2)), r*(R - r)/(2R + r))
@printf("r = %g; R = %g; a = %g; x = %g; y = %g\n", r, R, a, x, y)
plot([a, -a, 0, a], [-R - 2r, -R - 2r, R, -R - 2r], color=:red, lw=0.5)
circle(0, 0, R, :green)
circle2(x, y, r, :blue)
circle(0, -R - r, r, :blue)
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/2)
point(0, -R, "-R", :green, :center, :bottom, delta=delta/2)
point(x, y, "乙円:r,(x,y)", :blue, :center, delta=-delta/2)
point(0, -R - r, "乙円:r,(0,-R-r)", :blue, :center, delta=-delta/2)
end
end;
draw(1, true)
