算額（その198）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（92）
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)
問4 全円の中に大円，中円，小円が入っている。大円，中円，小円の径が与えられたとき，全円の径を求めよ。
図のように，全円，大円，中円，小円の半径を R, r1, r2, r3 とし，小円の中心座標を (x3, y3) とする。連立方程式を立て，R, x3, y3 について解く。
using SymPy
@syms r1::positive, r2::positive, r3::positive, R::positive,
x3::positive, y3::negative;
eq1 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = x3^2 + (R - 2r1 - r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (R, x3, y3));
println("R = ", res[1][1])
println("x3 = ", res[1][2])
println("y3 = ", res[1][3])
R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;
function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function draw(r1, r2, r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))
@printf("R = %.5f, x3 = %.5f, y3 = %.5f\n", R, x3, y3)
plot()
circle(0, 0, R)
circle(0, R - r1, r1, :blue)
circle(x3, y3, r3, :green)
circle(-x3, y3, r3, :green)
circle(0, R - 2r1 - r2, r2, :orange)
if more == true
point(0, R - r1, " R-r1", :blue)
point(0, R - 2r1 - r2, " R-2r1-r2", :orange)
point(x3, y3, "(x3,y3)", :green)
point(0.4r1, R - 0.6r1, "大円", :blue, mark=false)
point(0.4r2, R - 2r1 - 0.6r2, "中円", :orange, mark=false)
point(x3 + 0.4r3, y3 + 0.4r3, "小円", :green, mark=false)
point(0.75R, 0.75R, "全円", :red, mark=false)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
draw(6, 2, 4, false)
R = 10.80000, x3 = 6.00000, y3 = -3.20000
