算額（その984）
一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
外円の中に，甲円 1 個，乙円 2 個，丙円 2 個が入っている。外円，甲円，乙円の直径がそれぞれ 156 寸，52 寸，26 寸のとき，丙円の直径はいかほどか。
外円の半径と中心座標を R, (0, )
甲円の半径と中心座標を r1, (0, y1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (x3, y3)
とおいて，以下の連立方程式を解く。
注：丙円の中心は x 軸上にあるわけではない。y3 ≠ 0
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, y1::negative,
r2::positive, y2::positive,
r3::positive, x3::positive, y3::negative
eq1 = x3^2 + (y1 - y3)^2 - (r3 - r1)^2
eq2 = x3^2 + (y2 - y3)^2 - (r3 - r2)^2
eq3 = x3^2 + (y3 - r2 + R)^2 - (r2 + r3)^2
eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = (y2 - y1) - (r1 + r2);
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, y1, y2, x3, y3))[1];
println("r3: ", res[1])
println("y1: ", res[2])
println("y2: ", res[3])
println("x3: ", res[4])
println("y3: ", res[5])
r3: (-R + r2)*(r1 + r2)*(-2*R*r1 - R*r2 + r1*r2)/((R + r1)*(R*r1 + r2^2))
y1: -(-R + r1)*(-R + r1 + 2*r2)/(R + r1)
y2: -(R^2 - 3*R*r1 - 3*R*r2 + r1*r2)/(R + r1)
x3: 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-(-R + r2)*(-R + r1 + 2*r2)*(-2*R*r1 - R*r2 + r1*r2))*sqrt(r1 + r2)/((R + r1)*(R*r1 + r2^2))
y3: (-R^3*r1 + R^2*r1^2 + 3*R^2*r1*r2 + R*r1^2*r2 + R*r1*r2^2 + R*r2^3 - r1^2*r2^2 - r1*r2^3)/(R^2*r1 + R*r1^2 + R*r2^2 + r1*r2^2)
外円，甲円，乙円の直径がそれぞれ 156 寸, 52 寸, 26 寸のとき，丙円の直径は 105 寸である。
その他のパラメータは以下のとおりである。
r1 = 26; r2 = 13; R = 78
r3 = 52.5; y1 = -13; y2 = 26; x3 = 25.0998; y3 = -4.5
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r2) = (156, 52, 26) ./ 2
(r3, y1, y2, x3, y3) = ((-R + r2)*(r1 + r2)*(-2*R*r1 - R*r2 + r1*r2)/((R + r1)*(R*r1 + r2^2)), -(-R + r1)*(-R + r1 + 2*r2)/(R + r1), -(R^2 - 3*R*r1 - 3*R*r2 + r1*r2)/(R + r1), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-(-R + r2)*(-R + r1 + 2*r2)*(-2*R*r1 - R*r2 + r1*r2))*sqrt(r1 + r2)/((R + r1)*(R*r1 + r2^2)), (-R^3*r1 + R^2*r1^2 + 3*R^2*r1*r2 + R*r1^2*r2 + R*r1*r2^2 + R*r2^3 - r1^2*r2^2 - r1*r2^3)/(R^2*r1 + R*r1^2 + R*r2^2 + r1*r2^2))
@printf("外円，甲円，乙円の直径がそれぞれ %g, %g, %g のとき，丙円の直径は %g である。\n", 2R, 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; R = %g\n", r1, r2, R)
@printf("r3 = %g; y1 = %g; y2 = %g; x3 = %g; y3 = %g\n", r3, y1, y2, x3, y3)
plot()
circle(0, 0, R, :green)
circle(0, y1, r1)
circle(0, y2, r2, :blue)
circle(0, r2 - R, r2, :blue)
circle2(x3, y3, 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(0, y1, "甲円:r1\n(0,y1)", :brown, :center, delta=-delta/2)
point(0, y2, "乙円:r2\n(0,y2)", :blue, :center, delta=-delta/2)
point(0, r2 - R, "乙円:r2\n(0,r2-R)", :blue, :center, :bottom, delta=delta/2)
point(x3, y3, " 丙円:r3,(x3,y3)", :blue, :left, :vcenter)
point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
end
end;
