算額（その227）

算額（その227）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（217）長野県下高井郡山ノ内町 渋温泉医王殿 慶応元年(1865)

正方形の中に，四分円，大円，中円が各1個，小円が2個入っている。正方形の一辺の長さが与えられたとき，大円，中円，小円の径を求めよ。

円弧を一部とする円の半径，中心座標を r0, (r0, 0)
大円の半径，中心座標を r1，(r0 - r1, r1)
中円の半径，中心座標を r2，(r2, r0 - r2)
小円の半径，中心座標を r3，(r3, y3)
とおき，以下の方程式を r1, r2, r3, y3 について解く（それぞれの解は r0 を含む）。

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, y3::positive;

eq1 = 2(r0 - r2)^2 - (r0 + r2)^2 |> expand
eq2 = (r0 - r3)^2 + y3^2 - (r0 + r3)^2 |> expand
eq3 = (r2 - r3)^2 + (r0 - r2 - y3)^2 - (r2 + r3)^2 |> expand
eq4 = 2r1^2 - (r0 - r1)^2

res = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, y3))

   3-element Vector{NTuple{4, Sym}}:
    (r0*(-1 + sqrt(2)), r0*(3 - 2*sqrt(2)), r0/2, sqrt(2)*r0)
    (r0*(-1 + sqrt(2)), r0*(3 - 2*sqrt(2)), r0*(3 - 2*sqrt(2))/2, r0*(2 - sqrt(2)))
    (r0*(-1 + sqrt(2)), r0*(2*sqrt(2) + 3), r0*(2*sqrt(2) + 3)/2, r0*(sqrt(2) + 2))

2 番目の解が適切な解である。

r1 = r0*(√2 - 1)
r2 = r0*(3 - 2√2)
r3 = r0*(3 - 2√2)/2 =  r2 / 2
y3 = r0*(2 - √2)

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(r0, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, y3) = (r0*(-1 + sqrt(2)), r0*(3 - 2*sqrt(2)), r0*(3 - 2*sqrt(2))/2, r0*(2 - sqrt(2)))
   r1 = r0*(√2 - 1)
   r2 = r0*(3 - 2√2)
   r3 = r0*(3 - 2√2)/2  # =  r2 / 2
   y3 = r0*(2 - √2)
   @printf("r1 = %.7f,  r2 = %.7f;  r3 = %.7f; y3 = %.7f\n", r1, r2, r3, y3)
   println("r3/r2 = $(r3/r2)")
   plot([0, r0, r0, 0, 0], [0, 0, r0, r0, 0], color=:black, lw=0.5)
   circle(r0, 0, r0, beginangle=90, endangle=180)
   circle(r0 - r1, r1, r1, :blue)
   circle(r2, r0 - r2, r2, :green)
   circle(r3, y3, r3, :orange)
   circle(r0 - y3, r0 - r3, r3, :orange)
   if more == true
       point(r0, 0, "r0=$r0 ", :black, :right, :bottom)
       point(r0 - r1, r1, "(r0-r1,r1)", :blue, :center)
       point(r2, r0 - r2, "(r2,r0-r2)", :green, :center)
       point(r3, y3, "(r3,y3)", :orange, :center)
       point(r0 - y3, r0 - r3, "(ro-y3,r0-r3)", :orange, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;