算額（その60）

算額（その60）

石川県松任市若宮 若宮八幡宮 文政13年吉祥日
http://www.wasan.jp/isikawa/wakamiyahatiman1.html

石川県鳳珠郡穴水町 美麻奈比古神社 文政13年8月
http://www.wasan.jp/isikawa/mimanahiko.html

群馬県高崎市 祖師堂 明治16年
http://www.wasan.jp/gunma/sosido.html

外円の中に 8 個の小円と 1 個の大円がある。それぞれの径を求めよ。

外円の半径を 1 とし，図のように記号を定め，方程式を解く。

using SymPy
@syms r1::positive, x2::positive, y2::positive, r2::positive;

x2 = (1 - r2) / sqrt(Sym(2))
y2 = x2
eq1 = r1 + 2r2 - 1 |> expand
eq2 = x2^2 + y2^2 - (1 - r2)^2 |> expand
eq3 = x2^2 + (1 - r2 - y2)^2 - 4r2^2 |> expand;

res = solve([eq1, eq3], (r1, r2))

   2-element Vector{Tuple{Sym, Sym}}:
    (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
    (-4*sqrt(2) + 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 - sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))

解が 2 組得られるが，最初のものが適切な解である。

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;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   (r1, r2) = (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
   x2 = y2 = (1 - r2)/sqrt(2)
   println("r1 = $r1, x2 = $x2, y2 = $y2, r2 = $r2")
   circle(0, 0, 1, :black)
   circle(0, 0, r1, :brown)
   circle(x2, y2, r2)
   circle(x2, -y2, r2)
   circle(-x2, y2, r2)
   circle(-x2, -y2, r2)
   circle(0, 1-r2, r2, :green)
   circle(0, r2-1, r2, :green)
   circle(1-r2, 0, r2, :magenta)
   circle(r2-1, 0, r2, :magenta)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 1-r2, "1-r2 ", :green, :right)
       point(0, r1, "r1 ", :brown, :right)
       point(1-r2, 0, "1-r2 ", :green, :right)
       point(r1, 0, "r1 ", :brown, :right)
       point(x2, y2, "(x2,y2)", :red, :center)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.4464626921716892, x2 = 0.5114017891839755, y2 = 0.5114017891839755, r2 = 0.2767686539141554