算額(その44)
山形県鶴岡市羽黒町 出羽三山神社
http://www.wasan.jp/yamagata/haguro.html
複数の小さな円が大きな円と直線に接している。
以下のように記号を定め,小円の半径と中心座標の漸化式を求める。
using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive;
@syms x1::positive, x2::positive, x3::positive;
eq1 = (x3 - x2)^2 + (r3 - r2)^2 - (r3 + r2)^2
eq2 = x2^2 + (R + 2r1 - r2)^2 - (R + r2)^2
eq3 = x3^2 + (R + 2r1 - r3)^2 - (R + r3)^2
res = solve([eq1, eq2, eq3], (r3, x3))
2-element Vector{Tuple{Sym, Sym}}:
(-(-4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2)))/(4*(R + r1 - r2)), -sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2))
(-(-4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2)))/(4*(R + r1 - r2)), sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2))
大円の半径を 25,最小の円の半径を 1 として図を描く。
using Plots
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=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
nextcircle(r2, x2, R, r1) = (-(-4*R*r1 - 4*r1^2 + x2^2 - 2*x2*(-sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2)))/(4*(R + r1 - r2)),
-sqrt(r2)*sqrt(R + r1)*sqrt(4*R*r1 + 4*r1^2 - 4*r1*r2 + x2^2)/(-R - r1 + r2) + x2*(R + r1)/(R + r1 - r2));
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = 25
r1 = 1
x1 = 0
plot()
circle(0, R + 2r1, R)
circle(0, r1, r1, :magenta)
ri, xi = r1, x1
println("x1 = $xi, r1 = $ri")
for i in 2:7
ri, xi = nextcircle(ri, xi, R, r1)
println("x$i = $xi, r$i = $ri")
circle(xi, ri, ri, :blue)
circle(-xi, ri, ri, :blue)
more && i >= 6 && point(xi, ri, "(x$i,r$i)", :blue, :center)
end
if more
point(0, R+2r1, "R+2r1 ", :red, :right)
vline!([0], color=:black, lw=0.25)
end
hline!([0], color=:black, lw=0.25)
end;
x1 = 0, r1 = 1
x2 = 2.039607805437114, r2 = 1.04
x3 = 4.249182927993988, r3 = 1.1736111111111112
x4 = 6.860498981924838, r4 = 1.4525619834710746
x5 = 10.283736834136707, r5 = 2.0168773391709625
x6 = 15.436610653781944, r6 = 3.291239889196675
x7 = 25.06402103903699, r7 = 7.040434140820087