算額（その1631）
～落書き帳「○△□」～ 940. 『算法天生指南』巻之二（その12）
http://streetwasan.web.fc2.com/math20.05.03.html
キーワード：円3個，長方形
#Julia, #Julia, #SymPy, #算額, #和算, #数学
長方形の中に甲円，乙円，丙円を容れる。甲円の直径が 12 寸，乙円の直径が 9 寸のとき，丙円の直径はいかほどか。
算額（その1630）で半円だったものが円になっただけで，本質的には同じ問題である。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (x3, 2r1 - r3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, x1::positive, r2::positive,
r3::positive, x3::positive
eq1 = (x1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x3 - r2)^2 + (2r1 - r3 - r2)^2 - (r2 + r3)^2;
res = solve([eq1, eq2, eq3], (x1, r3, x3))[2]; # 2 of 2
# r3
res[2] |> println
r1^2/(4*r2)
乙円の半径 r3 は，甲円と乙円の半径 r1, r2 の関数 r1^2/(4*r2) で表される。
甲円の直径が 12 寸，乙円の直径が 9 寸のとき，丙円の直径は 4 寸である。
function draw(r1, r2, more)
直平 = 2r1
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, r3, x3) = ((2*r1^(3/2)*sqrt(r2) - 4*sqrt(r1)*r2^(3/2) + r1*r2 - 2*r2^2)/(r1 - 2*r2), r1^2/(4*r2), -sqrt(r1)*(r1 - 2*r2)/sqrt(r2) + r2)
@printf("甲円，乙円の直径が %g, %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
直長 = x1 + r1
plot([0, 直長, 直長, 0, 0], [0, 0, 直平, 直平, 0], color=:green, lw=0.5)
circle(x1, r1, r1)
circle(r2, r2, r2, :blue)
circle(x3, 2r1 - r3, r3, :magenta)
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, r2, "乙円:r2,(r2,r2)", :blue, :center, delta=-delta)
point(x3, 2r1 - r3, "丙円:r3,(x3,2r1-r3)", :magenta, :right, :vcenter, deltax=-2delta)
end
end;
draw(12/2, 9/2, true)
