算額（その495）
新潟県小千谷市 小千谷二荒神社 天保4年(1833)
涌田和芳，外川一仁:小千谷二荒神社の紛失算額，長岡工業高等専門学校研究紀要，第 51 巻，p.35-40，2015
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_51/51_35wakuta.pdf
2 個の等円が交わってできる領域に，甲円，乙円，丙円，丁円が入っている。甲円，乙円，丙円の直径がそれぞれ 18.2 寸，27．3 寸，9.1 寸のとき，丁円の直径はいかほどか。
等円の半径と中心座標を r0, (0, a), (0, -a)
甲円の半径と中心座標を r1, (x1, 0)
乙円の半径と中心座標を r2, (x2, 0)
丙円の半径と中心座標を r3, (x3, 0)
丁円の半径と中心座標を r4, (x4, 0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive,
r1::positive, x1::positive,
r2::positive, x2::negative,
r3::positive, x3::negative,
r4::positive, x4::negative;
@syms a, r0, r1, r2, r3, r4, x1, x2, x3, x4
x2 = x1 - r1 - r2
x3 = x2 - r2 - r3
x4 = x3 - r3 - r4
eq1 = x1^2 + a^2 - (r0 - r1)^2
eq2 = x2^2 + a^2 - (r0 - r2)^2
eq3 = x3^2 + a^2 - (r0 - r3)^2
eq4 = x4^2 + a^2 - (r0 - r4)^2
res = solve([eq1, eq2, eq3, eq4], (a, r0, x1, r4))
2-element Vector{NTuple{4, Sym}}:
(-2*r2*sqrt(r1*r3*(r1 + r2)*(r2 + r3))/(r1*r3 - r2^2), -r2*(r1 + r2)*(r2 + r3)/(r1*r3 - r2^2), (r1^2*r3 + r1*r2*r3 - r2^3 - r2^2*r3)/(r1*r3 - r2^2), r1*r3^2/(r1*r2 - r1*r3 + r2^2))
(2*r2*sqrt(r1*r3*(r1 + r2)*(r2 + r3))/(r1*r3 - r2^2), -r2*(r1 + r2)*(r2 + r3)/(r1*r3 - r2^2), (r1^2*r3 + r1*r2*r3 - r2^3 - r2^2*r3)/(r1*r3 - r2^2), r1*r3^2/(r1*r2 - r1*r3 + r2^2))
最初のものが適解である。
丁円の半径は r1*r3^2/(r1*r2 - r1*r3 + r2^2) であり，甲円，乙円，丙円の直径がそれぞれ 18.2 寸，27．3 寸，9.1 寸のとき，丁円の直径は 1.4 寸である。
(r1, r2, r3) = (182//20, 273//20, 91//20)
float(2*r1*r3^2//(r1*r2 - r1*r3 + r2^2))
1.4
丁円の直径 = 1.4; r0 = 39; x1 = 16.9; r4 = 0.7
using Plots
function draw(r1, r2, r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r0, x1, r4) = (-2*r2*sqrt(r1*r3*(r1 + r2)*(r2 + r3))/(r1*r3 - r2^2), -r2*(r1 + r2)*(r2 + r3)/(r1*r3 - r2^2), (r1^2*r3 + r1*r2*r3 - r2^3 - r2^2*r3)/(r1*r3 - r2^2), r1*r3^2/(r1*r2 - r1*r3 + r2^2))
x2 = x1 - r1 - r2
x3 = x2 - r2 - r3
x4 = x3 - r3 - r4
@printf("丁円の直径 = %g; r0 = %g; x1 = %g; r4 = %g\n", 2r4, r0, x1, r4)
x = sqrt(r0^2 - a^2)
θ = atand(a/x)
plot()
circle(0, a, r0, :red, beginangle=180+θ, endangle=360-θ)
circle(0, -a, r0, :red, beginangle=θ, endangle=180-θ)
circle(x1, 0, r1, :green)
circle(x2, 0, r2, :orange)
circle(x3, 0, r3, :magenta)
circle(x4, 0, r4, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(x1, 0, "甲円:r1\nx1", :green, :center, :bottom, delta=delta)
point(x2, 0, "乙円:r2\nx2", :orange, :center, :bottom, delta=delta)
point(x3, 0, "丙円:r3\nx3", :magenta, :center, :bottom, delta=delta)
point(x4, 0, " 丁円:r4\nx4", :blue, :center, :top, delta=-5delta)
point(0, a, " a", :red, :left, :vcenter)
point(0, r0 - a, " r0-a", :red, :left, :bottom, delta=delta)
else
plot!(showaxis=false)
end
end;
