算額（その44）

算額（その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