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