算額（その140）

算額（その140）

福島県白河市明神 境明神 万延元年(1860)9月
http://www.wasan.jp/fukusima/sakai1.html

長方形の中に，半円 2 個，天円 4 個，地円 2 個，人円 2 個が入っている。人円の径が 1 寸のとき，天円の径はいくらか。

半円，天円，地円，人円の半径 をそれぞれ r0, r1, r2, r3 = 1 とおく。天円，地円の中心の x 座標を x1, x2 とおく。方程式を立て，解く。

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, x1::positive, x2::positive;
r3 = 1
eq1 = x2^2 + r3^2 - (r2 + r3)^2  # 地円と人円が外接する
eq2 = x2^2 + (2r3 + r0)^2 - (r0 + r2)^2  # 半円と地円が外接する
eq3 = (x1 - x2)^2 + r1^2 - (r1 + r2)^2  # 天円と地円が外接する
eq4 = x1^2 + (2r3 + r0 - r1)^2 - (r0 + r1)^2  # 半円と天円が外接する
eq5 = 2r3 + r0 - 2r1;  # 半径の関係

res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, x1, x2))

   1-element Vector{NTuple{5, Sym}}:
    (6, 4, 14/5, 2*sqrt(21), 4*sqrt(21)/5)

天円の半径は 4，元の単位では天円の直径は 4 寸。
ちなみに半円の直径は 6 寸，地円の直径は 14/5 = 2寸8分。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r2, x1, x2) = res[1]
   @printf("半円の半径 = %.3f\n", r0)
   @printf("天円の半径 = %.3f\n", r1)
   @printf("地円の半径 = %.3f\n", r2)
   plot([x1 + r1, x1 + r1, -x1 - r1, -x1 - r1, x1 + r1], [-2r1, 2r1, 2r1, -2r1, -2r1])
   circle(0, 2r3 + r0, r0, beginangle=180, endangle=360)
   circle(0, -2r3 - r0, r0, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle(x1, -r1, r1, :blue)
   circle(-x1, r1, r1, :blue)
   circle(-x1, -r1, r1, :blue)
   circle(x2, 0, r2, :green)
   circle(-x2, 0, r2, :green)
   circle(0, r3, r3, :orange)
   circle(0, -r3, r3, :orange)
   if more == true
       point(0, 2r3 + r0, "2r3+r0 ", :red, :right)
       point(0, r3, "r3", :orange, :right)
       point(x2, 0, "x2", :green, :center)
       point(x2+r1, 0, "x2+r2 ", :green, :right)
       point(x1, 0, "x1", :blue, :center)
       point(x1, r1, "(x1,r1)", :blue, :center)
       point(x1 + r1, 0, "x1+r1", :blue, :right)
       point(x1 + r1, 2r1, "(x1+r1,2r1) ", :blue, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

draw(false)

   半円の半径 = 6.000
   天円の半径 = 4.000
   地円の半径 = 2.800