算額（その9）
大きな円の中に 4 種類の大きさの異なる円 A, B, C, D があり，図のように互いに接している（すべての円が相互に接しているわけではない）。
外側の円の半径を 1 としたとき，それぞれの円の半径を求めよ。
山形県鶴岡市大山 椙尾神社 文政元年8月
http://www.wasan.jp/yamagata/sugio.html
円 A, B, C, D の半径を r1, r2, r3, r4 とする。円 C, D の中心座標を (x3, y3), (x4, y4) とし，方程式を立てる。
using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive;
cos30d = sqrt(Sym(3))/2;
sin30d = 1//2;
x3 = (r1 + r3)cos30d;
y3 = (r1 + r3)sin30d;
x4 = (1 - r4)cos30d;
y4 = -(1 - r4)sin30d;
eq1 = (x4 + r1*cos30d)^2 + (r1*sin30d-y4)^2 - r4^2;
eq2 = x4^2 + (1 - r2 - y4)^2 - (r2 + r4)^2;
eq3 = (x3 - x4)^2 + (y3 - y4)^2 - (r4 - r3)^2;
この図形においては，r1, r2, r3, r4 は独立に決まるわけではなく，r1 を決めると残りの r2, r3, r4 が決まる。
eq1, eq2, eq3 を r2, r3, r4 について解くと，それぞれは r1 についての式になる。
res = solve([eq1, eq2, eq3], (r2, r3, r4))
1-element Vector{Tuple{Sym, Sym, Sym}}:
(-3*(r1 - 1)/(r1 + 7), -3*r1*(r1 - 1)/(7*r1 + 1), (r1 + 1)/2)
なお，算額では A と B の大きさが同じように見受けられるので，以下（r1 = r2）を解けば対応する図を得る。
solve(-3*(r1 - 1)/(r1 + 7) - r1)[1] |> println
-5 + 2*sqrt(7)
using Plots
function circle(ox, oy, r, color=:red)
θ = 0:0.01:2pi
x = r.*cos.(θ)
y = r.*sin.(θ)
plot!(ox .+ x, oy .+ y, color=color)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:center)
# scatter!([x], [y])
annotate!(x, y, text(string, 10, position, color))
end;
function draw(r1=2*sqrt(7) - 5; more=false)
pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3, r4) = (-3*(r1 - 1)/(r1 + 7), -3*r1*(r1 - 1)/(7*r1 + 1), (r1 + 1)/2)
cos30d = sqrt(Sym(3))/2;
sin30d = 1/2;
x3 = (r1 + r3)cos30d;
y3 = (r1 + r3)sin30d;
x4 = (1 - r4)*cos30d;
y4 = -(1 - r4)*sin30d;
println("$r1, $r2, $r3, $r4")
plot()
if more
point(0, 0, "A", :green)
point(0, 1 - r2, "B", :magenta)
point(x3, y3, "C", :brown)
point(0.8, 0, "D", :blue)
end
circle(0, 0, 1)
circle(0, 0, r1, :green)
circle(0, 1 - r2, r2, :magenta)
circle((1 - r2)cos30d, -(1 - r2)sin30d, r2, :magenta)
circle((-(1 - r2))cos30d, -(1 - r2)sin30d, r2, :magenta)
circle(x3, y3, r3, :brown)
circle(-x3, y3, r3, :brown)
circle(0, -r1-r3, r3, :brown)
circle(x4, y4, r4, :blue)
circle(-x4, y4, r4, :blue)
circle(0, 1-r4, r4, :blue)
end;
r1 = 0.26 のとき
0.26, 0.30578512396694213, 0.20468085106382977, 0.63
r1 = r2 = 2*sqrt(7) - 5 のとき
0.29150262212918143, 0.29150262212918104, 0.20377661238703057, 0.6457513110645907
