算額（その135）

算額（その135）

岐阜県大垣市西外側町 大垣八幡神社 天保年間
http://www.wasan.jp/gifu/ogakihatiman.html

外円に等弧を描き，天円 3 個，地円 1 個，人円 2 個が入っている。天径が与えられたとき，人径を求めよ。

他の問題に準じ，外円の半径を 1 とする。図のように記号を定め方程式を解く。

外円の中心座標と半径　　　　(0, 0), 1
右の天円の中心座標と半径　　(x1, y1), r1
地円の中心座標と半径　　　　(0, r2-1), r2
右の人円の中心座標と半径　　(x3, y3), r3
円の一部が等弧になる円の半径は外円と同じ 1 とし，中心座標を (√3/3, 1), (-√3/3, 1) とする。
この 2 つの円と外円の中心は，正三角形を構成する（正三角形である必要性はない。中心座標が異なっても，それぞれの円が互いに接することには変わりない）。

using SymPy

@syms x0, y0, r0, x1::positive, y1::positive, r1::positive, r2::positive, x3::positive, y3::negative, r3::positive;

y0 = 1
x0 = y0 / sqrt(Sym(3))
eq1 = x0^2 + (y0 - 1 + r1)^2 - (1 - r1)^2
eq2 = x0^2 + (y0 - r2 + 1)^2 - (1 + r2)^2
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - (1 - r1)^2
eq4 = (x0 - x3)^2 + (y0 - y3)^2 - (1 + r3)^2
eq5 = x3^2 + (y3 - r2 + 1)^2 - (r2 + r3)^2
eq6 = (x1 + x0)^2 + (y0 - y1)^2 - (1 + r1)^2
eq7 = x1^2 + y1^2 - (1 - r1)^2
eq8 = x3^2 + y3^2 - (1 - r3)^2;

res = solve([eq1, eq2, eq4, eq5, eq6, eq7, eq8], (r1, x1, y1, r2, r3, x3, y3))

   1-element Vector{NTuple{7, Sym}}:
    (1/3, sqrt(3)/3, 1/3, 5/9, 5/17 - 5*sqrt(2)/102, 5*sqrt(3)/102 + 55*sqrt(6)/204, 1/34 - 35*sqrt(2)/204)

天円　半径 1/3, 中心 (√3/3, 1/3)
地円　半径 5/9
人円　半径 (60 - 10*√2)/204, 中心 ((10√3 + 55*√6)/204, (6 - 35*√2)/204)

人円の直径は，天円の径の (60 - 10*√2)/68 倍である。

@syms x, y
eq11 = x^2 + y^2 - 1
eq12 = (x - x0)^2 + (y - y0)^2 - 1
res2 = solve([eq11, eq12], (x, y))
(xs, ys) = res2[1]
(xe, ye) = res2[2]
θs = float((asin(ys-y0) * 180 / pi).evalf())
θe = float((asin(ye-y0) * 180 / pi).evalf());

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 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(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, y1, r2, r3, x3, y3) = res[1] # [1.949028691169056, 0.6774004450055521, 0.5770838618531914, 0.11011396993484983, 0.8155454510489388, 0.8056721998501019, 0.3178625307754063, 0.1338912989224465]
   r0 = 1
   y0 = 1
   x0 = y0 / sqrt(Sym(3))
   @printf("天円の半径 = %.5f,  地円の半径 = %.5f,  人円の半径 = %.5f\n", r1, r2, r3)
   plot()
   circle(0, 0, 1, :green)
   circle(0, r2 - 1, r2, :blue)
   circle(0, 1 - r1, r1)
   circle(x1, y1, r1)
   circle(-x1, y1, r1)
   circle(x3, y3, r3, :orange)
   circle(-x3, y3, r3, :orange)
   if more == true
       @printf("r1 = %.5f;  x1 = %.5f;  y1 = %.5f\nr2 = %.5f\nr3 = %.5f;  x3 = %.5f;  y3 = %.5f\n\n", r1, x1, y1, r2, r3, x3, y3) 
       #thetas = 180 - round(Int, float(θs))
       #thetae = 360 + round(Int, float(θe))
       circle(x0, y0, r0, :magenta, lw=0.5, linestyle=:dash)
       #thetas = 180 - round(Int, float(θe))
       #thetae = 360 + round(Int, float(θs))
       circle(-x0, y0, r0, :magenta, lw=0.5, linestyle=:dash)
       point(x0, y0, " (x0,y0;1)", :magenta)
       point(-x0, y0, "(-x0,y0;1) ", :magenta, :right)
       point(0, 1 - r1, "天:1-r1 ", :red, :right)
       point(0, r2 - 1, "地:r2-1 ", :blue, :right)
       point(x1, y1, "(x1,y1;r1)", :red, :top)
       point(x3, y3, "人:(x3,y3;r3)", :orange, :top)
       point(0, 0, "", :green)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
   thetas = 180 - round(Int, θs)
   thetae = 360 + round(Int, θe)
   circle(x0, y0, r0, :green, beginangle=thetas, endangle=thetae, lw=0.5)
   thetas = 180 - round(Int, θe)
   thetae = 360 + round(Int, θs)
   circle(-x0, y0, r0, :green, beginangle=thetas, endangle=thetae, lw=0.5)
end;

   天円の半径 = 0.33333,  地円の半径 = 0.55556,  人円の半径 = 0.22479
   r1 = 0.33333;  x1 = 0.57735;  y1 = 0.33333
   r2 = 0.55556
   r3 = 0.22479;  x3 = 0.74531;  y3 = -0.21322