算額（その275）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（262）
長野県飯綱町 牟礼神社 明治31年(1898)
外球に円柱が内接している。甲球 1 個，乙球 2 個，丙球 2 個，丁球 1 個，戊球 1 個がそれぞれが隣の球（甲乙丙丁丙乙甲）の順に外接しており，外球に内接し円柱に外接している。また，外球を含めてすべての球の中心は同一平面上にある。
外球の直径が 85 寸，甲球の直径が 17 寸のとき，乙球の直径を求めよ。
問題は三次元で難しそうであるが，全ての球の中心が同一平面上にあるということから，上方からの透視図を考えると図のような二次元の問題であることがわかる。
それぞれの円の半径，中心座標を以下のようにする。
外円 r0, (0, 0) 既知 r0 = 85 // 2
甲円 r1, (0, r0 - r1) 既知 r1 = 17 // 2
乙円 r2, (x2, y2)
丙円 r3, (x3, y3)
丁円 r4, (0, r4 - r0)
戊円 r5, (0, r5 + 2r4 - r0)
以下の 8 本の方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive,
r4::positive, r5::positive, x2::positive, y2::positive,
x3::positive, y3::positive;
r0 = 85 // 2
r1 = 17 // 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2 |> expand
eq2 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq3 = x2^2 + (y2 - (r5 + 2r4 - r0))^2 - (r2 + r5)^2 |> expand
eq4 = x3^2 + (y3 - (r4 - r0))^2 - (r3 + r4)^2 |> expand
eq5 = x3^2 + (y3 - (r5 + 2r4 - r0))^2 - (r3 + r5)^2 |> expand
eq6 = x2^2 + y2^2 - (r0 - r2)^2 |> expand
eq7 = x3^2 + y3^2 - (r0 - r3)^2 |> expand
eq8 = r1 + r4 + r5 - r0 |> expand;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r2, r3, r4, r5, x2, y2, x3, y3))
1-element Vector{NTuple{8, Sym}}:
(10, 170/11, 85/4, 51/4, 10*sqrt(3), 55/2, 170*sqrt(3)/11, 85/22)
外円 = 42.500; 甲 = 8.500; 乙 = 10.000; 丙 = 15.455; 丁 = 21.250; 戊 = 12.750
x2 = 17.321; y2 = 27.500; x3 = 26.768; y3 = 3.864
すなわち，乙円（乙球）の直径は 20 寸である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r0 = 85 // 2
r1 = 17 // 2
(r2, r3, r4, r5, x2, y2, x3, y3) = res[1]
# (r2, r3, r4, r5, x2, y2, x3, y3) = [19, 20, 21, 22, 20, 20, 20, -20]
@printf("外円 = %.3f; 甲 = %.3f; 乙 = %.3f; 丙 = %.3f; 丁 = %.3f; 戊 = %.3f\n", r0, r1, r2, r3, r4, r5)
@printf("x2 = %.3f; y2 = %.3f; x3 = %.3f; y3 = %.3f\n", x2, y2, x3, y3)
plot()
circle(0, 0, r0, :black)
circle(0, r0 - r1, r1, :green)
circle(x2, y2, r2, :red)
circle(-x2, y2, r2, :red)
circle(x3, y3, r3, :blue)
circle(-x3, y3, r3, :blue)
circle(0, r4 - r0, r4, :magenta)
circle(0, r5 + 2r4 - r0, r5, :orange)
if more
point(r0, 0, "r0 ", :black, :right, :bottom)
point(0, r0 - r1, " 甲:r1\n r0-r1", :green, :left, :bottom)
point(x2, y2, " 乙:r2\n (x2,y2)", :red, :center, :bottom)
point(x3, y3, " 丙:r3\n (x3,y3)", :blue, :center, :bottom)
point(0, r4 - r0, " 丁:r4\n r4-r0", :magenta, :left, :bottom)
point(0, r5 + 2r4 - r0, " 戊:r5\n r5+2r4-r0", :orange, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
